#' Summarize MCMC Output (for one simulation run)
#'
#'A function that intakes MCMC samples and outputs a tibble that consists of posterior estimates with their fixed values at start in one tibble.
#'Other variable included in tibble is bias calculated for each parameter. Included parameters in the tibble are all global parameters, country-year estimates of sens, spec, vradj and gamma.truemat.vr.
#' @param simsetting Describes sim setting
#' @param runname Describes run setting
#' @param simnum Simulation number i out of Nsim
#' @param globalvars Vector of global variable names
#' @param output.dir Output directory name
#'
#' @return Posterior estimates of global variables and country-year parameters (sens, spec, vradj, gammatruematvr), calculate bias, and set fixed values from simulation code.
#' @export
#'
#' @examples
#' summarize_simulation(simsetting = "test", runname = "test", simnum = 1, globalvars = c("sensworld"), output.dir = output.dir)
summarize_simulation <-
function(simsetting,
runname,
simnum,
globalvars = globalvars,
output.dir) {
tib <-
readRDS(paste(output.dir, 'mcmctib_', simnum, '.RDS', sep = ""))
tidydat <- tib %>%
select(
contains("sens"),
contains("spec"),
contains("vradj.ct"),
contains("gamma.truematvr.ct")
)
if (ncol(tidydat) < 62) { #for testing purposes
#only and all selected parameter should be in this set
stop(
paste0(
"At least oneparameter of interest is missing. .\n",
"Check to make sure mcmctib_ is correct"
)
)
call. = FALSE
} else {
tidydat <- tidydat %>%
gather(key = "parameter", value = "value") %>%
group_by(parameter) %>%
summarize(
lower = quantile(value, .025),
median = quantile(value, .5),
upper = quantile(value, .975)
) %>%
mutate(sim = simnum, runsetting = runname
)
#Pulls in simlist with the set values, but this is a list so much call by name
truthdat <-
readRDS(paste(output.dir, "simlist_iter", simnum, ".RDS", sep = ""))
global_settings <-
readRDS(paste(output.dir, "global_settings.RDS", sep = ""))
globalvars <-
c("sensworld",
"specworld",
"rho.alphabeta",
"sigma.alpha",
"sigma.beta")
C <- truthdat$C
Nyears <- truthdat$Nyears
tidydat[["truth"]] <- NA
for (c in 1:C) {
for (t in 1:Nyears) {
tidydat[which(tidydat$parameter == paste("gamma.truematvr.ct[", c, ",", t , "]", sep = "")), "truth"] <- truthdat$GAMMATRUE.ct[c, t]
tidydat[which(tidydat$parameter == paste("sens.ct[", c, ",", t , "]" , sep = "")), "truth"] <- truthdat$SENS.ct[c, t]
tidydat[which(tidydat$parameter == paste("spec.ct[", c, ",", t , "]", sep = "")), "truth"] <- truthdat$SPEC[c, t]
tidydat[which(tidydat$parameter == paste("vradj.ct[", c, ",", t , "]", sep = "")), "truth"] <- truthdat$VRADJ.ct[c, t]
}
}
for (g in 1:length(globalvars)) {
tidydat[which(tidydat$parameter == globalvars[g]), "truth"] <-
global_settings[[globalvars[g]]]
}
tidydat <- tidydat %>%
mutate(bias = (truth - median),
coverage = NA)
for (i in 1:nrow(tidydat)) {
if ((tidydat$lower[i] < tidydat$truth[i]) == TRUE &&
(tidydat$upper[i] > tidydat$truth[i]) == TRUE) {
tidydat$coverage[i] <- 1
} else {
tidydat$coverage[i] <- 0
}
}
saveRDS(tidydat, file = paste(output.dir, "sumtib_", simnum, '.RDS', sep = ""))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.