Catch potential DIC calculation issue

5 files changed, 25 insertions, 19 deletions

diff --git a/R/autojags.R b/R/autojags.R
index 27dd870..15b0e9d 100644
--- a/R/autojags.R
+++ b/R/autojags.R
@@ -131,7 +131,7 @@ autojags <- function(data,inits=NULL,parameters.to.save,model.file,n.chains,n.ad
   samples <- order_samples(samples, parameters.to.save)
   
   #Convert rjags output to jagsUI form 
-  output <- process_output(samples, coda_only = codaOnly, quiet = !verbose)
+  output <- process_output(samples, coda_only = codaOnly, DIC, quiet = !verbose)
   if(is.null(output)){
     output <- list()
     output$samples <- samples
diff --git a/R/jags.R b/R/jags.R
index 7e58b58..9523575 100644
--- a/R/jags.R
+++ b/R/jags.R
@@ -64,7 +64,7 @@ jagsUI <- jags <- function(data,inits=NULL,parameters.to.save,model.file,n.chain
   samples <- order_samples(samples, parameters.to.save)
   
   #Convert rjags output to jagsUI form 
-  output <- process_output(samples, coda_only = codaOnly, quiet = !verbose)
+  output <- process_output(samples, coda_only = codaOnly, DIC, quiet = !verbose)
   if(is.null(output)){
     output <- list()
     output$samples <- samples
diff --git a/R/process_output.R b/R/process_output.R
index 4877b8b..72bb669 100644
--- a/R/process_output.R
+++ b/R/process_output.R
@@ -2,7 +2,7 @@
 #Process output master function
 #To generate backwards-compatible jagsUI output
 #WIP
-process_output <- function(mcmc_list, coda_only=NULL, quiet=FALSE){
+process_output <- function(mcmc_list, coda_only=NULL, DIC, quiet=FALSE){
   if(!quiet){cat('Calculating statistics.......','\n')} 
 
   tryCatch({
@@ -15,7 +15,7 @@ process_output <- function(mcmc_list, coda_only=NULL, quiet=FALSE){
     # Get final summary table
     sum_list <- list(summary = stat_summary_table(stats, coda_only))
     # DIC stuff
-    dic_list <- calc_DIC(mcmc_list)
+    dic_list <- calc_DIC(mcmc_list, DIC)
     
     # Bind it all together
     if(!quiet){cat('\nDone.','\n')}
@@ -149,6 +149,7 @@ calc_f <- function(values, mn){
 
 calc_Rhat <- function(mcmc_list){
   stopifnot(has_one_parameter(mcmc_list))
+  if(length(mcmc_list) == 1) return(NA)
   out <- try(coda::gelman.diag(mcmc_list, 
             autoburnin=FALSE, multivariate=FALSE)$psrf[1])
   if(inherits(out, "try-error") || !is.finite(out)) out <- NA
@@ -241,10 +242,10 @@ calc_stats <- function(mcmc_list, coda_only=NULL){
 
 #------------------------------------------------------------------------------
 #Calculate pD and DIC from deviance if it exists in output samples
-calc_DIC <- function(samples){ 
-  ind <- which_params('deviance', param_names(samples))
-  #if(is.null(ind)) return(c(pD=NA, DIC=NA))
-  if(is.null(ind)) return(NULL)
+calc_DIC <- function(samples, DIC){ 
+  if(!DIC | !("deviance" %in% param_names(samples))){
+    return(NULL)
+  }
 
   dev <- mcmc_to_mat(samples[,'deviance'])  
   #if(any(is.na(dev)) || any(is.infinite(dev))) return(c(pD=NA, DIC=NA))
diff --git a/R/update.R b/R/update.R
index ad65d27..f1d2db5 100644
--- a/R/update.R
+++ b/R/update.R
@@ -52,7 +52,7 @@ update.jagsUI <- function(object, parameters.to.save=NULL, n.adapt=NULL, n.iter,
   samples <- order_samples(samples, parameters)
     
   #Run process output
-  output <- process_output(samples, coda_only = codaOnly, quiet = !verbose)
+  output <- process_output(samples, coda_only = codaOnly, DIC, quiet = !verbose)
   if(is.null(output)){
     output <- list()
     output$samples <- samples
diff --git a/inst/tinytest/test_process_output.R b/inst/tinytest/test_process_output.R
index fae78d5..84ab2b2 100644
--- a/inst/tinytest/test_process_output.R
+++ b/inst/tinytest/test_process_output.R
@@ -16,7 +16,7 @@ param_names <- jagsUI:::param_names
 
 # test that process_output generates correct list of output--------------------
 samples <- readRDS('coda_samples.Rds')
-out <- process_output(samples, quiet=TRUE)
+out <- process_output(samples, DIC=TRUE, quiet=TRUE)
 expect_inherits(out,'list')
 expect_equal(length(out),15)
 expect_equal(names(out),c("sims.list", "mean", "sd", "q2.5", "q25", "q50", 
@@ -31,16 +31,19 @@ colnames(cs)[3:7] <- c("2.5%","25%","50%","75%","97.5%")
 cs <- cs[,c(1:7,10,11,8,9)]
 expect_equal(out$summary, cs)
 #Check error handling
-expect_message(result <- process_output("junk", quiet=TRUE))
+expect_message(result <- process_output("junk", DIC=TRUE, quiet=TRUE))
 expect_true(all(is.na(result)))
 
 #DIC=FALSE
-out2 <- process_output(samples[,-coda::nvar(samples)], quiet=TRUE)
+out2 <- process_output(samples[,-coda::nvar(samples)], DIC=FALSE, quiet=TRUE)
+expect_true(is.null(out2$DIC))
+expect_true(is.null(out2$pD))
+out2 <- process_output(samples[,-coda::nvar(samples)], DIC=TRUE, quiet=TRUE)
 expect_true(is.null(out2$DIC))
 expect_true(is.null(out2$pD))
 
 #Exclude parameters
-out3 <- process_output(samples, coda_only=c("alpha","kappa", "mu"),quiet=TRUE)
+out3 <- process_output(samples, coda_only=c("alpha","kappa", "mu"), DIC=TRUE, quiet=TRUE)
 expect_identical(names(out3$sims.list), names(out$sims.list))
 expect_identical(names(out3$mean), names(out$mean))
 expect_false(any(is.na(unlist(out3$mean))))
@@ -48,7 +51,7 @@ expect_identical(rownames(out3$summary), c("beta", "sigma", "deviance"))
 
 #test that process_output matches old jagsUI process.output--------------------
 old_all <- readRDS("old_jagsUI_output.Rds")
-new_po <- process_output(old_all$samples, quiet=TRUE)
+new_po <- process_output(old_all$samples, DIC=TRUE, quiet=TRUE)
 old_po <- readRDS("old_process_output.Rds")  
 expect_identical(new_po$summary, old_all$summary)
 new_po$summary <- NULL
@@ -227,16 +230,18 @@ expect_identical(ref_output, st_sub)
 
 # test that calculation of pD/DIC works----------------------------------------
 samples <- readRDS('coda_samples.Rds')
-expect_equal(calc_DIC(samples), c(pD=6.660906,DIC=40.712014), tol=1e-4)
+expect_equal(calc_DIC(samples, DIC=TRUE), c(pD=6.660906,DIC=40.712014), tol=1e-4)
 dev_ind <- which_params('deviance', param_names(samples))
 no_dev <- samples[,-coda::nvar(samples)]
-expect_true(is.null(calc_DIC(no_dev)))
+expect_true(is.null(calc_DIC(no_dev, DIC=TRUE)))
+expect_true(is.null(calc_DIC(no_dev, DIC=FALSE)))
+expect_true(is.null(calc_DIC(samples, DIC=FALSE)))
 samp_na <- samples
 ind <- which_params('deviance',param_names(samples))
 samp_na[[1]][1,ind] <- NA 
-expect_true(is.null(calc_DIC(samp_na)))
+expect_true(is.null(calc_DIC(samp_na, DIC=TRUE)))
 samp_inf <- samples
 samp_inf[[1]][1,ind] <- Inf
-expect_true(all(is.na(calc_DIC(samp_na))))
+expect_true(all(is.na(calc_DIC(samp_na, DIC=TRUE))))
 samp_inf[[1]][1,ind] <- -Inf
-expect_true(is.null(calc_DIC(samp_na)))
+expect_true(is.null(calc_DIC(samp_na, DIC=TRUE)))