# Clear workspace
rm(list = ls())
# Load packages
suppressPackageStartupMessages({
library("Matrix")
library("MASS")
library("OpenMx")
library("metaSEM")
library("here")
library("parallel")
library("pbapply")
})
The function sim_wrapper() simulates a single meta-analysis given the set of provided arguments. All required functions for the simulation are specified as internal functions within sim_wrapper() and are indicated by a leading “.”.
The entire simulation of multiple meta-analyses (i.e., for multiple replications) given a set of arguments is invoked by sim_run() and is run in parallel on multiple cores.
The SEM model specifies either a simple mediation model with a single mediator (M), a multiple mediator model with two parallel mediators (M, W), or a sequential mediator model with two successive mediators (M, W). All models include a single independent variable (X) and a single outcome (Y).
#
# Simulate a single meta-analysis
#
# @param rid A numeric identifier for the replication
# @param MONX A standardized regression weight of M on X
# @param YONX A standardized regression weight of Y on X
# @param YONM A standardized regression weight of Y on M (optional)
# @param WONX A standardized regression weight of W on X (optional)
# @param YONW A standardized regression weight of Y on W (optional)
# @param WONM A standardized regression weight of W on M (optional)
# @param rxx A vector with a minimum and maximum reliability
# for X (optional)
# @param rmm A vector with a minimum and maximum reliability
# for M (optional)
# @param rww A vector with a minimum and maximum reliability
# for W (optional)
# @param ryy A vector with a minimum and maximum reliability
# for Y (optional)
# @param tau A random effect (SD) for the population
# correlations (optional)
# Note: The same random effect is used for
# all correlations. Covariances between
# random effects are set to 0.
# @param k A vector with a minimum and maximum number of samples
# in the meta-analysis (optional)
# @param n A vector with a minimum and maximum sample size
# for each sample in the meta-analysis (optional)
# @param tau_max Maximum random effect (SD) (optional)
# @return A vector with the simulated SEM parameters,
# error codes, and summaries of the used
# settings
# @require MASS, Matrix, OpenMx, metaSEM
#
# If optional parameters are not set, they are randomly drawn from
# from reasonable distributions.
#
<- function(rid, MONX, YONX, YONM = NULL,
sim_wrapper WONX = NULL, YONW = NULL, WONM = NULL,
rxx = NULL, rmm = NULL, rww = NULL, ryy = NULL,
tau = NULL, k = NULL, n = NULL,
tau_max = .50) {
#
# Generate the implied correlation matrix from
# standardized regression weights of a mediation model
#
# Assumes a simple mediation model, unless WONX, YONW, or
# YONW are specified. If WONX and YONW are given, a two
# variable parallel mediation is assumed. If YONW and WONM
# are specified, a two variable sequential mediation is assumed
#
# @param MONX A standardized regression weight of M on X
# @param YONX A standardized regression weight of Y on X
# @param YONM A standardized regression weight of Y on M (optional)
# @param WONX A standardized regression weight of W on X (optional)
# @param YONW A standardized regression weight of Y on W (optional)
# @param WONM A standardized regression weight of W on M (optional)
# @return A 3x3 or 4x4 correlation matrix
# @require metaSEM
#
<- function(MONX, YONX, YONM = NULL,
.beta_to_cor WONX = NULL, YONW = NULL,
WONM = NULL) {
# Single mediator model
if (is.null(WONX) & is.null(YONW) & is.null(WONM) &
!is.null(YONM)) {
<- matrix(c(0, 0, 0,
Amat 0, 0,
MONX, 0), 3, 3, byrow = TRUE)
YONX, YONM, colnames(Amat) <- rownames(Amat) <- c("X", "M", "Y")
<- matrix(c(1, 0, 0,
Smat 0, 1 - MONX^2, 0,
0, 0, 1 - YONM^2 - YONX^2 - 2 * MONX * YONM * YONX),
3, 3, byrow = TRUE)
<- diag(3)
Imat
# Two parallel mediators model
else if (!is.null(WONX) & !is.null(YONM) & !is.null(YONW) &
} is.null(WONM)) {
<- matrix(c(0, 0, 0, 0,
Amat 0, 0, 0,
MONX, 0, 0, 0,
WONX, 0), 4, 4, byrow = TRUE)
YONX, YONM, YONW, colnames(Amat) <- rownames(Amat) <- c("X", "M", "W", "Y")
<- matrix(c(1, 0, 0, 0,
Smat 0, 1 - MONX^2, 0, 0,
0, 0, 1 - WONX^2, 0,
0, 0, 0, 1 - YONM^2 - YONW^2 - YONX^2 -
2 * MONX * YONM * YONX -
2 * WONX * YONW * YONX -
2 * MONX * WONX * YONM * YONW),
4, 4, byrow = TRUE)
<- diag(4)
Imat
# Two sequential mediators model
else if (is.null(WONX) & is.null(YONM) &
} !is.null(YONW) & !is.null(WONM)) {
<- matrix(c(0, 0, 0, 0,
Amat 0, 0, 0,
MONX, 0, WONM, 0, 0,
0, YONW, 0), 4, 4, byrow = TRUE)
YONX, colnames(Amat) <- rownames(Amat) <- c("X", "M", "W", "Y")
<- matrix(c(1, 0, 0, 0,
Smat 0, 1 - MONX^2, 0, 0,
0, 0, 1 - WONM^2, 0,
0, 0, 0, 1 - YONW^2 - YONX^2 -
2 * MONX * WONM * YONW * YONX),
4, 4, byrow = TRUE)
<- diag(4)
Imat
else {
}
stop("Please specify a valid set of arguments!")
}
<- solve(Imat - Amat) %*% Smat %*% t(solve(Imat - Amat))
rmat if (!metaSEM::is.pd(rmat))
stop("Correlation matrix is not positive definite. Check your arguments.")
return(rmat)
}
#
# Generate random heterogeneity estimates (SD)
#
# @param tau A fixed heterogeneity estimated (optional)
# @param times The number of heterogeneity estimates (optional)
# @param tau_max A maximum heterogeneity estimate (optional)
# @return A numeric vector
#
<- function(tau = NULL, times = 1, tau_max = .50) {
.get_tau
if (!is.null(tau)) return(rep(tau[1], times))
<- abs(rnorm(times, 0, .20))
tau repeat {
> tau_max] <- tau[tau > tau_max] / 2
tau[tau if (all(tau <= tau_max)) break
}return(round(tau, 2))
}
#
# Generate random sample sizes
#
# @param n A vector with a minimum and maximum sample size (optional)
# @param times The number of sample sizes (optional)
# @return A numeric vector
#
<- function(n = NULL, times = 1) {
.get_n
if (is.null(n)) n <- c(50, 900)
if (length(n) == 1) n <- c(n, n)
<- round(n)
n <- rgamma(times, 1.3, 1 / 250) + n[1]
ns repeat {
> n[2]] <- round(ns[ns > n[2]] / 2)
ns[ns < n[1]] <- n[1]
ns[ns if (all(ns <= n[2])) break
}return(round(ns))
}
#
# Generate random numbers of samples
#
# @param k A vector with a minimum and maximum number of samples (optional)
# @param times The number of number of samples (optional)
# @return A numeric vector
#
<- function(k = NULL, times = 1) {
.get_k
if (is.null(k)) k <- c(5, 50)
if (length(k) == 1) k <- c(k, k)
<- round(k)
k return(round(runif(times, k[1], k[2])))
}
#
# Generate random reliability estimates
#
# @param rxx A vector with a minimum and maximum reliability (optional)
# @param times The number of reliabilities (optional)
# @return A numeric vector
#
<- function(rxx = NULL, times = 1) {
.get_rel
if (is.null(rxx)) rxx <- c(.65, .95)
if (length(rxx) == 1) rxx <- c(rxx, rxx)
return(round(runif(times, rxx[1], rxx[2]), 2))
}
#
# Generate a RAM specification for a mediation model
# for metaSEM::tssem2() and metaSEM::osmasem()
#
# @param type Specify type of mediation model:
# - "simple" for simple mediation model
# - "parallel" for parallel mediation model
# - "sequential" for sequential mediation model
# @return A list with the RAM specification
# @require metaSEM
#
<- function(type = "simple") {
.get_model
if (type == "sequential") {
return(metaSEM::lavaan2RAM(
"M ~ X; W ~ M; Y ~ W + X; X ~~ 1*X",
obs.variables = c("X", "M", "W", "Y")
))
else if (type == "parallel") {
}
return(metaSEM::lavaan2RAM(
"M ~ X; W ~ X; Y ~ M + W + X; X ~~ 1*X",
obs.variables = c("X", "M", "W", "Y")
))
else if (type == "simple") {
}
return(metaSEM::lavaan2RAM(
"M ~ X; Y ~ M + X; X ~~ 1*X",
obs.variables = c("X", "M", "Y")
))
else {
}
stop("Unknown model!")
}
}
#
# Define OpenMx matrices for indirect effects
#
# @param model A list with the RAM specification of the mediation model
# as returned by .get_model()
# @return A list with OpenMx matrices
# @require OpenMx
#
<- function(model) {
.get_INDmatrix
# Simple mediation
if (ncol(model$A) == 3) {
<- list(
out IND = OpenMx::mxAlgebra(Amatrix[2, 1] * Amatrix[3, 2], name = "IND")
)
# Sequential mediation
else if (ncol(model$A) == 4 & model$A["W", "X"] == "0") {
}
<- list(
out IND = OpenMx::mxAlgebra(Amatrix[2, 1] * Amatrix[3, 2] * Amatrix[4, 3],
name = "IND")
)
# Parallel mediation
else {
}
<- list(
out IND1 = OpenMx::mxAlgebra(Amatrix[2, 1] * Amatrix[4, 2],
name = "IND1"),
IND2 = OpenMx::mxAlgebra(Amatrix[3, 1] * Amatrix[4, 3],
name = "IND2"),
IND = OpenMx::mxAlgebra(Amatrix[2, 1] * Amatrix[4, 2] +
3, 1] * Amatrix[4, 3],
Amatrix[name = "IND")
)
}
return(out)
}
#
# Draw a random population correlation matrix
# given an average population correlation matrix
# and a random effect
#
# Note: A non-positive definite correlation
# matrix is discarded and replaced by
# a valid correlation matrix.
#
# @param avg_mat The average population correlation
# matrix
# @param tau The random effect (SD)
# @return The population correlation matrix
# @require MASS, metaSEM, OpenMx
#
<- function(avg_mat, tau) {
.sim_pop_mat
<- (nrow(avg_mat)^2 - nrow(avg_mat)) / 2
p <- diag(tau^2, p, p)
Sigma <- OpenMx::vechs(avg_mat)
mu repeat {
<- MASS::mvrnorm(1, mu = mu, Sigma = Sigma)
rho_vec <- metaSEM::vec2symMat(rho_vec, diag = FALSE)
rho_mat if (metaSEM::is.pd(rho_mat)) break # & all(abs(rho_vec) < .95)
}return(rho_mat)
}
#
# (Dis)Attentuate correlations / (co)variances
#
# @param r A matrix of correlations / (co)variances
# @param rxx A vector of reliability estimates for X
# @param rmm A vector of reliability estimates for M
# @param rww A vector of reliability estimates for W (optional)
# @param ryy A vector of reliability estimates for Y
# @param adj A list with logical indicators for each variable
# showing whether to adjust it for unreliability (optional)
# @param invert A logical value indicating whether to
# disattenuate (FALSE) or attenuate (TRUE)
# the correlations
# @param return_att Return the attenuation factors instead of the
# correlations
# @param vcov A logical indicating whether to disattenuate the
# correlations (FALSE) or their (co)variances (TRUE)
# @return A matrix of (dis)attenuated correlations / (co)variances
#
<- function(x, rxx, rmm, ryy, rww = NULL,
.disattenuate adj = NULL, invert = FALSE,
return_att = FALSE, vcov = FALSE) {
# Determine which correlations to adjust for unreliability
if (is.null(adj) | !is.list(adj)) adj <- list()
for (i in c("rxx", "rmm", "rww", "ryy")) {
if (!(i %in% names(adj))) adj[[i]] <- TRUE
}
# Correlations / (co)variances
if (!is.matrix(x)) x <- rbind(x)
# Reliabilities
if ((!vcov && ncol(x) == 3) | (vcov && ncol(x) == 6)) {
<- cbind(rxx, rmm, ryy)
rel_vecs else {
} <- cbind(rxx, rmm, rww, ryy)
rel_vecs
}if (!adj$rxx) rel_vecs[, "rxx"] <- 1
if (!adj$rmm) rel_vecs[, "rmm"] <- 1
if (!adj$rww & ncol(rel_vecs) == 4) rel_vecs[, "rww"] <- 1
if (!adj$ryy) rel_vecs[, "ryy"] <- 1
# Adjustment factors
<- t(apply(rel_vecs, 1, \(y) {
adj_vecs <- OpenMx::vechs(cbind(y^.5) %*% y^.5)
a if (vcov) a <- OpenMx::vech(cbind(a) %*% a)
a
}))if (invert) adj_vecs <- 1 / adj_vecs # attenuate correlations
# Adjust correlations
if (!return_att) {
<- t(apply(x / adj_vecs, 1, \(y) {
out <- metaSEM::vec2symMat(y, diag = vcov)
y if (!metaSEM::is.pd(y)) y <- Matrix::nearPD(y, corr = !vcov)$mat
if (vcov) {
::vech(y)
OpenMxelse {
} ::vechs(y)
OpenMx
}
}))dimnames(out) <- dimnames(x)
# Return attenuation factor
else {
} <- adj_vecs
out
}
return(out)
}
<- function(...) { .disattenuate(invert = TRUE, ...) }
.attenuate
#
# Draw a random sample correlation matrix
# with attenuation given an average
# population correlation matrix and a
# random effect
#
# Note: A non-positive definite correlation
# matrix is discarded and replaced by
# a valid correlation matrix.
#
# @param avg_mat The average population correlation matrix
# @param tau A random effect (SD)
# @param n A vector with a minimum and maximum sample size (optional)
# @param rxx A vector with a minimum and maximum reliability
# for X (optional)
# @param rmm A vector with a minimum and maximum reliability
# for M (optional)
# @param rww A vector with a minimum and maximum reliability
# for W (optional)
# @param ryy A vector with a minimum and maximum reliability
# for Y (optional)
# @return A numeric vector with the lower half of
# the sample correlation matrix,
# reliabilities, and a sample size
# @require MASS, metaSEM, OpenMx
#
<- function(avg_mat, tau, n = NULL,
.sim_r rxx = NULL, rmm = NULL,
rww = NULL, ryy = NULL) {
# Draw population correlation matrix for present sample
<- .sim_pop_mat(avg_mat, tau)
rho_mat
# Get sample size
<- .get_n(n)
n
# Get reliabilities
<- .get_rel(rxx)
rxx <- .get_rel(rmm)
rmm <- .get_rel(rww)
rww <- .get_rel(ryy)
ryy
# Repeat draw for non-positive definite matrices
<- 0
iter repeat {
# Correlation for random sample
<- MASS::mvrnorm(n, mu = rep(0, nrow(rho_mat)), Sigma = rho_mat)
resp <- OpenMx::vechs(cor(resp))
r_vec if (nrow(rho_mat) == 3) {
names(r_vec) <- paste0("r_", c("XM", "XY", "MY"))
else if (nrow(rho_mat) == 4) {
} names(r_vec) <- paste0("r_", c("XM", "XW", "XY",
"MW", "MY", "WY"))
}<- .attenuate(r_vec, rxx = rxx, rmm = rmm,
r_vec ryy = ryy, rww = rww)
# Check if matrix is positive definite
<- metaSEM::vec2symMat(r_vec, diag = FALSE)
r_mat if (metaSEM::is.pd(r_mat)) break
<- iter + 1
iter
}
<- c(r_vec[1, ], n = n,
out rxx = rxx, rmm = rmm, ryy = ryy,
iter = iter)
if (nrow(rho_mat) > 3) out <- c(out, rww = rww)
return(out)
}
#
# Simulate correlations for multiple samples
# in a meta-analysis
#
# @param MONX A standardized regression weight of M on X
# @param YONX A standardized regression weight of Y on X
# @param YONM A standardized regression weight of Y on M (optional)
# @param WONX A standardized regression weight of W on X (optional)
# @param YONW A standardized regression weight of Y on W (optional)
# @param WONM A standardized regression weight of W on M (optional)
# @param k A vector with a minimum and maximum number of samples (optional)
# @param n A vector with a minimum and maximum sample size (optional)
# @param tau A random heterogeneity estimate (optional)
# @param rxx A vector with a minimum and maximum reliability
# for X (optional)
# @param rmm A vector with a minimum and maximum reliability
# for M (optional)
# @param rww A vector with a minimum and maximum reliability
# for W (optional)
# @param ryy A vector with a minimum and maximum reliability
# for Y (optional)
# @param tau_max Maximum random effet (SD) (optional)
# @return A matrix with correlations
#
<- function(MONX, YONX, YONM = NULL,
.sim_samples WONX = NULL, YONW = NULL, WONM = NULL,
k = NULL, n = NULL, tau = NULL,
rxx = NULL, rmm = NULL, rww = NULL,
ryy = NULL, tau_max = .50) {
# Number of samples
<- .get_k(k)
k
# Random effect for correlations
<- .get_tau(tau, tau_max = tau_max)
tau
# Average population correlation matrix
<- .beta_to_cor(MONX = MONX, YONX = YONX, YONM = YONM,
avg_mat WONX = WONX, YONW = YONW, WONM = WONM)
# Simulate correlations for samples
<- replicate(k, .sim_r(avg_mat = avg_mat, tau = tau, n = n,
out rxx = rxx, rmm = rmm, rww = rww,
ryy = ryy))
<- cbind(sid = seq_len(k), t(out), tau = tau)
out
return(out)
}
#
# First stage in two-stage MASEM
#
# The sample correlations are adjusted individually
# for unreliability before pooling them.
#
# @param obj A matrix with sample statistics for multiple samples
# as returned by .sim_samples()
# @param adj A list with logical indicators for each variable
# showing whether to adjust it for unreliability (optional)
# @return A list including the pooled correlation
# matrix, its variance-covariance matrix,
# the total sample size, the number of studies,
# and the mean reliabilities
# @require metaSEM
#
<- function(obj, adj = NULL) {
.sim_tsmasem1
# Sample correlations
<- list()
r_mats for (i in seq_len(nrow(obj))) {
<- metaSEM::vec2symMat(
r_mats[[i]] grepl("^r_", colnames(obj))],
obj[i, diag = FALSE
)
}<- metaSEM::list2matrix(r_mats, diag = FALSE)
r_vecs
# Sampling (co)variances of correlations
<- metaSEM::asyCov(
V_vecs x = r_mats,
n = obj[, "n"],
cor.analysis = TRUE,
as.matrix = TRUE,
acov = "weighted"
)
# Disattenuate sample correlations
<- obj[, grepl("^rww$", colnames(obj))]
rww if (length(rww) == 0) rww <- NULL
<- .disattenuate(r_vecs,
r_vecs rxx = obj[, "rxx"], rmm = obj[, "rmm"],
ryy = obj[, "ryy"], rww = rww,
adj = adj, vcov = FALSE)
# Disattentuate sampling (co)variances
<- .disattenuate(V_vecs,
V_vecs rxx = obj[, "rxx"], rmm = obj[, "rmm"],
ryy = obj[, "ryy"], rww = rww,
adj = adj, vcov = TRUE)
# Pool correlations
<-
fit1 try(
suppressWarnings(
::meta(
metaSEMy = r_vecs,
v = V_vecs,
RE.constraints = Diag(rep(paste0(0.1, "*Tau2"), ncol(r_vecs))),
RE.lbound = 1e-10,
silent = TRUE
)
), silent = TRUE
)
# Rerun in case of convergence failures
if (!inherits(fit1, "try-error") & !is.null(fit1$mx.fit$output$status$code)) {
<- 0
i while (!(fit1$mx.fit$output$status$code %in% c(0, 1)) & (i < 5)) {
<-
fit1 try(
suppressMessages(
::rerun(fit1, autofixtau2 = TRUE, extraTries = 10,
metaSEMsilent = TRUE)
)
)<- i + 1
i
}
}
# Return pooled results
if (inherits(fit1, "try-error") |
is.null(fit1$mx.fit$output$status$code)) {
<- list(status = 9)
out else {
} <- coef(fit1, select = "fixed")
pooledS <- metaSEM::vec2symMat(pooledS, diag = FALSE)
pooledS <- suppressWarnings(vcov(fit1, select = "fixed"))
aCov <- list(
out pooledS = pooledS,
aCov = aCov,
n = sum(obj[, "n"]),
k = nrow(obj),
rxx = weighted.mean(obj[, "rxx"], obj[, "n"]),
rmm = weighted.mean(obj[, "rmm"], obj[, "n"]),
ryy = weighted.mean(obj[, "ryy"], obj[, "n"]),
iter_sample = sum(obj[, "iter"]),
status = fit1$mx.fit$output$status$code
)if (!is.null(rww))
$rww <- weighted.mean(obj[, "rww"], obj[, "n"])
out
}return(out)
}
#
# Second stage in two-stage MASEM
#
# The pooled correlations and their asymptotic (co)variances
# are adjusted for unreliability before fitting the SEM.
#
# @param obj The results of the first step of MASEM
# as returned by .sim_tsmasem1()
# @param model The RAM specification for the SEM model
# as returned by metaSEM::lavaan2RAM()
# @param adj A list with logical indicators for each variable
# showing whether to adjust it for unreliability (optional)
# @return A data.frame with the parameter estimates
# @require metaSEM, OpenMx
#
<- function(obj, model, adj = NULL) {
.sim_tsmasem2
# Disattenuate pooled correlations
<- .disattenuate(OpenMx::vechs(obj$pooledS),
pooledS rxx = obj$rxx, rmm = obj$rmm,
ryy = obj$ryy, rww = obj$rww,
adj = adj, vcov = FALSE)
<- metaSEM::vec2symMat(pooledS, diag = FALSE)
pooledS
# Disattenuate asymptotic (co)variances
<- .disattenuate(OpenMx::vech(obj$aCov),
aCov rxx = obj$rxx, rmm = obj$rmm,
ryy = obj$ryy, rww = obj$rww,
adj = adj, vcov = TRUE)
<- metaSEM::vec2symMat(aCov, diag = TRUE)
aCov
# OpenMx matrices for indirect effects
<- .get_INDmatrix(model)
ind
# Fit SEM
<-
fit2 try(
suppressWarnings(
::wls(
metaSEMCov = pooledS,
aCov = aCov,
n = obj$n,
RAM = model,
cor.analysis = TRUE,
intervals.type = "LB",
mx.algebras = ind
)
), silent = TRUE
)
# Rerun in case of convergence failures
if (!inherits(fit2, "try-error") & !is.null(fit2$mx.fit$output$status$code)) {
<- 0
i while (!(fit2$mx.fit$output$status$code %in% c(0, 1)) & (i < 5)) {
<- suppressMessages(
fit2 ::rerun(fit2, extraTries = 10, intervals = TRUE, silent = TRUE)
metaSEM
) <- i + 1
i
}
}
# Return SEM results
if (inherits(fit2, "try-error") |
is.null(fit2$mx.fit$output$status$code) |
is.null(fit2$mx.fit$output$confidenceIntervals)) {
<- c(status = 9)
out else {
} <- fit2$mx.fit$output$confidenceIntervals
pars if (length(ind) == 3)
rownames(pars)[nrow(pars) - c(2, 1)] <- c("IND1", "IND2")
rownames(pars)[nrow(pars)] <- "IND"
<- c(status = fit2$mx.fit$output$status$code,
out iter_sample = obj$iter_sample)
for (i in rownames(pars)) {
paste0(i, "_est")] <- pars[i, "estimate"]
out[paste0(i, "_lb")] <- pars[i, "lbound"]
out[paste0(i, "_ub")] <- pars[i, "ubound"]
out[
}
}return(out)
}
#
# One-stage MASEM with adjustments for unreliability
#
#
# @param obj A matrix with sample statistics for multiple samples
# as returned by .sim_samples()
# @param model The RAM specification for the SEM model
# as returned by metaSEM::lavaan2RAM()
# @param adj A list with logical indicators for each variable
# showing whether to adjust it for unreliability (optional)
# @param individual A logical indicating whether to adjust
# the correlations individually (TRUE) or adjust the
# model implied correlation matrix (FALSE)
# @return A data.frame with the parameter estimates
# @require metaSEM, OpenMx
#
<- function(obj, model, adj = NULL, individual = TRUE) {
.sim_osmasem
# Sample correlations
<- list()
r_mats for (i in seq_len(nrow(obj))) {
<- metaSEM::vec2symMat(
r_mats[[i]] grepl("^r_", colnames(obj))],
obj[i, diag = FALSE
)
}<- metaSEM::list2matrix(r_mats, diag = FALSE)
r_vecs
# Sampling (co)variances of correlations
<- metaSEM::asyCov(
V_vecs x = r_mats,
n = obj[, "n"],
cor.analysis = TRUE,
as.matrix = TRUE,
acov = "weighted"
)
# Adjust correlations individually
<- obj[, grepl("^rww$", colnames(obj))]
rww if (length(rww) == 0) rww <- NULL
if (individual) {
# Disattenuate sample correlations
<- .disattenuate(r_vecs,
r_vecs rxx = obj[, "rxx"], rmm = obj[, "rmm"],
ryy = obj[, "ryy"], rww = rww,
adj = adj, vcov = FALSE)
# Disattentuate sampling (co)variances
<- .disattenuate(V_vecs,
V_vecs rxx = obj[, "rxx"], rmm = obj[, "rmm"],
ryy = obj[, "ryy"], rww = rww,
adj = adj, vcov = TRUE)
# Disattenuation factors for adjustments of implied correlations
else {
}
<- .disattenuate(r_vecs,
att_vecs rxx = weighted.mean(obj[, "rxx"], obj[, "n"]),
rmm = weighted.mean(obj[, "rmm"], obj[, "n"]),
ryy = weighted.mean(obj[, "ryy"], obj[, "n"]),
rww = weighted.mean(rww, obj[, "n"]),
adj = adj, vcov = FALSE, return_att = TRUE)
<- metaSEM::vec2symMat(att_vecs, diag = FALSE)
att_mat
}
# OpenMx matrices for indirect effects
<- .get_INDmatrix(model)
ind $CI <- OpenMx::mxCI(names(ind))
ind
# Generate OpenMx model
<- metaSEM::osmasem(
mx.model RAM = model,
RE.type = "Diag",
data = list(
data = cbind(r_vecs, V_vecs),
n = obj[, "n"],
obslabels = NULL,
ylabels = colnames(r_vecs),
vlabels = colnames(V_vecs)
),intervals.type = "LB",
mxModel.Args = ind,
run = FALSE
)
# Constrain random effects
$vecTau1$labels <- cbind(rep("Tau1", length(mx.model$vecTau1$labels)))
mx.model
# Adjust model implied correlations
if (!individual) {
$G <- metaSEM::as.mxMatrix(att_mat, name = "G")
mx.model$impliedR <- OpenMx::mxAlgebra(
mx.model%*% solve(Iden - Amatrix)) %&% Smatrix) * G,
((Fmatrix name = "impliedR"
)
}
# Fit SEM
<- list(mx.fit =
fit try(
suppressWarnings(
::mxRun(mx.model, intervals = TRUE, silent = TRUE,
OpenMxsuppressWarnings = TRUE)
), silent = TRUE
)
)class(fit) <- "osmasem"
# Rerun in case of convergence failures
if (!inherits(fit$mx.fit, "try-error") & !is.null(fit$mx.fit$output$status$code)) {
<- 0
i while (!(fit$mx.fit$output$status$code %in% c(0, 1)) & (i < 5)) {
<- suppressMessages(
fit ::rerun(fit, autofixtau2 = TRUE, extraTries = 10,
metaSEMintervals = TRUE, silent = TRUE)
) <- i + 1
i
}
}
# Return SEM results
if (inherits(fit$mx.fit, "try-error") |
is.null(fit$mx.fit$output$status$code) |
is.null(fit$mx.fit$output$confidenceIntervals)) {
<- c(status = 9)
out else {
} <- fit$mx.fit$output$confidenceIntervals
pars <- pars[grepl("^osmasem.(Amatrix|IND)", rownames(pars)), ]
pars <- c("X", "M", "W", "Y")
lbl if (ncol(r_vecs) < 4) lbl <- lbl[-3]
rownames(pars) <- c(paste0(rep(lbl, length(lbl)), "ON",
rep(lbl, each = length(lbl))),
names(ind)[-length(ind)])
<- pars[apply(pars, 1, \(x) { !all(x %in% c(0, NA)) }), ]
pars <- c(status = fit$mx.fit$output$status$code,
out iter_sample = sum(obj[, "iter"]))
for (i in rownames(pars)) {
paste0(i, "_est")] <- pars[i, "estimate"]
out[paste0(i, "_lb")] <- pars[i, "lbound"]
out[paste0(i, "_ub")] <- pars[i, "ubound"]
out[
}
}return(out)
}
#########################################################
# RAM specification for SEM model
<- ifelse(!is.null(WONM),
type "sequential",
ifelse(!is.null(WONX), "parallel", "simple"))
<- .get_model(type)
model
# Number of errors
<- list(masem1a = 0, masem2a = 0,
status masem1b = 0, masem2b = 0,
masem1c = 0, masem2c = 0,
masem0d = 0, masem0e = 0,
masem0f = 0)
# Reestimate models until all have converged
repeat {
# Calculate correlation matrices from randomly
# generated samples in the meta-analysis
<- .sim_samples(MONX = MONX, YONM = YONM, YONX = YONX,
simdat WONX = WONX, YONW = YONW, WONM = WONM,
k = k, n = n, tau = tau, rxx = rxx,
rmm = rmm, rww = rww, ryy = ryy,
tau_max = tau_max)
<- c()
out <- list()
stage1
# The meta-analysis is estimated six times:
# (a) TSMASEM without adjustments
# (b) TSMASEM with adjustments for all variables
# at the first stage
# (c) TSMASEM with adjustments for all variables
# at the second stage
# (d) OSMASEM without adjustments
# (e) OSMASEM with individual adjustments for all variables
# (f) OSMASEM with adjustments of implied correlations
# for all variables
# If a model results in an error (e.g., failed to converge),
# the entire estimation restarts with new random draws.
# Estimation error occurred
<- lapply(status, \(x) 0)
e
# (a) TSMASEM without adjustments
$masem_a <- .sim_tsmasem1(simdat,
stage1adj = list(rxx = FALSE, rmm = FALSE,
rww = FALSE, ryy = FALSE))
$masem1a <- stage1$masem_a$status
eif (stage1$masem_a$status <= 1) {
$masem_a <- .sim_tsmasem2(stage1$masem_a, model,
outadj = list(rxx = FALSE, rmm = FALSE,
rww = FALSE, ryy = FALSE))
$masem2a <- out$masem_a["status"]
e
}
# (b) TSMASEM with adjustments for all variables at the first stage
$masem_b <- .sim_tsmasem1(simdat,
stage1adj = list(rxx = TRUE, rmm = TRUE,
rww = TRUE, ryy = TRUE))
$masem1b <- stage1$masem_b$status
eif (stage1$masem_b$status <= 1) {
$masem_b <- .sim_tsmasem2(stage1$masem_b, model,
outadj = list(rxx = FALSE, rmm = FALSE,
rww = FALSE, ryy = FALSE))
$masem2b <- out$masem_b["status"]
e
}
# (c) TSMASEM with adjustments for all variables at the second stage
$masem1c <- e$masem1a
eif (stage1$masem_a$status <= 1) {
$masem_c <- .sim_tsmasem2(stage1$masem_a, model,
outadj = list(rxx = TRUE, rmm = TRUE,
rww = TRUE, ryy = TRUE))
$masem2c <- out$masem_c["status"]
e
}
# (d) OSMASEM without adjustments
$masem_d <- .sim_osmasem(simdat, model,
outadj = list(rxx = FALSE, rmm = FALSE,
rww = FALSE, ryy = FALSE))
$masem0d <- out$masem_d["status"]
e
# (i) OSMASEM with individual adjustments for all variables
$masem_e <- .sim_osmasem(simdat, model,
outadj = list(rxx = TRUE, rmm = TRUE,
rww = TRUE, ryy = TRUE),
individual = TRUE)
$masem0e <- out$masem_e["status"]
e
# (j) OSMASEM with adjustments of implied correlations
# for all variables
$masem_f <- .sim_osmasem(simdat, model,
outadj = list(rxx = TRUE, rmm = TRUE,
rww = TRUE, ryy = TRUE),
individual = FALSE)
$masem0f <- out$masem_f["status"]
e
# Record errors
for (i in names(e)) {
if (e[[i]] %in% 2:9) status[[i]] <- e[[i]]
}
# Repeat meta-analysis unless converged
if (all(sapply(e, \(x) (x %in% 0:1))))
break
}
# True indirect effect
<- switch(type,
IND simple = c(IND = MONX * YONM),
parallel = c(IND1 = MONX * YONM,
IND2 = WONX * YONW,
IND = MONX * YONM + WONX * YONW),
sequential = c(IND = MONX * WONM * YONW))
# Return results
<- unlist(status)
status names(status) <- paste0(names(status), ".status")
<- lapply(out, \(o) { o[-c(1, 2)] })
out <- c(id = rid,
out MONX = MONX, YONM = YONM, YONX = YONX,
WONX = WONX, YONW = YONW, WONM = WONM,
IND,1, "tau"],
simdat[k = stage1$masem_a$k,
n = stage1$masem_a$n,
rxx = stage1$masem_a$rxx,
rmm = stage1$masem_a$rmm,
rww = stage1$masem_a$rww,
ryy = stage1$masem_a$ryy,
unlist(out),
iter_samp = stage1$masem_a$iter_samp,
|>
status) unlist()
return(out)
}
#
# Run simulation for a given condition
#
# @param R The number of replications (optional)
# @param cores The number of cores to use for parallel processing
# (optional)
# @param seed A random seed (optional)
# @param ... Arguments passed to sim_wrapper()
# @require parallel, pbapply
#
<- function(R = 1000, cores = parallel::detectCores(),
sim_run seed = 9789, ...) {
# Don't run the simulation
if (isTRUE(TEST)) return(NULL)
# Create clusters
<- parallel::makeCluster(cores)
cl
# Set seed
::clusterSetRNGStream(cl = cl, iseed = seed)
parallel
# Run simulation
::clusterExport(cl, "sim_wrapper")
parallelinvisible(parallel::clusterEvalQ(
cl, library("OpenMx"); library("metaSEM");
{ library("MASS"); library("Matrix") }
))<-
data ::pbsapply(
pbapplyX = 1:R,
FUN = \(i) { sim_wrapper(rid = i, ...) },
cl = cl
)<- as.data.frame(t(data))
data
# Stop clusters
::stopCluster(cl)
parallel
return(data)
}
<- TRUE TEST
<- list()
simdat
# Low reliability
$simple1 <- sim_run(MONX = .30, YONM = .30, YONX = .00,
simdatrmm = c(.65, .75),
R = 1000, seed = 4312)
# High reliability
$simple2 <- sim_run(MONX = .30, YONM = .30, YONX = .00,
simdatrmm = c(.85, .95),
R = 1000, seed = 4312)
# Small indirect effect
$simple3 <- sim_run(MONX = .30, YONM = .075, YONX = .00,
simdatR = 1000, seed = 8214)
# Large indirect effect
$simple4 <- sim_run(MONX = .30, YONM = .675, YONX = .00,
simdatR = 1000, seed = 8214)
# Small a-path
$simple5 <- sim_run(MONX = .15, YONM = .60, YONX = .00,
simdatR = 1000, seed = 3419)
# Large a-path
$simple6 <- sim_run(MONX = .60, YONM = .15, YONX = .00,
simdatR = 1000, seed = 3419)
# Low reliability
$parallel1 <- sim_run(MONX = .30, YONM = .30, YONX = .00,
simdatWONX = .30, YONW = .30,
rmm = c(.65, .75),
rww = c(.65, .75),
R = 1000, seed = 4312)
# High reliability
$parallel2 <- sim_run(MONX = .30, YONM = .30, YONX = .00,
simdatWONX = .30, YONW = .30,
rmm = c(.85, .95),
rww = c(.85, .95),
R = 1000, seed = 4312)
# Small indirect effect
$parallel3 <- sim_run(MONX = .30, YONM = .075, YONX = .00,
simdatWONX = .30, YONW = .075,
R = 1000, seed = 8214)
# Large indirect effect
$parallel4 <- sim_run(MONX = .30, YONM = .675, YONX = .00,
simdatWONX = .30, YONW = .675,
R = 1000, seed = 8214)
# Small a-path
$parallel5 <- sim_run(MONX = .15, YONM = .60, YONX = .00,
simdatWONX = .15, YONW = .60,
R = 1000, seed = 3419)
# Large a-path
$parallel6 <- sim_run(MONX = .60, YONM = .15, YONX = .00,
simdatWONX = .60, YONW = .15,
R = 1000, seed = 3419)
# Low reliability
$sequential1 <- sim_run(MONX = .30, YONX = .00,
simdatWONM = .44, YONW = .30,
rmm = c(.65, .75),
rww = c(.65, .75),
R = 1000, seed = 4312)
# High reliability
$sequential2 <- sim_run(MONX = .30, YONX = .00,
simdatWONM = .44, YONW = .30,
rmm = c(.85, .95),
rww = c(.85, .95),
R = 1000, seed = 4312)
# Small indirect effect
$sequential3 <- sim_run(MONX = .30, YONX = .00,
simdatWONM = .44, YONW = .075,
R = 1000, seed = 8214)
# Large indirect effect
$sequential4 <- sim_run(MONX = .30, YONX = .00,
simdatWONM = .44, YONW = .675,
R = 1000, seed = 8214)
# Small a-path
$sequential5 <- sim_run(MONX = .15, YONX = .00,
simdatWONM = .44, YONW = .60,
R = 1000, seed = 3419)
# Large a-path
$sequential6 <- sim_run(MONX = .60, YONX = .00,
simdatWONM = .44, YONW = .15,
R = 1000, seed = 3419)
if (isFALSE(TEST))
saveRDS(simdat, file = here::here("data", "simdat.rds"))