R/run_mcmc.R

In Enpeterson/outputsim: Output Simulation Results

#' Run Bayesian Model for Se and Sp
#'
#'A function that takes in settings for a JAGS run, and parameter settings. Function writes model based on which parameters we are choosing to model, runs JAGS and saves tibbles of MCMC values per simulations.
#'The MCMC tibbles are saved to the output folder by simulation number out of Nsim simulations.
#'Last, there are traceplots and convergence information saved to output folder as part of coda package output.
#' @param runname Extension to output folder that describes run setting, ie alldata.
#' @param Nchains Number of chains ran by JAGS.
#' @param Niter Number of iterations per chain ran by JAGS.
#' @param i #Simulation number out of Nsim.
#' @param modelonegammatrue Indicator (TRUE or FALSE) if the model includes one parameter for gammatruematvr.
#' @param modelgammatruect Indicator (TRUE or FALSE) if the model includes unique country-year parametes for gammatruematvr.
#' @param simsetting Describes simulation setting.
#'
#' @return Takes in jags inputs from run_simulation. Writes and runs specified model. Outputs a tibble of MCMC samples (for all parnames listed below) for the given simulation.
#' @export
#'
#' @examples
#' run_mcmc(runname = "test", Nchains = 1, Niter = 10, i=1, modelonegammatrue = T, modelgammatruect = F, simsetting = "test")
#'
run_mcmc <- function(runname , Nchains, Niter,  i , modelonegammatrue , modelgammatruect, simsetting){

jagsdata <- readRDS(paste(output.dir, "simlist_iter", i, ".RDS", sep=""))


rnorm(1)
set.seed(1234)
parnames <-  c("rhovr.ct",
               "gamma_true",
                    "gamma.truematvr.ct",
                    "rho.ctb",
                    "sens.ct",
                    "spec.ct",
                    "lambda.ctk",
                    "vradj.ct",
                    "sensworld",
                    "specworld",
                    "rho.alphabeta",
                    "sigma.alpha",
                    "sigma.beta"
                    #"rrho.alphabeta",
                    #"rsigma.alpha",
                    #"rsigma.beta"
                    )


      write_model_to_run(output.dir = output.dir,
                         modelonegammatrue= modelonegammatrue,
                         modelgammatruect = modelgammatruect,
                         simsetting=simsetting) # Writes model that I use to model simulated data

      model <- jags.model( file = file.path(paste(output.dir, sep=""), "run_data_model.txt"),
                                           data = jagsdata,
                                           n.chains = Nchains,
                                          n.adapt = 3000)
      mcmc_samples <- coda.samples(model,
                              parnames,
                              n.iter = Niter)



      pdf(paste(output.dir, "traceplots_iter", i, ".pdf", sep=""))
      plot(mcmc_samples)
      dev.off()

      # pdf(paste(output.dir, "shrinkageplots_iter", i, ".pdf", sep=""))
      # gelman.plot(mcmc_samples)
      # dev.off()

     #mcmc.array <- as.array(mcmc_samples)
     #saveRDS(mcmc.array, paste(output.dir, "mcmc_array_iter",i, ".RDS", sep=""))

     mcmc.list <- combine.mcmc(mcmc_samples)
     mcmc_tib <- as.tibble(mcmc.list)
     mcmc_tib$runsetting <- simsetting
     mcmc_tib$sim <- i
     saveRDS(mcmc_tib, paste(output.dir, "mcmctib_", i, ".RDS", sep=""))



}#end function