Refactor process_output and related functions

And add associated tests.

48 files changed, 1211 insertions, 572 deletions

diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index 5c67090..f7a31f8 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -1,14 +1,14 @@
-# For help debugging build failures open an issue on the RStudio community with the 'github-actions' tag.
-# https://community.rstudio.com/new-topic?category=Package%20development&tags=github-actions
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+#
+# NOTE: This workflow is overkill for most R packages and
+# check-standard.yaml is likely a better choice.
+# usethis::use_github_action("check-standard") will install it.
 on:
   push:
-    branches:
-      - main
-      - master
+    branches: [main, master]
   pull_request:
-    branches:
-      - main
-      - master
+    branches: [main, master]
 
 name: R-CMD-check
 
@@ -27,62 +27,25 @@ jobs:
           - {os: ubuntu-20.04, r: 'oldrel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
 
     env:
-      R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
-      RSPM: ${{ matrix.config.rspm }}
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
 
-      - uses: r-lib/actions/setup-pandoc@v1
-
-      - name: Query dependencies
-        run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
-        shell: Rscript {0}
-
-      - name: Cache R packages
-        if: runner.os != 'Windows'
-        uses: actions/cache@v2
+      - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
-
-      - name: Install JAGS
-        run: |
-          sudo apt-get update -y
-          sudo apt-get install -y jags
-
-      - name: Install system dependencies
-        if: runner.os == 'Linux'
-        run: |
-          while read -r cmd
-          do
-            eval sudo $cmd
-          done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))')
-
-      - name: Install dependencies
-        run: |
-          remotes::install_deps(dependencies = TRUE)
-          remotes::install_cran("rcmdcheck")
-        shell: Rscript {0}
-
-      - name: Check
-        env:
-          _R_CHECK_CRAN_INCOMING_REMOTE_: false
-        run: rcmdcheck::rcmdcheck(args = c("--no-manual", "--ignore-vignettes", "--no-build-vignettes", "--as-cran"), build_args = c("--no-manual", "--no-build-vignettes"), error_on = "warning", check_dir = "check")
-        shell: Rscript {0}
+          extra-packages: any::rcmdcheck
+          needs: check
 
-      - name: Upload check results
-        if: failure()
-        uses: actions/upload-artifact@main
+      - uses: r-lib/actions/check-r-package@v2
         with:
-          name: ${{ runner.os }}-r${{ matrix.config.r }}-results
-          path: check
+          upload-snapshots: true
diff --git a/DESCRIPTION b/DESCRIPTION
index 9e3a7b4..8fe4d68 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: jagsUI
-Version: 1.5.2.9002
-Date: 2021-09-27
+Version: 1.5.3.9000
+Date: 2023-12-03
 Title: A Wrapper Around 'rjags' to Streamline 'JAGS' Analyses
 Authors@R: c(
     person("Ken", "Kellner", email="contact@kenkellner.com", role=c("cre","aut")),
@@ -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/autojags.R b/R/autojags.R
index d880dc9..15b0e9d 100644
--- a/R/autojags.R
+++ b/R/autojags.R
@@ -128,13 +128,12 @@ autojags <- function(data,inits=NULL,parameters.to.save,model.file,n.chains,n.ad
   date <- start.time
   
   #Reorganize JAGS output to match input parameter order
-  samples <- order.params(samples,parameters.to.save,DIC,verbose=verbose)
+  samples <- order_samples(samples, parameters.to.save)
   
   #Convert rjags output to jagsUI form 
-  output <- process.output(samples,DIC=DIC,codaOnly,verbose=verbose) 
+  output <- process_output(samples, coda_only = codaOnly, DIC, quiet = !verbose)
   if(is.null(output)){
     output <- list()
-    samples <- order.params(samples,parameters.to.save,DIC,verbose=verbose)
     output$samples <- samples
     output$model <- mod
     output$n.cores <- n.cores
@@ -143,10 +142,6 @@ autojags <- function(data,inits=NULL,parameters.to.save,model.file,n.chains,n.ad
   }
 
   #Add additional information to output list
-  
-  #Summary
-  output$summary <- summary.matrix(output,samples,n.chains,codaOnly)
-  
   output$samples <- samples
   output$modfile <- model.file
   #If user wants to save input data/inits
diff --git a/R/densityplot.R b/R/densityplot.R
index 0c6f987..b350f0a 100644
--- a/R/densityplot.R
+++ b/R/densityplot.R
@@ -23,7 +23,7 @@ densityplot <- function(x, parameters=NULL, layout=NULL, ask=NULL){
 param_density <- function(x, parameter, m_labels=FALSE){
   
   #Get samples
-  vals <- mcmc_to_mat(x$samples, parameter)
+  vals <- mcmc_to_mat(x$samples[,parameter])
   
   # Get bandwidth, one value for all chains
   bw <- mean(apply(vals, 2, stats::bw.nrd0))
diff --git a/R/getdim.R b/R/getdim.R
deleted file mode 100644
index 382b658..0000000
--- a/R/getdim.R
+++ /dev/null
@@ -1,34 +0,0 @@
-#This function gets the dimensions of non-scalar parameters 
-#for which the user has requested posterior distributions.
-
-get.dim <- function(params){
-  
-  #Get all unique parameters (i.e., collapse indexed non-scalars)
-  ps <- unique(sapply(strsplit(params, "\\["), "[", 1)) 
-  #Slice indexes from non-scalar parameter entries
-  test <- sapply(strsplit(params, "\\["), "[", 1)
-  
-  #Calculate dimension for each parameter i
-  dim <- lapply(ps, function(i){
-    
-    #Extract indices from each element j of parameter i
-    w <- params[test==i]
-    getinds <- lapply(w,FUN=function(j){
-      
-      w2 <- strsplit(j,'\\[')[[1]][2]
-      w3 <- strsplit(w2,"\\]")[[1]] 
-      w4 <- as.numeric(unlist(strsplit(w3,",")))
-      return(w4)
-      
-    })
-    
-    #Get max value from each dimension of i
-    collapsedinds <- do.call(rbind,getinds)
-    apply(collapsedinds,2,max)  
-        
-  })
-  
-  names(dim) = ps
-  dim
-  
-}
\ No newline at end of file
diff --git a/R/jags.R b/R/jags.R
index 3436e0e..9523575 100644
--- a/R/jags.R
+++ b/R/jags.R
@@ -61,15 +61,12 @@ jagsUI <- jags <- function(data,inits=NULL,parameters.to.save,model.file,n.chain
   if(parallel){mcmc.info$n.cores <- n.cores}
   
   #Reorganize JAGS output to match input parameter order
-  if(dim(samples[[1]])[2]>1){
-    samples <- order.params(samples,parameters.to.save,DIC,verbose=verbose)
-  }
+  samples <- order_samples(samples, parameters.to.save)
   
   #Convert rjags output to jagsUI form 
-  output <- process.output(samples,DIC=DIC,codaOnly,verbose=verbose)
+  output <- process_output(samples, coda_only = codaOnly, DIC, quiet = !verbose)
   if(is.null(output)){
     output <- list()
-    samples <- order.params(samples,parameters.to.save,DIC,verbose=verbose)
     output$samples <- samples
     output$model <- m
     output$n.cores <- n.cores
@@ -78,10 +75,6 @@ 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$samples <- samples
   output$modfile <- model.file
   #If user wants to save input data/inits
diff --git a/R/jagsbasic.R b/R/jagsbasic.R
index a086019..17d465a 100644
--- a/R/jagsbasic.R
+++ b/R/jagsbasic.R
@@ -54,7 +54,7 @@ jags.basic <- function(data,inits=NULL,parameters.to.save,model.file,n.chains,n.
   
   if(save.model){
   output <- list()
-  samples <- order.params(samples,parameters.to.save,DIC,verbose=verbose)
+  samples <- order_samples(samples, parameters.to.save)
   output$samples <- samples
   output$model <- m
   output$n.cores <- n.cores
diff --git a/R/mcmc_tools.R b/R/mcmc_tools.R
index 205631d..612fda5 100644
--- a/R/mcmc_tools.R
+++ b/R/mcmc_tools.R
@@ -1,40 +1,22 @@
-#Functions for manipulating and extracting info from mcmc.list-class objects
-#from package rjags/coda
+#------------------------------------------------------------------------------
+#Get names of parameters from an mcmc.list
+#If simplify=T, also drop brackets/indices
+param_names <- function(mcmc_list, simplify=FALSE){
+  out <- coda::varnames(mcmc_list)
+  if(simplify) out <- strip_params(out, unique=TRUE)
+  out
+}
 
-# This is a subset of the functions in mcmc_tools in devel version 1.5.1.9024
 
-###------------------------------------------------------------------------------
-#Remove brackets and indices from parameter names in mcmc.list
+#------------------------------------------------------------------------------
 strip_params <- function(params_raw, unique=FALSE){
-  params_strip <- sapply(strsplit(params_raw,'[', fixed=T),'[',1)
+  params_strip <- sapply(strsplit(params_raw,'[', fixed=TRUE),'[',1)
   if(unique) return( unique(params_strip) )
   params_strip
 }
-#------------------------------------------------------------------------------
 
-###------------------------------------------------------------------------------
-#Identify which columns in mcmc.list object correspond to a given
-#parameter name (useful for non-scalar parameters)
-which_params <- function(param, params_raw){
-  params_strip <- strip_params(params_raw)
-  if( ! param %in% params_strip ){
-    return(NULL)
-  }
-  which(params_strip == param)
-}
-#------------------------------------------------------------------------------
 
-###------------------------------------------------------------------------------
-#Get names of parameters from an mcmc.list
-#If simplify=T, also drop brackets/indices
-param_names <- function(mcmc_list, simplify=FALSE){
-  raw <- colnames(mcmc_list[[1]])
-  if(!simplify) return(raw)
-  strip_params(raw, unique=T)
-}
 #------------------------------------------------------------------------------
-
-###------------------------------------------------------------------------------
 #Match parameter name to scalar or array versions of parameter name
 match_params <- function(params, params_raw){
   unlist(lapply(params, function(x){
@@ -43,28 +25,58 @@ match_params <- function(params, params_raw){
     params_raw[which_params(x, params_raw)]
     }))
 }
+
+
 #------------------------------------------------------------------------------
+#Reorder output samples from coda to match input parameter order
+order_samples <- function(samples, params){
+  tryCatch({
+    matched <- match_params(params, param_names(samples))
+    if("deviance" %in% param_names(samples) & ! "deviance" %in% matched){
+      matched <- c(matched, "deviance")
+    }
+    samples[,matched,drop=FALSE]
+  }, error = function(e){ 
+    message(paste0("Caught error re-ordering samples:\n",e,"\n"))
+    samples
+  })
+}
 
-###------------------------------------------------------------------------------
-#Subset cols of mcmc.list (simple version of [.mcmc.list method)
-select_cols <- function(mcmc_list, col_inds){
-  out <- lapply(1:length(mcmc_list), FUN=function(x){
-            mcmc_element <- mcmc_list[[x]][,col_inds,drop=FALSE]
-            attr(mcmc_element,'mcpar') <- attr(mcmc_list[[x]], 'mcpar')
-            class(mcmc_element) <- 'mcmc'
-            mcmc_element
-            })
-  class(out) <- 'mcmc.list'
-  out
+
+#------------------------------------------------------------------------------
+#Identify which columns in mcmc.list object correspond to a given
+#parameter name (useful for non-scalar parameters)
+which_params <- function(param, params_raw){
+  params_strip <- strip_params(params_raw)
+  if( ! param %in% params_strip ){
+    return(NULL)
+  } 
+  which(params_strip == param)
 }
+
+
 #------------------------------------------------------------------------------
+mcmc_to_mat <- function(mcmc_list){
+  stopifnot(coda::nvar(mcmc_list) == 1)
+  matrix(unlist(mcmc_list), 
+         nrow=coda::niter(mcmc_list), ncol=coda::nchain(mcmc_list))
+}
 
-###------------------------------------------------------------------------------
-#Convert one parameter in mcmc.list to matrix, n_iter * n_chains
-mcmc_to_mat <- function(samples, param){
-  psamples <- select_cols(samples, param)
-  n_chain <- length(samples)
-  n_iter <- nrow(samples[[1]])
-  matrix(unlist(psamples), nrow=n_iter, ncol=n_chain)
+#------------------------------------------------------------------------------
+#Extract index values inside brackets from a non-scalar parameter
+#param is the "base" name of the parameter and params_raw is a vector of 
+#strings that contain brackets
+get_inds <- function(param, params_raw){
+  inds_raw <- sub(paste(param,'[',sep=''),'', params_raw,fixed=T)
+  inds_raw <- sub(']','', inds_raw, fixed=T)
+  inds_raw <- strsplit(inds_raw,',',fixed=T)
+  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/orderparams.R b/R/orderparams.R
deleted file mode 100644
index f50c5c6..0000000
--- a/R/orderparams.R
+++ /dev/null
@@ -1,19 +0,0 @@
-
-order.params <- function(samples,parameters.to.save,DIC,verbose=TRUE){
-  
-  params <- colnames(samples[[1]])
-  params <- params[order(match(sapply(strsplit(params, "\\["), "[", 1),
-                               sapply(strsplit(parameters.to.save, "\\["), "[", 1)))]
-  
-  if(DIC&&('deviance'%in%params)){
-    params <- c(params[params!='deviance'],'deviance')
-  } else if (DIC&&!('deviance'%in%params)){
-    if(verbose){warning('JAGS did not monitor deviance.')}
-    DIC <- FALSE
-  } 
-  
-  samples <- samples[,params, drop=FALSE]
-  
-  return(samples)
-  
-}
diff --git a/R/populate.R b/R/populate.R
deleted file mode 100644
index 15e1a26..0000000
--- a/R/populate.R
+++ /dev/null
@@ -1,36 +0,0 @@
-
-populate <- function(input,dim,simslist=FALSE,samples=NULL){
-  
-  if(!simslist){
-    
-    charinds <- sub(".*\\[(.*)\\].*", "\\1", names(input), perl=TRUE) 
-  
-    fill <- array(NA,dim=dim)
-    
-    for (i in 1:length(input)){
-      
-      ind <- lapply(strsplit(charinds[i], ','), as.integer)[[1]]
-      fill[matrix(ind,1)] <- input[i]
-    
-    }
-  } else {
-    
-    charinds <- sub(".*\\[(.*)\\].*", "\\1", colnames(input), perl=TRUE) 
-    
-    fill <- array(NA,dim=c(samples,dim))
-    
-    for (i in 1:length(charinds)){
-      
-      #ind <- lapply(strsplit(charinds[i], ','), as.integer)[[1]]
-      
-      eval(parse(text=paste('fill[','1:',samples,',',charinds[i],']','<- input[,i]',sep="")))
-      
-    }
-  }
-
-  return(fill)
-  
-}
-
-
-
diff --git a/R/ppcheck.R b/R/ppcheck.R
index 8229b7b..aae4408 100644
--- a/R/ppcheck.R
+++ b/R/ppcheck.R
@@ -11,8 +11,8 @@ pp.check <- function(x, observed, simulated,
     stop("Simulated parameter not found in output")
   }
 
-  obs <- c(mcmc_to_mat(x$samples, observed))
-  sim <- c(mcmc_to_mat(x$samples, simulated))
+  obs <- c(mcmc_to_mat(x$samples[,observed]))
+  sim <- c(mcmc_to_mat(x$samples[,simulated]))
 
   bpval <- mean(sim > obs)
   plotrange <- range(obs, sim)
diff --git a/R/print.R b/R/print.R
index 1d3030d..bb593ec 100644
--- a/R/print.R
+++ b/R/print.R
@@ -1,76 +1,83 @@
 
 print.jagsUI <- function(x,digits=3,...){
-  
+
+  mc <- x$mcmc.info
+
+  # Header
   #bugs.format=TRUE prints a nearly exact replica of WinBUGS-style output
-  
-  #Header
-  if(!x$bugs.format){
-    cat('JAGS output for model \'',x$modfile,'\', generated by jagsUI.','\n',sep="")
-    cat('Estimates based on',x$mcmc.info$n.chains,'chains of',x$mcmc.info$n.iter,'iterations,\n') 
-    if(all(x$mcmc.info$sufficient.adapt)){cat('adaptation =',mean(x$mcmc.info$n.adapt),'iterations (sufficient),\n')
-      } else{cat('adaptation =',mean(x$mcmc.info$n.adapt),'iterations (possibly insufficient),\n')}
-    cat('burn-in = ',x$mcmc.info$n.burnin,' iterations and thin rate = ',x$mcmc.info$n.thin,',','\n',sep="")
-    cat('yielding',x$mcmc.info$n.samples,'total samples from the joint posterior.','\n')
-    if(!x$parallel){cat('MCMC ran for ',x$mcmc.info$elapsed.mins,' minutes at time ',paste(x$run.date),'.\n','\n',sep="")
-    } else{cat('MCMC ran in parallel for ',x$mcmc.info$elapsed.mins,' minutes at time ',paste(x$run.date),'.\n','\n',sep="")}
-  } else{
+  if(x$bugs.format){
     cat('Inference for Bugs model at \'',x$modfile,'\', fit using JAGS,','\n',sep="")
-    cat(x$mcmc.info$n.chains,'chains, each with',x$mcmc.info$n.iter,'iterations (first ',x$mcmc.info$n.burnin,'discarded), n.thin =',x$mcmc.info$n.thin) 
-    cat('\nn.sims =',x$mcmc.info$n.samples,'iterations saved','\n')    
-  }
-  
-  #Organize columns
-  if(!x$bugs.format){
-    if(x$mcmc.info$n.chains!=1){y = x$summary[,c(1,2,3,5,7,10,11,8,9)]
-    } else {y = x$summary[,c(1,2,3,5,7,8,9)]}
-    z <-  as.data.frame(round(as.matrix(y),digits))
-    if(is.vector(y)){
-      z <- as.data.frame(t(z))
-      row.names(z) <- rownames(x$summary)
-    }
-    z[,6] <- z[,6]==1
+    cat(mc$n.chains,'chains, each with',mc$n.iter,'iterations (first ',mc$n.burnin,'discarded), n.thin =',mc$n.thin) 
+    cat('\nn.sims =',mc$n.samples,'iterations saved','\n')  
   } else {
-    if(x$mcmc.info$n.chains!=1){y = x$summary[,c(1:9)]
-    } else {y = x$summary[,c(1:7)]}
-    z <- as.data.frame(round(as.matrix(y),digits))
-    if(is.vector(y)){
-      z <- as.data.frame(t(z))
-      row.names(z) <- rownames(x$summary)
+    cat('JAGS output for model \'',x$modfile,'\', generated by jagsUI.','\n',sep="")
+    cat('Estimates based on',mc$n.chains,'chains of',mc$n.iter,'iterations,\n') 
+    if(all(mc$sufficient.adapt)){
+      cat('adaptation =',mean(mc$n.adapt),'iterations (sufficient),\n')
+    } else{
+      cat('adaptation =',mean(mc$n.adapt),'iterations (possibly insufficient),\n')
+    }
+    cat('burn-in = ',mc$n.burnin,' iterations and thin rate = ',mc$n.thin,',','\n',sep="")
+    cat('yielding',mc$n.samples,'total samples from the joint posterior.','\n')
+    if(!x$parallel){
+      cat('MCMC ran for ',mc$elapsed.mins,' minutes at time ',paste(x$run.date),'.\n','\n',sep="")
+    } else{
+      cat('MCMC ran in parallel for ',mc$elapsed.mins,' minutes at time ',paste(x$run.date),'.\n','\n',sep="")
     }
   }
-  
-  #print the output
-  print(z)
-  
+
+  # Table
+  col_idx <- 1:ncol(x$summary)
+  if(x$bugs.format){
+    col_idx <- 1:9
+    if(mc$n.chains == 1) col_idx <- 1:7
+  } else {
+    col_idx <- c(1,2,3,5,7,10,11,8,9)
+    if(mc$n.chains == 1) col_idx <- c(1,2,3,5,7,10,11)
+  }
+  tab <- x$summary[,col_idx,drop=FALSE]
+  tab <- as.data.frame(round(tab, digits))
+  if(!x$bugs.format) tab[,"overlap0"] <- tab[,"overlap0"] == 1
+  print(tab)
+
   #Print Rhat/n.eff information if necessary
-  if(x$mcmc.info$n.chains>1){
-    if(!x$bugs.format){
-      if(max(unlist(x$Rhat),na.rm=TRUE)>1.1){cat('\n**WARNING** Rhat values indicate convergence failure.','\n')
-      }else{cat('\nSuccessful convergence based on Rhat values (all < 1.1).','\n')}
-      cat('Rhat is the potential scale reduction factor (at convergence, Rhat=1).','\n')
-      cat('For each parameter, n.eff is a crude measure of effective sample size.','\n')
-    } else {
+  if(mc$n.chains>1){
+    if(x$bugs.format){
       cat('\nFor each parameter, n.eff is a crude measure of effective sample size,','\n')
       cat('and Rhat is the potential scale reduction factor (at convergence, Rhat=1).','\n')
+    } else {
+      rhats <- unlist(tab[,"Rhat"])
+      if(any(is.na(rhats))){
+        cat('\n**WARNING** Some Rhat values could not be calculated.')
+      }
+      if(all(is.na(rhats))){
+        cat('\n')
+      } else if (max(rhats, na.rm=TRUE) > 1.1){
+        cat('\n**WARNING** Rhat values indicate convergence failure.','\n')
+      } else {
+        cat('\nSuccessful convergence based on Rhat values (all < 1.1).','\n')
+      }
+      cat('Rhat is the potential scale reduction factor (at convergence, Rhat=1).','\n')
+      cat('For each parameter, n.eff is a crude measure of effective sample size.','\n')
     }
   }
-  
+
   #Print overlap0/f statistic info
   if(!x$bugs.format){
-  cat('\noverlap0 checks if 0 falls in the parameter\'s 95% credible interval.\n')
-  cat('f is the proportion of the posterior with the same sign as the mean;\n')
-  cat('i.e., our confidence that the parameter is positive or negative.\n')
+    cat('\noverlap0 checks if 0 falls in the parameter\'s 95% credible interval.\n')
+    cat('f is the proportion of the posterior with the same sign as the mean;\n')
+    cat('i.e., our confidence that the parameter is positive or negative.\n')
   }
   
   #Print DIC info
   if(x$calc.DIC & !is.null(x$pD)){
-    if(!x$bugs.format){
-      cat('\nDIC info: (pD = var(deviance)/2)','\npD =',round(x$pD,1),'and DIC =',round(x$DIC,digits),'\n')
-      cat('DIC is an estimate of expected predictive error (lower is better).\n')
-    } else {
+    if(x$bugs.format){
       cat('\nDIC info (using the rule, pD = var(deviance)/2)','\npD =',round(x$pD,1),'and DIC =',round(x$DIC,digits),'\n')
       cat('DIC is an estimate of expected predictive error (lower deviance is better).\n')
-    }
+    } else {
+      cat('\nDIC info: (pD = var(deviance)/2)','\npD =',round(x$pD,1),'and DIC =',round(x$DIC,digits),'\n')
+      cat('DIC is an estimate of expected predictive error (lower is better).\n')
+    }  
   }
     
     
diff --git a/R/process_output.R b/R/process_output.R
new file mode 100644
index 0000000..a17ff48
--- /dev/null
+++ b/R/process_output.R
@@ -0,0 +1,258 @@
+#------------------------------------------------------------------------------
+#Process output master function
+#To generate backwards-compatible jagsUI output
+process_output <- function(mcmc_list, coda_only=NULL, DIC, quiet=FALSE){
+  if(!quiet){cat('Calculating statistics.......','\n')} 
+
+  tryCatch({
+    if(DIC == -999) stop("Throwing error for testing purposes", call.=FALSE)
+    # 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, DIC)
+    
+    # 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"))
+    NULL
+  })
+}
+
+
+#------------------------------------------------------------------------------
+#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"), drop=FALSE]
+  # 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))
+  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
+  out
+}
+
+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, 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))
+  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/R/processoutput.R b/R/processoutput.R
deleted file mode 100644
index 1aeaca7..0000000
--- a/R/processoutput.R
+++ /dev/null
@@ -1,181 +0,0 @@
-
-process.output <- function(x,DIC,params.omit,verbose=TRUE) {
-
-out <- tryCatch({
-
-if(verbose){cat('Calculating statistics.......','\n')}
-
-# Get parameter names
-params <- colnames(x[[1]])
-
-#Get number of chains
-m <- length(x)
-
-#Collapse mcmc.lists into matrix
-mat = do.call(rbind,x)
-
-#Get # of iterations / chain
-n <- dim(mat)[1] / m
-
-#Get parameter dimensions
-dim <- get.dim(params)
-
-#Create new parameter name vectors to handle non-scalar params
-expand <- sapply(strsplit(params, "\\["), "[", 1)
-params.simple <- unique(sapply(strsplit(params, "\\["), "[", 1))
-
-#Functions for statistics
-qs <- function(x,y){as.numeric(quantile(x,y, na.rm=TRUE))}
-#Overlap 0 function
-ov <- function(x){findInterval(0,sort(c(qs(x,0.025),qs(x,0.975))))==1}
-#f function (proportion of posterior with same sign as mean)
-gf <- function(x){if(mean(x, na.rm=TRUE)>=0){mean(x>=0, na.rm=TRUE)}else{mean(x<0, na.rm=TRUE)}}
-#n.eff function
-calcneff <- function(x,n,m){
-  xp <- matrix(x,nrow=n,ncol=m)
-  xdot <- apply(xp,2,mean, na.rm=TRUE)
-  s2 <- apply(xp,2,var, na.rm=TRUE)
-  W <- mean(s2)
-
-  #Non-degenerate case
-  if (is.na(W)){
-	  n.eff <- NA
-  } else if ((W > 1.e-8) && (m > 1)) {
-    B <- n*var(xdot)
-    sig2hat <- ((n-1)*W + B)/n
-    n.eff <- round(m*n*min(sig2hat/B,1),0)
-    #Degenerate case
-  } else {
-    n.eff <- 1
-  }
-  n.eff
-}
-
-#Gelman diag function
-gd <- function(i,hold){
-  r <- try(gelman.diag(hold[,i], autoburnin=FALSE)$psrf[1], silent=TRUE)
-  if(inherits(r, "try-error") || !is.finite(r)) {
-    r <- NA
-  }
-  return(r)
-}
-
-#Make blank lists
-sims.list <- means <- rhat <- n.eff <- se <- as.list(rep(NA,length(params.simple)))
-q2.5 <- q25 <- q50 <- q75 <- q97.5 <- overlap0 <- f <- as.list(rep(NA,length(params.simple)))
-names(sims.list) <- names(means) <- names(rhat) <- names(n.eff) <- params.simple
-names(se) <- names(q2.5) <- names(q25) <- names(q50) <- names(q75) <- names(q97.5) <- params.simple
-names(overlap0) <- names(f) <- params.simple
-
-#This function modifies objects in global environment (output is discarded)
-#Calculates statistics for each parameter
-calc.stats <- function(i){
-
-  #If parameter is not a scalar (e.g. vector/array)
-  if(!is.na(dim[i][1])){
-
-    #Get all samples
-    sims.list[[i]] <<- mat[,expand==i,drop=FALSE]
-
-	#if every iteration is NA, don't do anything else
-	if(all(is.na(sims.list[[i]]))){return(NA)}
-
-    #If more than 1 chain, calculate rhat
-    #Done separately for each element of non-scalar parameter to avoid errors
-    if(m > 1 && (!i%in%params.omit)){
-      hold <- x[,expand==i,drop=FALSE]
-      nelements <- sum(expand==i)
-      rhat.vals <- sapply(1:nelements,gd,hold=hold)
-      names(rhat.vals) <- colnames(hold[[1]])
-      rhat[[i]] <<- populate(rhat.vals,dim[[i]])
-    } else if (m == 1){
-      hold <- x[,expand==i]
-      rhat[[i]] <<- array(NA,dim=dim[[i]])
-    }
-
-    #Calculate other statistics
-    ld <- length(dim(sims.list[[i]]))
-    means[[i]] <<- populate(colMeans(sims.list[[i]], na.rm=TRUE),dim[[i]])
-    if(!i%in%params.omit){
-    se[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),sd),dim=dim[[i]])
-    q2.5[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),qs,0.025),dim=dim[[i]])
-    q25[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),qs,0.25),dim=dim[[i]])
-    q50[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),qs,0.5),dim=dim[[i]])
-    q75[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),qs,0.75),dim=dim[[i]])
-    q97.5[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),qs,0.975),dim=dim[[i]])
-    overlap0[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),ov),dim=dim[[i]])
-    f[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),gf),dim=dim[[i]])
-    n.eff[[i]] <<- populate(apply(sims.list[[i]],c(2:ld),calcneff,n,m),dim=dim[[i]])
-    }
-
-    sims.list[[i]] <<- populate(sims.list[[i]],dim=dim[[i]],simslist=T,samples=dim(mat)[1])
-
-  #If parameter is a scalar
-  } else {
-
-    if(m > 1 && (!i%in%params.omit)){rhat[[i]] <<- gelman.diag(x[,i],autoburnin=FALSE)$psrf[1]}
-
-    sims.list[[i]] <<- mat[,i]
-
-	if(all(is.na(sims.list[[i]]))){return(NA)}
-
-    means[[i]] <<- mean(sims.list[[i]], na.rm=TRUE)
-    if(!i%in%params.omit){
-    se[[i]] <<- sd(sims.list[[i]], na.rm=TRUE)
-    q2.5[[i]] <<- qs(sims.list[[i]],0.025)
-    q25[[i]] <<- qs(sims.list[[i]],0.25)
-    q50[[i]] <<- qs(sims.list[[i]],0.5)
-    q75[[i]] <<- qs(sims.list[[i]],0.75)
-    q97.5[[i]] <<- qs(sims.list[[i]],0.975)
-    overlap0[[i]] <<- ov(sims.list[[i]])
-    f[[i]] <<- gf(sims.list[[i]])
-    n.eff[[i]] <<- calcneff(sims.list[[i]],n,m)}
-  }
-
-}
-
-#Actually run function(nullout not used for anything)
-nullout <- sapply(params.simple,calc.stats)
-
-#Warn user if at least one Rhat value was NA
-rhat.sub <- unlist(rhat)[!is.na(unlist(means))]
-if(NA%in%rhat.sub&&verbose){
-  options(warn=1)
-  warning('At least one Rhat value could not be calculated.')
-  options(warn=0,error=NULL)
-}
-
-#Do DIC/pD calculations if requested by user
-if(DIC & 'deviance' %in% params){
-  dev <- matrix(data=mat[,'deviance'],ncol=m,nrow=n)
-  pd <- numeric(m)
-  dic <- numeric(m)
-  for (i in 1:m){
-    pd[i] <- var(dev[,i])/2
-    dic[i] <- mean(dev[,i]) + pd[i]
-  }
-  pd <- mean(pd)
-  dic <- mean(dic)
-
-  #Return this list if DIC/pD requested
-  if(verbose){cat('\nDone.','\n')}
-  return(list(sims.list=sims.list,mean=means,sd=se,q2.5=q2.5,q25=q25,q50=q50,q75=q75,q97.5=q97.5,overlap0=overlap0,
-              f=f,Rhat=rhat,n.eff=n.eff,pD=pd,DIC=dic))
-} else {
-  #Otherwise return list without pD/DIC
-  if(verbose){cat('\nDone.','\n')}
-  return(list(sims.list=sims.list,mean=means,sd=se,q2.5=q2.5,q25=q25,q50=q50,q75=q75,q97.5=q97.5,overlap0=overlap0,
-              f=f,Rhat=rhat,n.eff=n.eff))
-}
-
-}, error = function(cond){
-  message('Calculating statistics failed with the following error:')
-  message(cond)
-  message('\nOutput falling back to class jagsUIbasic\n')
-  return(NULL)
-  }
-)
-return(out)
-}
-
-
diff --git a/R/summarymatrix.R b/R/summarymatrix.R
deleted file mode 100644
index 1fc3d06..0000000
--- a/R/summarymatrix.R
+++ /dev/null
@@ -1,55 +0,0 @@
-
-summary.matrix <- function(output,samples,n.chains,codaOnly){
-  
-  hold <- unlist(output$mean[!names(output$mean)%in%codaOnly])
-  toremove <- which(!is.na(hold))
-
-  #Get sorted names
-  sort.names = c()
-  for (i in 1:length(output$mean)){
-    if(length(output$mean[[i]])>1){
-      raw.ind <- which(output$mean[[i]]==output$mean[[i]],arr.ind=T)
-      if(is.matrix(raw.ind)){
-        ind <- apply(raw.ind,1,paste,collapse=',')
-      } else {
-        ind <- raw.ind
-      }
-
-      newnames <- paste(names(output$mean)[i],'[',ind,']',sep='')
-      sort.names <- c(sort.names,newnames)
-
-    } else {
-      sort.names <- c(sort.names,names(output$mean[i]))
-    }
-  }
-
-
-  rnames <- colnames(samples[[1]])
-  sorted.order <- order(match(rnames,sort.names))
-  rnames <- rnames[sorted.order]
-  
-  cleanup <- function(input,codaOnly){
-
-    out.raw <- unlist(input[!names(input)%in%codaOnly])
-    
-    out <- out.raw[toremove]
-    return(out)
-  }
-  
-  y = data.frame(cleanup(output$mean,codaOnly),cleanup(output$sd,codaOnly),
-                 cleanup(output$q2.5,codaOnly),cleanup(output$q25,codaOnly),
-                 cleanup(output$q50,codaOnly),cleanup(output$q75,codaOnly),
-                 cleanup(output$q97.5,codaOnly),
-                 cleanup(output$Rhat,codaOnly),cleanup(output$n.eff,codaOnly),
-                 cleanup(output$overlap0,codaOnly),cleanup(output$f,codaOnly))
- 
-  p <- rnames
-  expand <- sapply(strsplit(p, "\\["), "[", 1)  
-  row.names(y) = p[!expand%in%codaOnly]
-  names(y) = c('mean','sd','2.5%','25%','50%','75%','97.5%','Rhat','n.eff','overlap0','f')
-  if(n.chains==1){
-    y = y[,-c(8,9)]
-  }
-  
-  return(as.matrix(y))
-}
diff --git a/R/traceplot.R b/R/traceplot.R
index 4a7c8bd..755d3aa 100644
--- a/R/traceplot.R
+++ b/R/traceplot.R
@@ -25,7 +25,7 @@ traceplot <- function(x, parameters=NULL, Rhat_min=NULL,
 param_trace <- function(x, parameter, m_labels=FALSE){
 
   #Get samples and Rhat values
-  vals <- mcmc_to_mat(x$samples, parameter)
+  vals <- mcmc_to_mat(x$samples[, parameter])
   Rhat <- sprintf("%.3f",round(x$summary[parameter, 'Rhat'],3))
 
   #Draw plot
diff --git a/R/translateparams.R b/R/translateparams.R
deleted file mode 100644
index e9ad7ce..0000000
--- a/R/translateparams.R
+++ /dev/null
@@ -1,59 +0,0 @@
-
-translate.params <- function(x,params.sub){
-
-
-params = colnames(x$samples[[1]])
-
-params.simple.sub = unique(sapply(strsplit(params.sub, "\\["), "[", 1))
-params.simple <- unique(sapply(strsplit(params, "\\["), "[", 1))
-n = length(params.simple.sub)
-
-if(sum(params.simple.sub%in%params.simple)!=n){stop('One or more specified parameters are not in model output./n')}
-
-
-params.sub.1 <- sapply(strsplit(params.sub, "\\]"), "[", 1)
-params.2 <- sapply(strsplit(params.sub.1, "\\["), "[", 2)
-expand <- sapply(strsplit(params, "\\["), "[", 1)
-
-dim = get.dim(params)
-
-gen.samp.mat <- function(x){
-  out = x
-  for(i in 1:length(x)){
-    if(!is.na(x[[i]][1])){
-      if(length(x[[i]])>1){
-        out[[i]] = array(params[expand==names(x)[i]],dim=x[[i]])
-    }
-    if(length(x[[i]])==1){
-      out[[i]] = params[expand==names(x)[i]]
-    }
-    
-  } else {out[[i]] = NA}
-}
-return(out)
-}
-
-mats = gen.samp.mat(dim)
-
-mats.sub = mats[params.simple.sub]
-
-index=1
-params.new = character()
-for (i in 1:length(params.sub)){
-
-  if(!is.na(mats.sub[i])||!is.na(params.2[i])){
-    if(params.sub[i]==params.simple.sub[i]){
-      st = paste('mats.sub$',params.simple.sub[i],"[]",sep="")
-    } else {
-      st = paste('mats.sub$',params.simple.sub[i],"[",params.2[i],']',sep="")
-    }
-    ind = eval(parse(text=st))
-    params.new[index:(index+length(ind)-1)] = ind
-    index = index+length(ind)
-  } else {
-    params.new[index]=params.sub[i]
-    index=index+1
-  }
-}
-return(params.new)
-}
diff --git a/R/update.R b/R/update.R
index 4576426..f1d2db5 100644
--- a/R/update.R
+++ b/R/update.R
@@ -49,10 +49,10 @@ update.jagsUI <- function(object, parameters.to.save=NULL, n.adapt=NULL, n.iter,
   date <- start.time
   
   #Reorganize JAGS output to match input parameter order
-  samples <- order.params(samples,parameters,DIC,verbose=verbose)
+  samples <- order_samples(samples, parameters)
     
   #Run process output
-  output <- process.output(samples,DIC=DIC,codaOnly,verbose=verbose)
+  output <- process_output(samples, coda_only = codaOnly, DIC, quiet = !verbose)
   if(is.null(output)){
     output <- list()
     output$samples <- samples
@@ -62,9 +62,6 @@ update.jagsUI <- function(object, parameters.to.save=NULL, n.adapt=NULL, n.iter,
     return(output)
   }
  
-  #Summary
-  output$summary <- summary.matrix(output,samples,object$mcmc.info$n.chains,codaOnly)
-  
   #Save other information to output object
   output$samples <- samples
   
diff --git a/R/updatebasic.R b/R/updatebasic.R
index 5b5528d..0a15d17 100644
--- a/R/updatebasic.R
+++ b/R/updatebasic.R
@@ -46,7 +46,7 @@ update.jagsUIbasic <- function(object, parameters.to.save=NULL, n.adapt=NULL, n.
     m <- rjags.output$m    
   }
   
-  samples <- order.params(samples,parameters,DIC)
+  samples <- order_samples(samples, parameters)
   
   end.time <- Sys.time() 
   time <- round(as.numeric(end.time-start.time,units="mins"),digits=3)
diff --git a/R/whiskerplot.R b/R/whiskerplot.R
index 5b2bda9..8b09bf1 100644
--- a/R/whiskerplot.R
+++ b/R/whiskerplot.R
@@ -20,7 +20,7 @@ whiskerplot <- function(x,parameters,quantiles=c(0.025,0.975),
   #Calculate means and CIs
   post_stats <- sapply(parameters, 
                        function(i){
-                          sims <- mcmc_to_mat(x$samples, i)
+                          sims <- mcmc_to_mat(x$samples[, i])
                           c(mean(sims,na.rm=TRUE), 
                             stats::quantile(sims,na.rm=TRUE,quantiles))
                         })
diff --git a/inst/tinytest/autojags_ref.Rds b/inst/tinytest/autojags_ref.Rds
new file mode 100644
index 0000000..3fe83c3
--- /dev/null
+++ b/inst/tinytest/autojags_ref.Rds
diff --git a/inst/tinytest/autojags_ref_codaonly.Rds b/inst/tinytest/autojags_ref_codaonly.Rds
new file mode 100644
index 0000000..0f731e1
--- /dev/null
+++ b/inst/tinytest/autojags_ref_codaonly.Rds
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
new file mode 100644
index 0000000..3d25cc5
--- /dev/null
+++ 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_ref_update.Rds b/inst/tinytest/jagsbasic_ref_update.Rds
new file mode 100644
index 0000000..9472277
--- /dev/null
+++ b/inst/tinytest/jagsbasic_ref_update.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
new file mode 100644
index 0000000..6947a6e
--- /dev/null
+++ 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
new file mode 100644
index 0000000..ff6bc92
--- /dev/null
+++ b/inst/tinytest/one_sample.Rds
diff --git a/inst/tinytest/reference_codaOnly.Rds b/inst/tinytest/reference_codaOnly.Rds
new file mode 100644
index 0000000..4081a59
--- /dev/null
+++ b/inst/tinytest/reference_codaOnly.Rds
diff --git a/inst/tinytest/reference_noDIC.Rds b/inst/tinytest/reference_noDIC.Rds
new file mode 100644
index 0000000..74ae4fe
--- /dev/null
+++ b/inst/tinytest/reference_noDIC.Rds
diff --git a/inst/tinytest/reference_parsorder.Rds b/inst/tinytest/reference_parsorder.Rds
new file mode 100644
index 0000000..8ebc9c3
--- /dev/null
+++ b/inst/tinytest/reference_parsorder.Rds
diff --git a/inst/tinytest/reference_parsorder_noDIC.Rds b/inst/tinytest/reference_parsorder_noDIC.Rds
new file mode 100644
index 0000000..f177143
--- /dev/null
+++ b/inst/tinytest/reference_parsorder_noDIC.Rds
diff --git a/inst/tinytest/test_autojags.R b/inst/tinytest/test_autojags.R
new file mode 100644
index 0000000..d8b6b84
--- /dev/null
+++ b/inst/tinytest/test_autojags.R
@@ -0,0 +1,55 @@
+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')     
+
+nul <- capture.output(
+            out <- autojags(data = data, inits = inits, parameters.to.save = params,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.burnin=50, 
+            iter.increment=10, n.thin = 2, verbose=FALSE))
+ref <- readRDS("autojags_ref.Rds")
+
+# Remove time/date based elements
+out$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out[-c(17,18,21)], ref[-c(17,18,21)])
+
+
+# codaOnly---------------------------------------------------------------------
+nul<- capture.output(
+            out <- autojags(data = data, inits = inits, parameters.to.save = params,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.burnin=50, 
+            iter.increment=10, n.thin = 2, verbose=FALSE, codaOnly=c("mu")))
+ref <- readRDS("autojags_ref_codaonly.Rds")
+
+# Remove time/date based elements
+out$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out[-c(17,18,21)], ref[-c(17,18,21)])
+
+# Check recovery after process_output errors-----------------------------------
+# Setting DIC to -999 forces process_output to error for testing
+expect_message(nul<- capture.output(
+            out <- autojags(data = data, inits = inits, parameters.to.save = params,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.burnin=50, 
+            iter.increment=10, n.thin = 2, verbose=FALSE, codaOnly=c("mu"), DIC=-999)))
+expect_inherits(out, "jagsUIbasic")
+expect_equal(coda::varnames(out$samples), 
+             c("alpha","beta", "sigma", paste0("mu[",1:16,"]"),"deviance"))
+expect_equal(names(out), c("samples", "model"))
diff --git a/inst/tinytest/test_jags.R b/inst/tinytest/test_jags.R
new file mode 100644
index 0000000..9937f1e
--- /dev/null
+++ b/inst/tinytest/test_jags.R
@@ -0,0 +1,216 @@
+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(data = data, inits = inits, parameters.to.save = params,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 1000,
+            n.burnin = 500, n.thin = 2, verbose=FALSE)
+
+# Used below
+mu2_est <- out$mean$mu[2]
+
+ref <- readRDS("longley_reference_fit.Rds")
+
+# Remove time/date based elements
+out$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out[-c(17,18,21)], ref[-c(17,18,21)])
+
+# Plots
+pdf(NULL)
+tp <- traceplot(out, ask=FALSE)
+dev.off()
+expect_true(is.null(tp))
+
+pdf(NULL)
+dp <- densityplot(out, ask=FALSE)
+dev.off()
+expect_true(is.null(dp))
+
+pdf(NULL)
+wp <- whiskerplot(out, "mu")
+dev.off()
+expect_true(is.null(wp))
+
+pdf(NULL)
+pp <- pp.check(out, "alpha", "beta")
+dev.off()
+expect_equal(pp, 0)
+
+# codaOnly---------------------------------------------------------------------
+out <- jags(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 = 1, codaOnly=c("mu", "sigma"), verbose=FALSE)
+ref <- readRDS("reference_codaOnly.Rds")
+
+out$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out[-c(17,18,21)], ref[-c(17,18,21)])
+
+# DIC = FALSE------------------------------------------------------------------
+out <- jags(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 = 1, DIC=FALSE, verbose=FALSE)
+expect_false(out$calc.DIC)
+
+ref <- readRDS("reference_noDIC.Rds")
+out$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out[-c(15,16,19)], ref[-c(15,16,19)])
+
+# Reordered parameter names----------------------------------------------------
+pars_new <- c("mu", "sigma", "alpha", "beta")
+out <- jags(data = data, inits = inits, parameters.to.save = pars_new,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 100,
+            n.burnin = 50, n.thin = 1, verbose=FALSE)
+ref <- readRDS("reference_parsorder.Rds")
+
+out$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out[-c(17,18,21)], ref[-c(17,18,21)])
+
+# Reordered parameter names and no DIC-----------------------------------------
+pars_new <- c("mu", "sigma", "alpha", "beta")
+out <- jags(data = data, inits = inits, parameters.to.save = pars_new,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 100,
+            n.burnin = 50, n.thin = 1, DIC = FALSE, verbose=FALSE)
+ref <- readRDS("reference_parsorder_noDIC.Rds")
+
+out$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out[-c(15,16,19)], ref[-c(15,16,19)])
+
+# Single parameter saved-------------------------------------------------------
+pars_new <- c("alpha")
+out <- jags(data = data, inits = inits, parameters.to.save = pars_new,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 100,
+            n.burnin = 50, n.thin = 1, DIC = FALSE, verbose=FALSE)
+expect_equal(nrow(out$summary), 1)
+expect_equal(ncol(out$samples[[1]]), 1)
+
+# Another example from Github issues
+modfile2 <- tempfile()
+writeLines("
+model{
+  for (i in 1:n) {
+    y[i] ~ dinterval(t[i], c[i, ])
+    t[i] ~ dweib(v, lambda)
+  }
+  v      ~ dunif(0, 10)
+  lambda ~ dunif(0, 10)
+}
+", con=modfile2)
+
+dataList <- list(n = 7L, c = structure(c(5, 7.5, 8, 9, 8, 8, 8, 15, 15, 15,
+15, 10, 10, 10), dim = c(7L, 2L)), y = c(1, 1, 1, 1, 1, 1, 1))
+out <- jags(data = dataList, parameters.to.save = c("v", "lambda"),
+            model.file =modfile2 ,n.chains = 3,n.adapt = 100, n.iter = 100+100,
+            n.burnin = 100, n.thin = 5, verbose=FALSE)
+expect_equal(rownames(out$summary), c("v", "lambda","deviance"))
+
+
+# No non-codaOnly parameters---------------------------------------------------
+
+out <- jags(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 = 1, verbose=FALSE, DIC=FALSE, 
+            codaOnly = params)
+expect_equal(nrow(out$summary), 0)
+
+# Check recovery after process_output errors-----------------------------------
+# Setting DIC to -999 forces process_output to error for testing
+expect_message(out <- jags(data = data, inits = inits, 
+                parameters.to.save = c("alpha","beta"),
+                model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 100,
+                n.burnin = 50, n.thin = 1, verbose=FALSE, DIC=-999))
+expect_inherits(out, "jagsUIbasic")
+expect_equal(coda::varnames(out$samples), c("alpha","beta", "deviance"))
+expect_equal(names(out), c("samples", "model"))
+
+# Single chain and single iteration--------------------------------------------
+out <- jags(data = data, inits = inits, parameters.to.save = params,
+            model.file = modfile, n.chains = 1, n.adapt = 100, n.iter = 100,
+            n.burnin = 99, n.thin = 1, DIC = FALSE, verbose=FALSE)
+expect_true(all(is.na(out$summary[,"sd"])))
+expect_true(all(is.na(out$summary[,"Rhat"])))
+expect_true(all(is.na(out$summary[,"n.eff"])))
+expect_true(all(out$summary["alpha",3:7] == out$summary["alpha",3]))
+
+
+# Single parameter slice-------------------------------------------------------
+set.seed(123)
+pars_new <- c("mu[2]")
+out <- jags(data = data, inits = inits, parameters.to.save = pars_new,
+            model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 1000,
+            n.burnin = 500, n.thin = 2, DIC = FALSE, verbose=FALSE)
+expect_equal(nrow(out$summary), 1)
+expect_equal(ncol(out$samples[[1]]), 1)
+expect_equal(out$mean$mu, mu2_est)
+
+# Ragged arrays----------------------------------------------------------------
+set.seed(123)
+
+# Should trigger creation of a bunch of missing values in mu[,2] in output
+modfile <- tempfile()
+writeLines("
+model{
+  for (i in 1:n){ 
+    employed[i] ~ dnorm(mu[i,1], tau)     
+    mu[i,1] <- alpha + beta*gnp[i] 
+  }
+
+  mu[1,2] <- 1
+
+  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(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 = 1, verbose=FALSE)
+
+expect_equal(dim(out$mean$mu), c(16, 2))
+expect_true(all(is.na(out$mean$mu[2:16,2])))
+expect_equal(dim(out$sims.list$mu), c(150, 16, 2))
+expect_equal(nrow(out$summary), 21)
+expect_equal(rownames(out$summary)[20], "mu[1,2]")
+
+# When no stochastic nodes (and deviance is not calculated)--------------------
+library(jagsUI)
+
+data <- list(a = 1, b = Inf)
+
+modfile <- tempfile()
+writeLines("
+model{
+  x <- a * b
+} ", con=modfile)
+
+expect_warning(out <- jags(data=data, parameters.to.save="x", model.file=modfile,
+                      n.chains=3, n.iter=10, n.burnin=5, n.adapt=10, DIC=TRUE, verbose=FALSE))
+
+expect_equal(coda::varnames(out$samples), c("x"))
+expect_true(is.na(out$Rhat$x))
+expect_true(is.null(out$pD))
+expect_true(is.null(out$DIC))
diff --git a/inst/tinytest/test_jagsbasic.R b/inst/tinytest/test_jagsbasic.R
new file mode 100644
index 0000000..1e03fc4
--- /dev/null
+++ b/inst/tinytest/test_jagsbasic.R
@@ -0,0 +1,50 @@
+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)
+
+# Update
+out2 <- update(out, n.iter=100, n.thin = 2, verbose=FALSE)
+expect_equal(nrow(out2$samples[[1]]), 50)
+
+ref <- readRDS('jagsbasic_ref_update.Rds')
+expect_identical(names(out2), names(ref))
+out2$model <- ref$model
+expect_equal(out2, ref)
diff --git a/inst/tinytest/test_mcmc_tools.R b/inst/tinytest/test_mcmc_tools.R
new file mode 100644
index 0000000..8e5f9cb
--- /dev/null
+++ b/inst/tinytest/test_mcmc_tools.R
@@ -0,0 +1,91 @@
+# test that param_names returns correct names---------------------------------
+param_names <- jagsUI:::param_names
+samples <- readRDS('coda_samples.Rds')
+expect_equal(param_names(samples),
+    c("alpha", "beta", "sigma", "mu[1]", "mu[2]", "mu[3]", "mu[4]",
+    "mu[5]", "mu[6]", "mu[7]", "mu[8]", "mu[9]", "mu[10]", "mu[11]",
+    "mu[12]", "mu[13]", "mu[14]", "mu[15]", "mu[16]", "kappa[1,1,1]",
+    "kappa[2,1,1]", "kappa[1,2,1]", "kappa[2,2,1]", "kappa[1,1,2]",
+    "kappa[2,1,2]", "kappa[1,2,2]", "kappa[2,2,2]", "deviance"))
+expect_equal(param_names(samples,simplify=T),
+             c('alpha','beta','sigma','mu','kappa','deviance'))
+
+
+# test that strip_params removes brackets and indices--------------------------
+strip_params <- jagsUI:::strip_params
+params_raw <- c('alpha','beta[1]','beta[2]','gamma[1,2]','kappa[1,2,3]')
+expect_equal(strip_params(params_raw),
+               c('alpha','beta','beta','gamma','kappa'))
+expect_equal(strip_params(params_raw,unique=T),
+               c('alpha','beta','gamma','kappa'))
+
+
+# test that match_param identifies correct set of params-----------------------
+match_params <- jagsUI:::match_params
+params_raw <- c('alpha','beta[1]','beta[2]','gamma[1,1]','gamma[3,1]')
+expect_equal(match_params('alpha', params_raw),'alpha')
+expect_equal(match_params('beta', params_raw), c('beta[1]','beta[2]'))
+expect_equal(match_params('gamma[1,1]', params_raw), 'gamma[1,1]')
+expect_true(is.null(match_params('fake',params_raw)))
+expect_equal(match_params(c('alpha','beta'),params_raw),
+               c('alpha','beta[1]','beta[2]'))
+expect_equal(match_params(c('alpha','fake','beta'),params_raw),
+               c('alpha','beta[1]','beta[2]'))
+
+# test that order_samples works correctly--------------------------------------
+order_samples <- jagsUI:::order_samples
+samples <- readRDS('coda_samples.Rds')
+new_order <- c('beta','mu','alpha')
+out <- order_samples(samples, new_order)
+expect_equal(class(out), 'mcmc.list')
+expect_equal(length(out),length(samples))
+expect_equal(lapply(out,class),lapply(samples,class))
+expect_equal(param_names(out),c('beta',paste0('mu[',1:16,']'),'alpha', "deviance"))
+expect_equal(dim(out[[1]]), c(30,19))
+expect_equal(as.numeric(out[[1]][1,1:2]), 
+               c(0.03690717, 59.78175), tol=1e-4)
+expect_equal(order_samples(samples, 'beta'),
+             order_samples(samples, c('beta','fake')))
+expect_message(order_samples('fake','beta'))
+expect_message(test <- order_samples('fake','beta'))
+expect_equal(test, 'fake')
+one_param <- samples[, 'alpha',drop=FALSE]
+expect_equal(order_samples(one_param,'alpha'),one_param)
+expect_equal(dim(order_samples(one_param, 'beta')[[1]]),c(30,0))
+new_order <- c('deviance', 'beta','mu','alpha')
+out <- order_samples(samples, new_order)
+expect_equal(param_names(out),c('deviance', 'beta',paste0('mu[',1:16,']'),'alpha'))
+
+# test that which_params gets param col indices--------------------------------
+which_params <- jagsUI:::which_params
+params_raw <- c('alpha','beta[1]','beta[2]','gamma[1,1]','gamma[3,1]')
+expect_equal(which_params('alpha',params_raw),1)
+expect_equal(which_params('beta',params_raw),c(2,3))
+expect_equal(which_params('gamma',params_raw),c(4,5))
+expect_null(which_params('kappa',params_raw))
+
+
+# test that mcmc_to_mat converts properly--------------------------------------
+mcmc_to_mat <- jagsUI:::mcmc_to_mat
+samples <- readRDS('coda_samples.Rds')
+mat <- mcmc_to_mat(samples[, 'alpha'])
+expect_true(inherits(mat, 'matrix'))
+expect_equal(dim(mat),c(nrow(samples[[1]]),length(samples)))
+expect_equal(mat[,1],as.numeric(samples[[1]][,'alpha'])) 
+one_sample <- readRDS('one_sample.Rds')
+mat <- mcmc_to_mat(one_sample[, 'alpha'])
+expect_equal(dim(mat), c(1,3))
+
+# test that get_inds extracts indices-----------------------------------------
+get_inds <- jagsUI:::get_inds
+params_raw <- c('beta[1]','beta[2]')  
+expect_equal(get_inds('beta',params_raw),matrix(c(1,2)))
+params_raw <- c('gamma[1,1]','gamma[2,1]','gamma[1,3]')
+expect_equal(get_inds('gamma',params_raw),
+               matrix(c(1,1,2,1,1,3),ncol=2,byrow=T))
+params_raw <- c('kappa[1,1,1]','kappa[2,1,1]','kappa[1,3,1]')
+expect_equal(get_inds('kappa',params_raw),
+               matrix(c(1,1,1,2,1,1,1,3,1),ncol=3,byrow=T))
+params_raw <- 'alpha'
+expect_warning(test <- get_inds('alpha',params_raw)[1,1])
+expect_true(is.na(test))
diff --git a/inst/tinytest/test_print.R b/inst/tinytest/test_print.R
new file mode 100644
index 0000000..2f765ed
--- /dev/null
+++ b/inst/tinytest/test_print.R
@@ -0,0 +1,48 @@
+out <- readRDS("longley_reference_fit.Rds")
+
+# Check standard output
+pr <- capture.output(print(out))
+expect_true(grepl(" 750 ", pr[5]))
+expect_true(grepl("Rhat", pr[8]))
+expect_true(grepl("n.eff", pr[8]))
+expect_true(grepl(" 51.876 ", pr[9]))
+expect_true(grepl("Successful", pr[30]))
+expect_true(grepl("overlap0", pr[34]))
+expect_true(grepl("DIC", pr[38]))
+
+# Rounded
+pr2 <- capture.output(print(out, digits=2))
+expect_equal(substr(pr2[9], 10,15), "51.88 ")
+
+# With bad Rhat
+out2 <- out
+out2$summary[1,"Rhat"] <- 1.2
+pr2 <- capture.output(print(out2))
+expect_true(grepl("convergence failure", pr2[30]))
+
+# With an NA rhat
+out2 <- out
+out2$summary[1,"Rhat"] <- NA
+pr2 <- capture.output(print(out2))
+expect_true(grepl("WARNING", pr2[30]))
+
+# With all NA rhats
+out2 <- out
+out2$summary[,"Rhat"] <- NA
+pr2 <- capture.output(print(out2))
+expect_true(grepl("WARNING", pr2[30]))
+
+# With 1 chain
+out2 <- out
+out2$mcmc.info$n.chains <- 1
+pr2 <- capture.output(print(out2))
+expect_false(grepl("Rhat", pr2[8]))
+expect_false(grepl("n.eff", pr2[8]))
+expect_true(grepl("overlap0", pr2[30]))
+
+# bugs.format
+out2 <- out
+out2$bugs.format <- TRUE
+pr2 <- capture.output(print(out2))
+expect_true(grepl("Bugs", pr2[1]))
+expect_true(grepl("and Rhat", pr2[27]))
diff --git a/inst/tinytest/test_process_output.R b/inst/tinytest/test_process_output.R
new file mode 100644
index 0000000..edc4de5
--- /dev/null
+++ b/inst/tinytest/test_process_output.R
@@ -0,0 +1,266 @@
+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, 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", 
+                            "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", DIC=TRUE, quiet=TRUE))
+expect_true(all(is.na(result)))
+
+#DIC=FALSE
+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"), 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))))
+expect_identical(rownames(out3$summary), c("beta", "sigma", "deviance"))
+
+# Check progress messages
+co <- capture.output(out <- process_output(samples, DIC=TRUE, quiet=FALSE))
+expect_identical(co, c("Calculating statistics....... ", "", "Done. "))
+                
+
+# Unexpected error happens during process_output-------------------------------
+
+# Here one of the arguments is missing
+expect_message(out_fail <- process_output(samples, quiet=TRUE))
+expect_true(is.null(out_fail)) # result is NULL
+
+#test that process_output matches old jagsUI process.output--------------------
+old_all <- readRDS("old_jagsUI_output.Rds")
+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
+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 calculation of Rhat is correct-------------------------------------
+samples <- readRDS('coda_samples.Rds')
+alpha <- samples[,"alpha"]
+expect_equal(calc_Rhat(alpha), 1.003831, tol=1e-4)
+expect_error(calc_Rhat(samples))
+expect_equal(calc_Rhat(alpha[1]), NA)
+alpha[[1]][1] <- Inf
+expect_equal(calc_Rhat(alpha), NA)
+
+
+# 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
+expect_message(out <- calc_param_stats(alpha_one))
+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, 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, 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, DIC=TRUE)))
+samp_inf <- samples
+samp_inf[[1]][1,ind] <- Inf
+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, DIC=TRUE)))
diff --git a/inst/tinytest/test_update.R b/inst/tinytest/test_update.R
new file mode 100644
index 0000000..60b05c3
--- /dev/null
+++ b/inst/tinytest/test_update.R
@@ -0,0 +1,69 @@
+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(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)
+
+out2 <- update(out, n.iter=100, n.thin=2, verbose=FALSE)
+ref <- readRDS("update_ref.Rds")
+expect_equal(out2$mcmc.info$n.iter, 200)
+expect_equal(out2$mcmc.info$n.samples, 150)
+expect_equal(nrow(out2$samples[[1]]), 50)
+
+# Remove time/date based elements
+out2$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out2[-c(17,19,21)], ref[-c(17,19,21)])
+
+# codaOnly---------------------------------------------------------------------
+out2 <- update(out, n.iter=100, n.thin=2, verbose=FALSE, codaOnly='mu')
+ref <- readRDS("update_ref_codaonly.Rds")
+
+out2$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out2[-c(17,19,21)], ref[-c(17,19,21)])
+
+# Different saved parameters---------------------------------------------------
+out2 <- update(out, n.iter=100, n.thin=2, verbose=FALSE, 
+               parameters.to.save=c('beta', 'alpha'))
+ref <- readRDS("update_ref_diffsaved.Rds")
+
+out2$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out2[-c(17,19,21)], ref[-c(17,19,21)])
+
+# DIC = FALSE------------------------------------------------------------------
+out2 <- update(out, n.iter=100, n.thin=2, verbose=FALSE, 
+               parameters.to.save=c('alpha'), DIC=FALSE)
+ref <- readRDS("update_ref_noDIC.Rds")
+
+expect_false(out2$calc.DIC)
+
+out2$mcmc.info$elapsed.mins <- ref$mcmc.inf$elapsed.mins
+expect_identical(out2[-c(15,17,19)], ref[-c(15,17,19)])
+
+# Check recovery after process_output errors-----------------------------------
+# Setting DIC to -999 forces process_output to error for testing
+expect_message(out2 <- update(out, n.iter=100, n.thin=2, verbose=FALSE, 
+               parameters.to.save=c('alpha'), DIC=-999))
+expect_inherits(out2, "jagsUIbasic")
+expect_equal(coda::varnames(out2$samples), c("alpha","deviance"))
+expect_equal(names(out2), c("samples", "model"))
diff --git a/inst/tinytest/update_ref.Rds b/inst/tinytest/update_ref.Rds
new file mode 100644
index 0000000..2563e73
--- /dev/null
+++ b/inst/tinytest/update_ref.Rds
diff --git a/inst/tinytest/update_ref_codaonly.Rds b/inst/tinytest/update_ref_codaonly.Rds
new file mode 100644
index 0000000..d382e06
--- /dev/null
+++ b/inst/tinytest/update_ref_codaonly.Rds
diff --git a/inst/tinytest/update_ref_diffsaved.Rds b/inst/tinytest/update_ref_diffsaved.Rds
new file mode 100644
index 0000000..d47393c
--- /dev/null
+++ b/inst/tinytest/update_ref_diffsaved.Rds
diff --git a/inst/tinytest/update_ref_noDIC.Rds b/inst/tinytest/update_ref_noDIC.Rds
new file mode 100644
index 0000000..fc16681
--- /dev/null
+++ b/inst/tinytest/update_ref_noDIC.Rds
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{
diff --git a/tests/tinytest.R b/tests/tinytest.R
new file mode 100644
index 0000000..6c6da60
--- /dev/null
+++ b/tests/tinytest.R
@@ -0,0 +1,4 @@
+if ( requireNamespace("tinytest", quietly=TRUE) ){
+	tinytest::test_package("jagsUI")
+}
+