Use older Rds format, use process_input in jags

19 files changed, 548 insertions, 4 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 9e3a7b4..75d1e6c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -16,6 +16,7 @@ Imports:
     grDevices,
     graphics,
     utils
+Suggests: tinytest
 SystemRequirements: JAGS (http://mcmc-jags.sourceforge.net)
 Description: A set of wrappers around 'rjags' functions to run Bayesian analyses in 'JAGS' (specifically, via 'libjags').  A single function call can control adaptive, burn-in, and sampling MCMC phases, with MCMC chains run in sequence or in parallel. Posterior distributions are automatically summarized (with the ability to exclude some monitored nodes if desired) and functions are available to generate figures based on the posteriors (e.g., predictive check plots, traceplots). Function inputs, argument syntax, and output format are nearly identical to the 'R2WinBUGS'/'R2OpenBUGS' packages to allow easy switching between MCMC samplers.
 License: GPL-3
diff --git a/R/jags.R b/R/jags.R
index 3436e0e..33d5e2f 100644
--- a/R/jags.R
+++ b/R/jags.R
@@ -66,7 +66,8 @@ jagsUI <- jags <- function(data,inits=NULL,parameters.to.save,model.file,n.chain
   }
   
   #Convert rjags output to jagsUI form 
-  output <- process.output(samples,DIC=DIC,codaOnly,verbose=verbose)
+  #output <- process.output(samples,DIC=DIC,codaOnly,verbose=verbose)
+  output <- process_output(samples, coda_only = codaOnly, quiet = !verbose)
   if(is.null(output)){
     output <- list()
     samples <- order.params(samples,parameters.to.save,DIC,verbose=verbose)
@@ -80,7 +81,7 @@ jagsUI <- jags <- function(data,inits=NULL,parameters.to.save,model.file,n.chain
   #Add additional information to output list
   
   #Summary
-  output$summary <- summary.matrix(output,samples,n.chains,codaOnly)
+  #output$summary <- summary.matrix(output,samples,n.chains,codaOnly)
  
   output$samples <- samples
   output$modfile <- model.file
diff --git a/R/mcmc_tools.R b/R/mcmc_tools.R
index 24c565c..01813fd 100644
--- a/R/mcmc_tools.R
+++ b/R/mcmc_tools.R
@@ -70,3 +70,10 @@ get_inds <- function(param, params_raw){
   inds <- as.integer(unlist(inds_raw))
   matrix(inds, byrow=T, ncol=length(inds_raw[[1]]))
 }
+
+
+#------------------------------------------------------------------------------
+# Check if mcmc.list has only one parameter (one column)
+has_one_parameter <- function(mcmc_list){
+  coda::nvar(mcmc_list) == 1
+}
diff --git a/R/process_output.R b/R/process_output.R
new file mode 100644
index 0000000..d468f93
--- /dev/null
+++ b/R/process_output.R
@@ -0,0 +1,254 @@
+#------------------------------------------------------------------------------
+#Process output master function
+#To generate backwards-compatible jagsUI output
+#WIP
+process_output <- function(mcmc_list, coda_only=NULL, quiet=FALSE){
+  if(!quiet){cat('Calculating statistics.......','\n')} 
+
+  tryCatch({
+    # Get the sims.list
+    sims <- list(sims.list = sims_list(mcmc_list))
+    # Calculate all stats
+    stats <- calc_stats(mcmc_list, coda_only)
+    # Convert them into stat arrays
+    stats_list <- all_stat_arrays(stats, coda_only)
+    # Get final summary table
+    sum_list <- list(summary = stat_summary_table(stats, coda_only))
+    # DIC stuff
+    dic_list <- calc_DIC(mcmc_list)
+    
+    # Bind it all together
+    if(!quiet){cat('\nDone.','\n')}
+    c(sims, stats_list, dic_list, sum_list)
+  }, error = function(e) {
+    message(paste0("Processing output failed with this error:\n",e,"\n"))
+    list(sims.list=NA, pD=NA, DIC=NA, summary=NA)
+  })
+}
+
+
+#------------------------------------------------------------------------------
+#Fill an array from vector using matching array indices
+fill_array <- function(data_vector, indices){
+  out <- array(NA, dim=apply(indices,2,max))
+  out[indices] <- data_vector
+  out
+}
+
+
+#------------------------------------------------------------------------------
+#Extract the posterior of a parameter and organize it into an array
+get_posterior_array <- function(parameter, samples){
+  
+  tryCatch({
+    #Subset output columns matching parameter
+    col_inds <- which_params(parameter, param_names(samples))
+    posterior_raw <- do.call(rbind, samples[,col_inds,drop=FALSE])
+  
+    #If parameter is scalar, return it now
+    if( ncol(posterior_raw) == 1 ){ return(as.vector(posterior_raw)) }
+
+    #If parameter is array, get indices
+    ind_raw <- get_inds(parameter, colnames(posterior_raw))
+    ndraws <- nrow(posterior_raw)
+    ind_array <- cbind(1:ndraws, ind_raw[rep(1:nrow(ind_raw), each=ndraws),])
+
+    #Create, fill, return output object
+    fill_array(as.vector(posterior_raw), ind_array)
+  }, error = function(e) {
+    message(paste0("Caught error when creating sims.list array for '",
+                   parameter,"':\n",e,"\n"))
+    NA
+  })
+}
+
+
+#------------------------------------------------------------------------------
+#Get sims list
+sims_list <- function(samples){
+  params <- param_names(samples)
+  sapply(strip_params(params, unique=TRUE), get_posterior_array, 
+                      samples, simplify=FALSE) 
+}
+
+
+#------------------------------------------------------------------------------
+#Extract stats for a parameter and organize into appropriately-sized array
+get_stat_array <- function(parameter, stat, model_summary){
+  
+  tryCatch({
+    #Subset vector of stats for parameter
+    row_ind <- which_params(parameter, rownames(model_summary))
+    stat_vector <- model_summary[row_ind, stat]
+
+    #If parameter is scalar, return it now
+    if( length(stat_vector) == 1 ) return(stat_vector)
+
+    #If parameter is array, get indices
+    ind_array <- get_inds(parameter, names(stat_vector))
+
+    #Create, fill, return output object
+    fill_array(stat_vector, ind_array)
+  }, error = function(e) {
+    message(paste0("Caught error when creating stat array for '",
+                   parameter,"':\n",e,"\n"))
+    NA
+  })
+}
+
+
+#------------------------------------------------------------------------------
+#Compile all stats for all parameters into list of lists
+all_stat_arrays <- function(summary_stats, coda_only){
+  
+  stat_array_list <- function(stat, summary_stats){
+    params <- strip_params(rownames(summary_stats), unique=TRUE)
+    sapply(params, function(x){
+      # If the parameter is in coda_only and the stat is not the mean, return NA
+      if(x %in% coda_only & stat != "mean") return(NA)
+      # Otherwise return the stat array for that parameter and stat
+      get_stat_array(x, stat, summary_stats)
+    }, simplify=FALSE)
+  }
+  # Do this for all stats
+  out <- sapply(colnames(summary_stats), stat_array_list, summary_stats,
+                simplify=FALSE)
+  
+  # Convert overlap0 to logical to match old jagsUI code
+  out$overlap0 <- lapply(out$overlap0, function(x) x == 1)
+  out
+}
+
+
+#------------------------------------------------------------------------------
+# Convert stats into summary table in original jagsUI format
+# For backwards compatibility
+stat_summary_table <- function(stats, coda_only){
+  # Move overlap 0 and f to the end of the table
+  stats <- stats[,c("mean", "sd", "q2.5", "q25", "q50", "q75", "q97.5",
+                    "Rhat", "n.eff", "overlap0", "f")]
+  # Rename the quantile columns
+  colnames(stats)[3:7] <- c("2.5%", "25%", "50%", "75%", "97.5%")
+  # Remove rows marked as coda_only
+  keep_rows <- ! strip_params(rownames(stats)) %in% coda_only
+  stats[keep_rows,,drop=FALSE]
+}
+
+
+#------------------------------------------------------------------------------
+#Determine if 95% credible interval of parameter overlaps 0
+overlap_0 <- function(lower, upper){
+  as.numeric(!(lower <= 0) == (upper < 0))
+}
+
+#Calculate proportion of posterior with same sign as mean
+calc_f <- function(values, mn){
+  if(mn >= 0) return(mean(values>=0,na.rm=TRUE))
+  mean(values<0, na.rm=TRUE)
+}
+
+calc_Rhat <- function(mcmc_list){
+  stopifnot(has_one_parameter(mcmc_list))
+  coda::gelman.diag(mcmc_list)$psrf[1]
+}
+
+mcmc_to_mat <- function(mcmc_list){
+  stopifnot(has_one_parameter(mcmc_list))
+  matrix(unlist(mcmc_list), 
+         nrow=coda::niter(mcmc_list), ncol=coda::nchain(mcmc_list))
+}
+
+# Based on R2WinBUGS code
+calc_neff <- function(mcmc_list){
+  niter <- coda::niter(mcmc_list)
+  nchain <- coda::nchain(mcmc_list)
+  mcmc_mat <- mcmc_to_mat(mcmc_list)
+  
+  xdot <- apply(mcmc_mat, 2, mean, na.rm=TRUE)
+  s2 <- apply(mcmc_mat, 2, var, na.rm=TRUE)
+  W <- mean(s2)
+
+  #Non-degenerate case
+  if(is.na(W)){
+	  n_eff <- NA
+  } else if ((W > 1.e-8) && (nchain > 1)) {
+    B <- niter * var(xdot)
+    sig2hat <- ((niter-1)*W + B)/ niter
+    n_eff <- round(nchain * niter * min(sig2hat/B,1),0)
+  } else {
+    #Degenerate case
+    n_eff <- 1
+  }
+  n_eff
+}
+
+#Calculate series of statistics for one parameter
+#Takes an mcmc.list as input
+calc_param_stats <- function(mcmc_list, coda_only){
+  stopifnot(has_one_parameter(mcmc_list))
+  values <- unlist(mcmc_list)
+  stat_names <- c('mean','sd','q2.5','q25','q50','q75','q97.5',
+                  'overlap0','f','Rhat','n.eff')
+
+  fallback <- sapply(stat_names, function(x) NA)
+  if(any(is.infinite(values)) | all(is.na(values))){
+    return(fallback)
+  }
+
+  #Handle any unexpected errors during calculation
+  tryCatch({
+    # If the parameter is in codaOnly, return only the mean
+    mn <- mean(values, na.rm=TRUE)
+    if(coda_only){
+      fallback['mean'] <- mn
+      return(fallback)
+    }
+    # Otherwise calculate all stats
+    quants <- stats::quantile(values, c(0.025, 0.25, 0.5, 0.75, 0.975), na.rm=TRUE)
+    out <- c(mn, 
+             stats::sd(values,na.rm=TRUE), 
+             quants, 
+             overlap0 = overlap_0(quants[1], quants[5]),
+             calc_f(values, mn),
+             calc_Rhat(mcmc_list), 
+             calc_neff(mcmc_list))
+    names(out) <- stat_names
+    out
+  }, error = function(e) {
+    message(paste0('Caught error when calculating stats:\n',e,'\n'))
+    fallback
+  })
+}
+
+
+#------------------------------------------------------------------------------
+#Calculate statistics for all parameters in posterior and organize into matrix
+#Takes mcmc.list as input
+calc_stats <- function(mcmc_list, coda_only=NULL){
+  params <- param_names(mcmc_list)
+  coda_only <- strip_params(params) %in% coda_only 
+
+  out <- sapply(1:length(params), function(i){
+    calc_param_stats(mcmc_list[,i], coda_only[i])
+  })
+  colnames(out) <- params
+  t(out)
+}
+
+
+#------------------------------------------------------------------------------
+#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)
+
+  dev <- mcmc_to_mat(samples[,'deviance'])  
+  #if(any(is.na(dev)) || any(is.infinite(dev))) return(c(pD=NA, DIC=NA))
+  if(any(is.na(dev)) || any(is.infinite(dev))) return(NULL)
+
+  pd <- apply(dev,2,FUN=function(x) stats::var(x)/2)
+  dic <- apply(dev,2,mean) + pd
+  
+  c(pD=mean(pd),DIC=mean(dic))
+}
diff --git a/inst/tinytest/calc_stats_output.Rds b/inst/tinytest/calc_stats_output.Rds
new file mode 100644
index 0000000..263f9fa
--- /dev/null
+++ b/inst/tinytest/calc_stats_output.Rds
diff --git a/inst/tinytest/coda_samples.Rds b/inst/tinytest/coda_samples.Rds
index db70638..3d25cc5 100644
--- a/inst/tinytest/coda_samples.Rds
+++ b/inst/tinytest/coda_samples.Rds
diff --git a/inst/tinytest/jagsbasic_ref_saved.Rds b/inst/tinytest/jagsbasic_ref_saved.Rds
new file mode 100644
index 0000000..1aa7e55
--- /dev/null
+++ b/inst/tinytest/jagsbasic_ref_saved.Rds
diff --git a/inst/tinytest/jagsbasic_reference_fit.Rds b/inst/tinytest/jagsbasic_reference_fit.Rds
new file mode 100644
index 0000000..a430aaa
--- /dev/null
+++ b/inst/tinytest/jagsbasic_reference_fit.Rds
diff --git a/inst/tinytest/longley_reference_fit.Rds b/inst/tinytest/longley_reference_fit.Rds
index 328069b..6947a6e 100644
--- a/inst/tinytest/longley_reference_fit.Rds
+++ b/inst/tinytest/longley_reference_fit.Rds
diff --git a/inst/tinytest/old_jagsUI_output.Rds b/inst/tinytest/old_jagsUI_output.Rds
new file mode 100644
index 0000000..8982dc1
--- /dev/null
+++ b/inst/tinytest/old_jagsUI_output.Rds
diff --git a/inst/tinytest/old_process_output.Rds b/inst/tinytest/old_process_output.Rds
new file mode 100644
index 0000000..1060be9
--- /dev/null
+++ b/inst/tinytest/old_process_output.Rds
diff --git a/inst/tinytest/one_sample.Rds b/inst/tinytest/one_sample.Rds
index 5274d8c..ff6bc92 100644
--- a/inst/tinytest/one_sample.Rds
+++ b/inst/tinytest/one_sample.Rds
diff --git a/inst/tinytest/reference_codaOnly.Rds b/inst/tinytest/reference_codaOnly.Rds
index 23c06ed..4081a59 100644
--- a/inst/tinytest/reference_codaOnly.Rds
+++ b/inst/tinytest/reference_codaOnly.Rds
diff --git a/inst/tinytest/reference_noDIC.Rds b/inst/tinytest/reference_noDIC.Rds
index 085dad8..74ae4fe 100644
--- a/inst/tinytest/reference_noDIC.Rds
+++ b/inst/tinytest/reference_noDIC.Rds
diff --git a/inst/tinytest/reference_parsorder.Rds b/inst/tinytest/reference_parsorder.Rds
index dda02ca..8ebc9c3 100644
--- a/inst/tinytest/reference_parsorder.Rds
+++ b/inst/tinytest/reference_parsorder.Rds
diff --git a/inst/tinytest/reference_parsorder_noDIC.Rds b/inst/tinytest/reference_parsorder_noDIC.Rds
index f6d4d5a..f177143 100644
--- a/inst/tinytest/reference_parsorder_noDIC.Rds
+++ b/inst/tinytest/reference_parsorder_noDIC.Rds
diff --git a/inst/tinytest/test_jagsbasic.R b/inst/tinytest/test_jagsbasic.R
new file mode 100644
index 0000000..e235722
--- /dev/null
+++ b/inst/tinytest/test_jagsbasic.R
@@ -0,0 +1,41 @@
+set.seed(123)
+
+data(longley)
+data <- list(gnp=longley$GNP, employed=longley$Employed, n=nrow(longley))
+
+modfile <- tempfile()
+writeLines("
+model{
+  for (i in 1:n){ 
+    employed[i] ~ dnorm(mu[i], tau)     
+    mu[i] <- alpha + beta*gnp[i] 
+  }
+  alpha ~ dnorm(0, 0.00001)
+  beta ~ dnorm(0, 0.00001)
+  sigma ~ dunif(0,1000)
+  tau <- pow(sigma,-2)
+}", con=modfile)
+
+inits <- function(){  
+  list(alpha=rnorm(1,0,1),beta=rnorm(1,0,1),sigma=runif(1,0,3))  
+}
+params <- c('alpha','beta','sigma', 'mu')     
+
+out <- jags.basic(data = data, inits = inits, parameters.to.save = params,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 100,
+            n.burnin = 50, n.thin = 2, verbose=FALSE)
+
+ref <- readRDS("jagsbasic_reference_fit.Rds")
+
+expect_identical(out, ref)
+
+# Saved model and reordered parameter names
+params <- c('beta', 'alpha', 'sigma', 'mu')     
+out <- jags.basic(data = data, inits = inits, parameters.to.save = params,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 100,
+            n.burnin = 50, n.thin = 2, verbose=FALSE, save.model=TRUE)
+ref <- readRDS("jagsbasic_ref_saved.Rds")
+
+expect_identical(names(out), names(ref))
+out$model <- ref$model
+expect_identical(out, ref)
diff --git a/inst/tinytest/test_process_output.R b/inst/tinytest/test_process_output.R
new file mode 100644
index 0000000..fae78d5
--- /dev/null
+++ b/inst/tinytest/test_process_output.R
@@ -0,0 +1,242 @@
+process_output <- jagsUI:::process_output
+calc_stats <- jagsUI:::calc_stats
+fill_array <- jagsUI:::fill_array
+get_posterior_array <- jagsUI:::get_posterior_array
+sims_list <- jagsUI:::sims_list
+get_stat_array <- jagsUI:::get_stat_array
+all_stat_arrays <- jagsUI:::all_stat_arrays
+overlap_0 <- jagsUI:::overlap_0
+calc_f <- jagsUI:::calc_f
+calc_Rhat <- jagsUI:::calc_Rhat
+calc_neff <- jagsUI:::calc_neff
+calc_param_stats <- jagsUI:::calc_param_stats
+calc_DIC <- jagsUI:::calc_DIC
+which_params <- jagsUI:::which_params
+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)
+expect_inherits(out,'list')
+expect_equal(length(out),15)
+expect_equal(names(out),c("sims.list", "mean", "sd", "q2.5", "q25", "q50", 
+                            "q75", "q97.5","overlap0", "f","Rhat","n.eff",
+                            "pD", "DIC", "summary"))
+expect_true(all(sapply(out$overlap, is.logical)))
+expect_inherits(out$sims.list,'list')
+expect_equal(out$sims.list, sims_list(samples))
+expect_equal(length(out$sims.list),length(param_names(samples,simplify=T)))
+cs <- calc_stats(samples)
+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_true(all(is.na(result)))
+
+#DIC=FALSE
+out2 <- process_output(samples[,-coda::nvar(samples)], 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)
+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))))
+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)
+old_po <- readRDS("old_process_output.Rds")  
+expect_identical(new_po$summary, old_all$summary)
+new_po$summary <- NULL
+expect_identical(new_po, old_po)
+
+# test that fill_array works properly------------------------------------------
+dat <- 1:10
+indices <- matrix(c(rep(1:5,2),rep(1:2,each=5)),ncol=2)
+out_mat <- matrix(dat,ncol=2)
+expect_equal(fill_array(dat,indices),out_mat)
+dat <- 1:20
+indices <- matrix(c(rep(1:5,4),
+                      rep(rep(1:2,each=5),2),
+                      rep(1:2,each=10)),ncol=3)
+out_arr <- array(dat,dim=c(5,2,2))
+expect_equal(fill_array(dat,indices),out_arr)
+
+
+# test that get_posterior_array output structure is correct--------------------
+samples <- readRDS('coda_samples.Rds')
+n_chains <- length(samples)
+n_samples <- nrow(samples[[1]])
+out <- get_posterior_array('alpha',samples)
+expect_equal(class(out),'numeric')
+expect_equal(length(out),n_samples*n_chains)
+out <- get_posterior_array('mu',samples)
+expect_true(inherits(out,'matrix'))
+expect_equal(dim(out),c(n_samples*n_chains,16))
+out <- get_posterior_array('kappa',samples)
+expect_equal(class(out),'array')
+expect_equal(dim(out),c(n_samples*n_chains,2,2,2))
+#Check error handling
+expect_message(result <- get_posterior_array('fake', samples))
+expect_true(is.na(result))
+
+
+# test that sims_list generates correct sims.list------------------------------
+samples <- readRDS('coda_samples.Rds')
+out <- sims_list(samples)
+expect_equal(class(out),'list')
+expect_equal(names(out),param_names(samples,T))
+expect_equal(sapply(out, function(x) class(x)[1]),
+             c(alpha = "numeric", beta = "numeric", sigma = "numeric",
+               mu = "matrix", kappa = "array", deviance = "numeric"))
+expect_equal(sapply(out,dim), list(alpha = NULL, beta = NULL, sigma = NULL, 
+                                  mu = c(90L, 16L), 
+                                  kappa = c(90L, 2L, 2L, 2L), 
+                                  deviance = NULL))
+
+
+# test that get_stat_array generates correct array-----------------------------
+samples <- readRDS('coda_samples.Rds')
+sum <- calc_stats(samples)
+mean_alpha <- get_stat_array('alpha','mean',sum)
+expect_equal(mean_alpha, 51.90933, tol=1e-4)
+f_mu <- get_stat_array('mu','f',sum)
+expect_equal(f_mu, structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+                                        1, 1, 1, 1, 1), .Dim = 16L))
+sd_kappa <- get_stat_array('kappa', 'sd', sum)
+expect_equal(sd_kappa, structure(c(3.24041434459656, 2.99177153512039, 
+                                   3.2912326732425, 3.35708890821258, 
+                                   3.01174961734256, 3.34886787628231, 
+                                   2.97520825307743, 3.16364214294695), 
+                                 .Dim = c(2L, 2L, 2L)))
+#Check error handling
+expect_message(result <- get_stat_array('fake','mean',sum))
+expect_true(is.na(result))
+
+#Exclude
+ex <- calc_stats(samples, coda_only="mu")
+mu_rows <- ex[grepl("mu", rownames(ex)),]
+expect_equal(nrow(mu_rows), 16)
+expect_equal(ex[,"mean"], sum[,"mean"])
+expect_true(all(is.na(mu_rows[,"sd"])))
+
+
+# test that all_stat_arrays makes a list of stat arrays for all params---------
+samples <- readRDS('coda_samples.Rds')
+sumstat <- calc_stats(samples, coda_only=NULL)
+sal <- all_stat_arrays(sumstat, coda_only=NULL)
+expect_inherits(sal, 'list')
+expect_equal(names(sal),colnames(sumstat))
+expect_true(all(sapply(sal,class)=='list'))
+expect_equal(sal$mean$alpha, 51.90933, tol=1e-4)
+expect_equal(dim(sal$sd$kappa), c(2,2,2))
+
+#With some parameters excluded
+sum_sub <- calc_stats(samples, coda_only=c('alpha','mu','kappa'))
+sal_sub <- all_stat_arrays(sum_sub, coda_only=c("alpha","mu","kappa"))
+expect_equal(names(sal_sub),colnames(sumstat))
+expect_false(any(is.na(unlist(sal_sub$mean))))
+expect_equal(dim(sal_sub$mean$kappa), c(2,2,2))
+expect_true(is.na(sal_sub$sd$kappa))
+
+
+# test that overlap0 calculation is correct------------------------------------
+expect_equal(overlap_0(-2.5, 3.5), 1)
+expect_equal(overlap_0(1.5, 3.5), 0)
+expect_equal(overlap_0(-3,0), 1)
+expect_equal(overlap_0(0, 2.5), 1)
+
+
+# test that calculation of f statistic is correct------------------------------
+set.seed(123)
+test <- c(runif(10,-3,-1),runif(20,1,3))
+test <- matrix(test, nrow=10)
+expect_equal(calc_f(test, mean(test)), 2/3 )
+set.seed(123)
+test <- c(runif(10,-10,-8),runif(20,1,2))
+test <- matrix(test, nrow=10)
+expect_equal(calc_f(test, mean(test)), 1/3 )
+
+test[1,1] <- NA
+expect_equal(round(calc_f(test, mean(test,na.rm=T)),4), 0.3103)
+
+
+# test that all stats for one parameter calculated correctly-------------------
+samples <- readRDS('coda_samples.Rds')
+ps <- calc_param_stats(samples[,'alpha'], FALSE)
+expect_equal(length(ps), 11)
+expect_equivalent(ps["mean"], mean(unlist(samples[,"alpha"])))
+expect_equivalent(ps["sd"], sd(unlist(samples[,"alpha"])))
+expect_equivalent(ps["q2.5"], quantile(unlist(samples[,"alpha"]), 0.025))
+expect_equivalent(ps["Rhat"], coda::gelman.diag(samples[,"alpha"])$psrf[1])
+expect_equivalent(ps, c(51.90933,0.89475,50.31759,51.17975,51.9649,52.46763,
+                   53.6077087,0,1,1.0038308,90), tol=1e-4)
+#Test if inf value is present
+alpha_inf <- samples[,"alpha"]
+alpha_inf[[1]][1] <- Inf
+expect_equivalent(calc_param_stats(alpha_inf, FALSE), rep(NA,11))
+#Test if NA is present
+alpha_na <- samples[,"alpha"]
+alpha_na[[1]][1] <- NA
+expect_equivalent(calc_param_stats(alpha_na, FALSE), 
+             c(51.9180331,0.89598,50.309243,51.198956,51.98700,52.47200,
+               53.60957,0,1,NA,90),tol=1e-4)
+#Test if all NA
+alpha_na[[1]][] <- NA
+alpha_na[[2]][] <- NA
+alpha_na[[3]][] <- NA
+expect_equivalent(calc_param_stats(alpha_na, FALSE), rep(NA, 11))
+#Test if one row
+alpha_one <- samples[,"alpha"]
+alpha_one[[1]] <- coda::mcmc(alpha_one[[1]][1])
+alpha_one[[2]] <- coda::mcmc(alpha_one[[2]][1])
+alpha_one[[3]] <- coda::mcmc(alpha_one[[3]][1])
+expect_equivalent(calc_param_stats(alpha_one, FALSE),
+             c(51.870939,0.8998954,51.15826,51.36934,51.6038732,52.239005,
+               52.8106251, 0, 1, NA, NA), tol=1e-4)
+#Test if error
+alpha_er <- samples[,"alpha"]
+alpha_er[[1]][1] <- 'a'
+expect_warning(out <- calc_param_stats(alpha_er, TRUE))
+expect_true(all(is.na(out)))
+expect_true(length(out) == 11)
+
+
+# test that stats for all parameters are calculated by calc_stats--------------
+samples <- readRDS('coda_samples.Rds')
+st <- calc_stats(samples)
+expect_equal(dim(st), c(length(param_names(samples)), 11))
+expect_equal(rownames(st), param_names(samples))
+expect_equal(colnames(st), c('mean','sd','q2.5','q25','q50','q75','q97.5',
+                             'overlap0','f','Rhat','n.eff'))
+ref_output <- readRDS('calc_stats_output.Rds')
+expect_equal(st, ref_output)
+
+
+# test that calculating stats for a subset of parameters works-----------------
+samples <- readRDS('coda_samples.Rds') 
+ref_output <- readRDS('calc_stats_output.Rds')
+ref_output[c(1,4:19),2:11] <- NA
+st_sub <- calc_stats(samples, coda_only=c('alpha','mu'))
+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)
+dev_ind <- which_params('deviance', param_names(samples))
+no_dev <- samples[,-coda::nvar(samples)]
+expect_true(is.null(calc_DIC(no_dev)))
+samp_na <- samples
+ind <- which_params('deviance',param_names(samples))
+samp_na[[1]][1,ind] <- NA 
+expect_true(is.null(calc_DIC(samp_na)))
+samp_inf <- samples
+samp_inf[[1]][1,ind] <- Inf
+expect_true(all(is.na(calc_DIC(samp_na))))
+samp_inf[[1]][1,ind] <- -Inf
+expect_true(is.null(calc_DIC(samp_na)))
diff --git a/man/jags_View.Rd b/man/jags_View.Rd
index e69e5a8..ed8221a 100644
--- a/man/jags_View.Rd
+++ b/man/jags_View.Rd
@@ -18,8 +18,6 @@
 
 \description{
 Show an R object in a separate, spreadsheet-style window via a call to \code{\link[utils]{View}}.
-
-This function replaces the \code{View} method for class \code{jagsUI}, which caused problems with RStudio.
 }
 
 \author{