#' @title Run MCMC
#' @description Conduct MCMC simulations using JAGS
#' @param jags.model specify which module to use
#' @param params define parameters to observe, Default: NULL
#' @param name.list list of names
#' @param data.list list of data
#' @param initial.list initial values for analysis, Default: list()
#' @param run.contrasts logical, indicating whether or not to run contrasts, Default: FALSE
#' @param use.contrast choose from "between", "within" and "mixed". Between compare groups at different conditions. Within compare a group at different conditions. Mixed compute all comparisons, Default: "between",
#' @param contrasts define contrasts to use for analysis (defaults to all) , Default: NULL
#' @param custom.contrast define contrasts for custom models , Default: NULL
#' @param run.ppp logical, indicating whether or not to conduct ppp analysis, Default: FALSE
#' @param k.ppp run ppp for every kth length of MCMC chains, Default: 10
#' @param n.data sample size for each parameter
#' @param credible.region summarize uncertainty by defining a region of most credible values (e.g., 95 percent of the distribution), Default: 0.95
#' @param save.data logical, indicating whether or not to save data, Default: FALSE
#' @param ROPE define range for region of practical equivalence (e.g., c(-0.05 , 0.05), Default: NULL
#' @param merge.MCMC logical, indicating whether or not to merge MCMC chains, Default: FALSE
#' @param run.diag logical, indicating whether or not to run diagnostics, Default: FALSE
#' @param param.diag define parameters to use for diagnostics, default equals all parameters, Default: NULL
#' @param sep symbol to separate data (e.g., comma-delimited), Default: ','
#' @param monochrome logical, indicating whether or not to use monochrome colors, else use \link[bfw]{DistinctColors}, Default: TRUE
#' @param plot.colors range of color to use, Default: c("#495054", "#e3e8ea")
#' @param graphic.type type of graphics to use (e.g., pdf, png, ps), Default: 'pdf'
#' @param plot.size size of plot, Default: '15,10'
#' @param scaling scale size of plot, Default: 100
#' @param plot.aspect aspect of plot, Default: NULL
#' @param vector.graphic logical, indicating whether or not visualizations should be vector or raster graphics, Default: FALSE
#' @param point.size point size used for visualizations, Default: 12
#' @param font.type font type used for visualizations, Default: 'serif'
#' @param one.file logical, indicating whether or not visualizations should be placed in one or several files, Default: TRUE
#' @param ppi define pixel per inch used for visualizations, Default: 300
#' @param units define unit of length used for visualizations, Default: 'in'
#' @param layout define a layout size for visualizations, Default: 'a4'
#' @param layout.inverse logical, indicating whether or not to inverse layout (e.g., landscape) , Default: FALSE
#' @param ... further arguments passed to or from other methods
#' @return list containing MCMC chains , MCMC chains as matrix , summary of MCMC, list of name used, list of data, the jags model, running time of analysis and names of saved files
#' @seealso
#' \code{\link[runjags]{runjags.options}},\code{\link[runjags]{run.jags}}
#' \code{\link[parallel]{detectCores}}
#' \code{\link[coda]{as.mcmc.list}},\code{\link[coda]{varnames}}
#' \code{\link[plyr]{rbind.fill}}
#' \code{\link[stats]{cor}},\code{\link[stats]{cov}},\code{\link[stats]{sd}}
#' \code{\link[MASS]{mvrnorm}}
#' \code{\link[utils]{write.table}}
#' @rdname RunMCMC
#' @export
#' @importFrom runjags runjags.options run.jags
#' @importFrom parallel detectCores
#' @importFrom coda as.mcmc.list varnames
#' @importFrom stats cor cov sd
#' @importFrom MASS mvrnorm
#' @importFrom utils write.table
RunMCMC <- function(jags.model,
params = NULL,
name.list,
data.list,
initial.list = list(),
run.contrasts = FALSE,
use.contrast = "between",
contrasts = NULL,
custom.contrast = NULL,
run.ppp = FALSE,
k.ppp = 10,
n.data,
credible.region = 0.95,
save.data = FALSE,
ROPE = NULL,
merge.MCMC = FALSE,
run.diag = FALSE,
param.diag = NULL,
sep = ",",
monochrome = TRUE,
plot.colors = c("#495054", "#e3e8ea"),
graphic.type = "pdf",
plot.size = "15,10",
scaling = 100,
plot.aspect = NULL,
vector.graphic = FALSE,
point.size = 12,
font.type = "serif",
one.file = TRUE,
ppi = 300,
units = "in",
layout = "a4",
layout.inverse = FALSE,
...
) {
# Name list
project.name <- name.list$project.name
project.data <- name.list$project.data
project.dir <- name.list$project.dir
model.type <- name.list$model.type
model.name <- name.list$model.name
job.title <- name.list$job.title
job.names <- name.list$job.names
job.group <- name.list$job.group
saved.steps = name.list$saved.steps
thinned.steps = name.list$thinned.steps
adapt.steps = name.list$adapt.steps
burnin.steps = name.list$burnin.steps
jags.seed <- name.list$jags.seed
jags.method <- name.list$jags.method
jags.chains <- name.list$jags.chains
# Add jags seed to initial list
initial.list <- c(initial.list, .RNG.seed = jags.seed)
# Number of samples
n.samples <- ceiling(saved.steps / jags.chains)
# Display start time
start.time <- Sys.time()
cat(format(start.time,"\nStarted at %d.%m.%Y - %H:%M:%S\n"))
# conduct JAGS
if ( is.null(merge.MCMC) ) {
# Initializing model
jags.data <- runjags::run.jags( method = jags.method ,
model= jags.model ,
monitor = params ,
data = data.list ,
inits = initial.list ,
n.chains = jags.chains ,
adapt = adapt.steps ,
burnin = burnin.steps ,
sample = n.samples ,
thin = thinned.steps ,
summarise = FALSE ,
plots = FALSE
)
# Generate random samples from the posterior distribution of monitored parameters
# Resulting data.MCMC object has these indices:
# data.MCMC[[ chainIdx ]][ stepIdx , paramIdx ]
data.MCMC <- coda::as.mcmc.list(jags.data)
# inspect merged MCMC chains
} else {
data.MCMC <- merge.MCMC
}
# Clean up working memory
RemoveGarbage("jags.data")
# Treat results as matrix for further examination
matrix.MCMC <- as.matrix(data.MCMC)
# Append bracket (ie., [1] ) for naming purposes
colnames(matrix.MCMC) <- unlist(lapply(colnames(matrix.MCMC), function (x) if(regexpr('\\[', x)[1]<0) paste0(x,"[1]") else x ))
coda::varnames(data.MCMC) <- unlist(lapply(coda::varnames(data.MCMC), function (x) if(regexpr('\\[', x)[1]<0) paste0(x,"[1]") else x ))
# Diagnostics
if (run.diag) {
cat("\nConducting diagnostics. Please wait.\n")
diag.start.time <- Sys.time()
if (length(param.diag)) {
param.diag <- grep(paste(TrimSplit(param.diag),collapse="|"),
coda::varnames(data.MCMC), value = TRUE)
} else {
param.diag <- coda::varnames(data.MCMC)
}
diag.length <- length(param.diag)
diag.plots <- lapply(1:diag.length, function(i) {
x <- param.diag[i]
if (stats::sd(matrix.MCMC[,x]) > 0) {
diag <- DiagMCMC(data.MCMC = data.MCMC,
par.name = x,
job.names = job.names,
job.group = job.group,
credible.region = credible.region
)
} else {
diag <- (paste0(x , " appears to be a constant. No available diagnostics."))
}
ETA(diag.start.time , i , diag.length)
return (diag)
} )
# Adding names to diagnostics
names(diag.plots) <- param.diag
if (save.data) {
cat("\nAdjusting and saving diagnostic plots. Please wait.\n")
# Save plots as PowerPoint, Default is raster graphics.
## Change vector.graphic to TRUE if needed (not recommended)
ParsePlot(diag.plots[lapply(diag.plots,length)>1],
project.dir = paste0(project.dir,"Diagnostics/"),
project.name = project.name,
graphic.type = graphic.type,
plot.size = plot.size,
save.data = save.data,
vector.graphic = vector.graphic,
point.size = point.size,
font.type = font.type,
one.file = one.file,
layout = layout,
)
}
}
# Add contrasts for selected models
if (run.contrasts) {
q.levels <- data.list$q.levels
defined.contrast <- if (length(custom.contrast)) tolower(custom.contrast) else tolower(model.type)
cat("\nComputing contrasts. Please wait.\n")
# Createsum to zero contrasts for metric and nominal models
if (defined.contrast == "metric" | defined.contrast == "nominal") {
sum.zero <- SumToZero(q.levels, matrix.MCMC, contrasts)
}
# Add odds (and odds-ratios/cohen's d) for nominal models
if (defined.contrast == "nominal") {
# Create expected and observed nominal and proportions data
expected.observed <- MatrixCombn(matrix.MCMC,
"o,o,e,e",
"NULL,p,NULL,p",
q.levels,
rm.last = FALSE,
row.means = FALSE)
# Remove expected and observed columns from MCMC matrix
matrix.MCMC <- matrix.MCMC[ , !colnames(matrix.MCMC) %in% colnames(expected.observed)]
contrast.type <- c("b","o")
contrast.data <- list(sum.zero, expected.observed)
}
# Add effect size cohen's d for metric model
if (defined.contrast == "metric") {
# Create mean difference data
mean.diff <- MatrixCombn(matrix.MCMC, "m,s", q.levels = q.levels, rm.last = FALSE)
# Remove mean difference columns from MCMC matrix
matrix.MCMC <- matrix.MCMC[ , !colnames(matrix.MCMC) %in% colnames(mean.diff)]
contrast.type <- c("b","m")
contrast.data <- list(sum.zero, mean.diff)
}
# Run contrasts
contrasts <- do.call(cbind, lapply(1:length(contrast.type), function (i) {
RunContrasts(contrast.type[[i]],
q.levels,
use.contrast,
contrasts,
contrast.data[[i]],
job.names)
}))
}
# Add R^2 if regression
if (model.type == "Regression") {
x <- data.list$x
n.x <- data.list$n.x
y <- data.list$y
q <- data.list$q
n <- data.list$n
r.squared <- do.call(cbind,lapply(1:q, function (i) {
zbeta <- matrix.MCMC[,grep(paste0("^zbeta$|^zbeta\\[",i),colnames(matrix.MCMC))]
# Compute R^2
cor.data <- stats::cor( x[,1:n.x[i]] , y )
r.squared <- zbeta %*% matrix( cor.data , ncol=1 )
# Compute robust R^2
cor.data <- stats::cor( x[,1:n.x[i]] , y , method="spearman")
rob.r.squared <- zbeta %*% matrix( cor.data , ncol=1 )
# Compute adjusted R^2
adjusted.r <- as.matrix(1 - (1-r.squared) * ( (n-1) / (n -(q+1)) ))
# Compute adjusted robust R^2
adjusted.rob.r <- as.matrix(1 - (1-rob.r.squared) * ( (n-1) / (n -(q+1)) ))
# Add names
colnames(r.squared) <- if (q>1) sprintf("R^2 (block: %s)",i) else "R^2"
colnames(rob.r.squared) <- if (q>1) sprintf("Robust R^2 (block: %s)",i) else "Robust R^2"
colnames(adjusted.r) <- if (q>1) sprintf("Adjusted R^2 (block: %s)",i) else "Adjusted R^2"
colnames(adjusted.rob.r) <- if (q>1) sprintf("Adjusted Robust R^2 (block: %s)",i) else "Adjusted Robust R^2"
# Return as matrix
cbind(r.squared , rob.r.squared , adjusted.r , adjusted.rob.r)
}))
}
# Add PPP if SEM/CFA
if (run.ppp) {
y <- data.list$y
lat <- data.list$lat
factor.seq <- data.list$factor.seq
n <- data.list$n
ppp.start.time <- Sys.time()
pppv <- 0
jags.ppp <- sample(1:nrow(matrix.MCMC), nrow(matrix.MCMC)/k.ppp)
cov.mat <- stats::cov(y)
cat("\nComputing PPP-value, please wait (it may take some time).\n")
PPP <- lapply(1:length(jags.ppp), function (i) {
x <- matrix.MCMC[jags.ppp[i],]
# Epsilon/Error variance matrix
eps.pos <- grep("error", names(x))
eps.length <- length(eps.pos)
eps.matrix <- diag(x[eps.pos],eps.length)
# Rho/Latent covariance matrix
rho.pos <- grep("cov", names(x))
rho.mat <- matrix(x[rho.pos], nrow=lat, ncol=lat, byrow=TRUE)
# Lambda/Loading matrix
lam.pos <- grep("lam", names(x))
lam.length <- length(lam.pos)
lambda <- do.call( rbind, lapply(1:lam.length, function (j) c(factor.seq[j], x[lam.pos][j]) ) )
lam.matrix <- matrix( unlist ( lapply(1:lat, function (j) {
if (j < max(lat)) {
c ( lambda[ lambda[,1] == j , 2], rep(0,lam.length) )
} else {
lambda[ lambda[,1] == j , 2]
}
} ) ), lam.length, lat)
# Calculate predicted covariance matrix
pred.sigma <- lam.matrix %*% rho.mat %*% t(lam.matrix) + eps.matrix
# Chi-square discrepancy (predicted data from model vs. covariance)
pred.fit <- (n - 1) * (log(det(pred.sigma)) + sum(diag(solve(pred.sigma) %*% cov.mat )) - log(det(cov.mat)) - eps.length)
# Simulate covariances from model
sim.sigma <- cov ( MASS::mvrnorm(n=n, mu=rep(0,eps.length), Sigma=pred.sigma, empirical=FALSE) )
# Chi-square discrepancy (simulated data from model vs. simulated covariance)
sim.fit <- (n - 1) * (log(det(pred.sigma)) + sum(diag(solve(pred.sigma) %*% sim.sigma)) - log(det(sim.sigma)) - eps.length)
# Compute PPP-value
pppv <<- if (pred.fit < sim.fit) pppv + 1 else if (pred.fit == sim.fit) pppv + 0.5 else pppv
PPP <- as.numeric(pppv / i)
ETA(ppp.start.time , i , length(jags.ppp))
# Create matrix with chi-square, discrepancy between predicted and simulated data and PPP
PPP <- as.data.frame(t(
c("Fit (Predicted)" = pred.fit,
"Fit (Simulated)" = sim.fit,
"Fit (Discrepancy)" = (pred.fit-sim.fit),
"PPP" = PPP)
))
return (PPP)
} )
if (requireNamespace("plyr", quietly = TRUE)) {
PPP <- plyr::rbind.fill(PPP)
} else {
PPP <- do.call(rbind,PPP)
}
}
# create list of matrices to summarize
list.MCMC <- list(
matrix = matrix.MCMC,
PPP = if (exists("PPP")) PPP,
r.squared = if (exists("r.squared")) r.squared,
sum.zero = if (exists("sum.zero")) sum.zero,
count.data = if (exists("expected.observed")) expected.observed,
mean.difference = if (exists("mean.diff")) mean.diff,
contrasts = if (exists("contrasts")) contrasts
)
# Remove empty plots
list.MCMC <- Filter(length, list.MCMC)
summary.MCMC <- do.call(rbind,lapply(1:length(list.MCMC), function (k) {
# Find params from MCMC list
params <- colnames(list.MCMC[[k]])
# Create final posterior parameter indicies
summary.start.time <- Sys.time()
summary.cat <- sprintf("\nSummarizing data for each parameter in %s. Please wait.\n" ,
gsub("\\."," ",names(list.MCMC)[[k]]) )
cat(summary.cat)
do.call(rbind,lapply(1:length(params), function(i) {
summary <- SumMCMC( par = list.MCMC[[k]][, params[i]] ,
par.names = params[i],
job.names = job.names,
job.group = job.group,
n.data = n.data,
credible.region = credible.region,
ROPE = ROPE
)
ETA(summary.start.time , i , length(params))
return (summary)
}))
}))
# Sort summary by rownames
if (nrow(summary.MCMC)>1) summary.MCMC <- summary.MCMC[order(rownames(summary.MCMC)) ,]
# Display completion and running time
stop.time <- Sys.time()
total.time <- capture.output(difftime(stop.time, start.time))
cat(format(stop.time,"\nCompleted at %d.%m.%Y - %H:%M:%S\n"),
gsub("Time difference of","Running time:",total.time),"\n", sep="")
# Create MCMC and summary list (without chain information)
final.MCMC <- list( raw.MCMC = data.MCMC,
matrix.MCMC = do.call(cbind,list.MCMC),
summary.MCMC = summary.MCMC,
name.list = name.list,
initial.list = initial.list,
data.list = data.list,
jags.model = jags.model,
run.time = c(start.time,stop.time)
)
# If applicable, add diag plots
if (run.diag & !save.data) final.MCMC <- c(final.MCMC, diag = list(diag.plots))
# Save data if requested
if (save.data) {
cat("\nSaving data. Please wait.\n")
# If requested, save summary as csv
data.file.name <- RemoveSpaces(paste0(project.dir,project.name,".csv"))
utils::write.table(summary.MCMC, data.file.name, sep = sep)
# Save final MCMC
MCMC.file.name <- paste0(project.dir,"MCMC/",project.name,".rds")
saveRDS( final.MCMC , file = MCMC.file.name, compress="bzip2")
# Create meta data
meta.final.MCMC <- list(summary.MCMC = summary.MCMC,
name.list = name.list,
initial.list = initial.list,
data.list = data.list,
jags.model = jags.model,
run.time = c(start.time,stop.time)
)
# Save meta data
MCMC.file.name <- sub('-', '-META-', MCMC.file.name)
saveRDS( meta.final.MCMC , file = MCMC.file.name, compress="bzip2")
# Append to final MCMC list
final.MCMC <- c(final.MCMC,
data.file.name = data.file.name,
MCMC.file.name = MCMC.file.name)
}
# Clear up working memory
RemoveGarbage("data.MCMC,list.MCMC")
return (final.MCMC)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.