#' 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.