Simulate High-Dimensional Mediation Data

Generates a data set from the mediation models studied in the article, for a continuous, binary, or count outcome, under the homogeneous or contrasting (heterogeneous) mediation patterns. The mediators carry an AR(1) correlation structure, a handful are truly active, and the rest are null, so the data exercise exactly the situation the power-enhanced tests are designed for. This is the generator behind the package's example data sets and behind pe_power_curve().

Usage

simulate_mediation_data(
  n,
  p,
  outcome = c("continuous", "binary", "count"),
  pattern = c("homogeneous", "contrasting", "heterogeneous"),
  c1 = 1,
  c2 = 0.5,
  q = 1L,
  d = 0L,
  alpha_m = NULL,
  tau = NULL,
  rho = 0.5,
  sigma_y = 0.5,
  alpha_z = NULL,
  Gamma_z = NULL,
  seed = NULL
)

Arguments

n: Number of observations.
p: Number of candidate mediators.
outcome: Outcome type: "continuous", "binary", or "count".
pattern: Mediation pattern: "homogeneous" (active effects share a sign) or "contrasting" (active effects of opposite sign). "heterogeneous" is accepted as a synonym for "contrasting". Ignored when alpha_m is supplied directly. The presets are the article's. With the default tau, the contrasting presets for the continuous and binary outcomes cancel to a total indirect effect of zero, while the count preset, (0, 0, 0, 0.8, -0.7), gives a total indirect effect of -0.03 times c1.
c1: Scale on the exposure-on-mediator coefficients \(\Gamma_x = c_1 \tau\). c1 = 0 makes every mediator inactive (the global null). Default 1.
c2: The direct effect \(\alpha_x\) of each exposure on the outcome. Default 0.5.
q: Number of exposures. Default 1. For q > 1 you must supply tau as a q by p matrix.
d: Number of confounders. Default 0 (none). When positive, confounders are drawn standard normal and enter both the mediator and the outcome model through Gamma_z and alpha_z, whose coefficients are themselves drawn from U(-0.5, 0.5) unless supplied. The article's "with confounders" scenario is d = 5.
alpha_m: Optional numeric vector of p finite mediator-on-outcome coefficients, overriding the pattern preset.
tau: Optional base pattern for \(\Gamma_x\), finite numbers. By default it is the article's (0.1, 0.2, 0.3, 0.4, 0.5, noise), for q = 1, with the noise loadings drawn from N(0, 0.1^2) for a continuous outcome and N(0, 0.01^2) for a binary or count one, the scales of the article's respective Monte Carlo designs. Required when q > 1, as a q by p matrix (one row per exposure) or a vector of its values in column order; when q = 1 a vector of length p serves.
rho: AR(1) correlation parameter for the mediator noise. Default 0.5.
sigma_y: Outcome noise standard deviation for a continuous outcome. Default 0.5.
alpha_z, Gamma_z: Optional confounder coefficients: a vector of d finite numbers and a d by p matrix (when d = 1, a vector of length p serves for Gamma_z). When d > 0 and either is left NULL, it is drawn from U(-0.5, 0.5), the article's Monte Carlo design, afresh on every call. Gamma_z is then drawn as d values recycled across the p mediator columns, reproducing the article's code: each confounder loads equally on every mediator, so the confounders enter the mediators as a single common factor. Supply Gamma_z yourself for d by p independent loadings. Both are used only when d > 0; supplying either with d = 0, where it would be ignored, is an error.
seed: Optional integer seed. When supplied it is set for the duration of the call and the caller's random number generator state is restored on exit. NULL (the default) sets no seed: the draws come from the session's random number stream, which the call advances as any random function does.

Value

A list with the simulated data and the truth used to generate it: X (n by q), M (n by p), Y (length n), Z (n by d, or NULL), the coefficient values alpha_m, Gamma_x, alpha_x, the confounder coefficients alpha_z and Gamma_z that were used (NULL when d = 0, so a drawn pair is recoverable), the total indirect effect beta (Gamma_x times alpha_m), and active_mediators (the indices of the truly active mediators, those nonzero in both paths), together with n, p, q, outcome, and pattern.

For a count outcome beta is the value before the overflow guard described under Details. When that guard fires, it rescales the log-mean by a factor the sample itself determines, so the coefficients realized in that replication, and with them the total indirect effect, are beta times that factor rather than beta.

Details

One exposure column (or \(q\) of them) is drawn standard normal. The mediators follow \(M = X \Gamma_x + Z \Gamma_z + \varepsilon_m\) with \(\varepsilon_m \sim N(0, \Sigma)\), \(\Sigma\) the AR(1) covariance mediation_ar1_cov() with parameter rho. The exposure-on-mediator coefficient is \(\Gamma_x = c_1 \tau\), where the base pattern \(\tau\) has small increasing loadings on the first five mediators and random noise loadings on the rest, so c1 scales the overall exposure-to-mediator signal (and c1 = 0 gives the global null of no mediation). The outcome is generated from its model with mediator coefficients \(\alpha_m\) (the preset for the chosen pattern, or a user-supplied vector) and direct effect \(\alpha_x = c_2\):

continuous: \(Y = M\alpha_m + X\alpha_x + Z\alpha_z + \varepsilon_y\), \(\varepsilon_y \sim N(0, \sigma_y^2)\);
binary: \(Y \sim \mathrm{Bernoulli}(\mathrm{logit}^{-1}( M\alpha_m + X\alpha_x + Z\alpha_z))\);
count: \(Y \sim \mathrm{Poisson}(\exp(\eta))\), where the log-mean \(\eta\) passes through the article's overflow guard: if its largest value exceeds 5, the whole vector is mapped affinely onto the interval from -5 to 5 before exponentiating, and otherwise it is left alone. The guard tests only the maximum, so a run of very negative values passes through untouched.

With d > 0 the confounders are standard normal and enter the mediators as \(Z \Gamma_z\) and the outcome as \(Z \alpha_z\), with both coefficient sets drawn from U(-0.5, 0.5) unless given. They leave the total indirect effect alone, since that runs through \(\Gamma_x\) and \(\alpha_m\), but they correlate the mediators with each other and with the outcome, which is what makes adjusting for them necessary.

Every design argument is checked before any random number is drawn: a coefficient vector or matrix of the wrong length or shape, or one holding a missing or infinite value, stops with a message naming it.

References

Yu, X., & Kelley, K. (in press). Power Enhancement in High-Dimensional Heterogeneous Mediation Analysis. Journal of the American Statistical Association.

Author

Xiufan Yu and Ken Kelley

Examples

set.seed(113)
d <- simulate_mediation_data(n = 100, p = 50, outcome = "continuous",
                             pattern = "contrasting")
dim(d$M)
#> [1] 100  50
d$active_mediators       # truly active mediators
#> [1] 1 2 3 4
round(d$beta, 10)        # total indirect effect, zero up to rounding error
#> [1] 0