simulate_estimate: Simulate data and estimate a model repeatedly

In covsim: VITA, IG, PLSIM and VM Simulation for Given Covariance and Marginals

Simulate data and estimate a model repeatedly

Description

Runs a Monte Carlo study: for every calibrated distribution in dist and every sample size in n, draws reps samples with simulate, fits model to each sample, and collects what extractor returns into one data frame.

Usage

simulate_estimate(
  dist,
  model,
  n,
  reps,
  extractor,
  ...,
  fit = "lavaan",
  cores = 1,
  seed = NULL
)

Arguments

dist	A calibration object (from calibrate_ig, calibrate_plsim, calibrate_vm or calibrate_vita), or a named list of such objects to compare several distributions.
model	Model syntax (lavaan syntax for fit = "lavaan" and fit = "magmaan"), or whatever a user-supplied fit function expects. Variable names must match the column names of sigma.target used in calibration.
n	Vector of sample sizes.
reps	Number of replications per distribution and sample size.
extractor	A function of one fitted object that returns what is to be kept from that replication, as a named vector or list of scalars. Each element becomes a column. The researcher decides what to extract; see the examples for parameter estimates, standard errors, p-values, fit measures and test statistics with lavaan and magmaan. Use function(fit) list(fit = fit) to keep the fitted objects themselves (as a list column).
...	Further arguments passed to lavaan::sem() or magmaan::magmaan(), for example estimator = "MLR", se = "robust.sem", start = "simple" or control = list(iter.max = 500). Not used with a user-supplied fit function, which can set such options itself.
fit	"lavaan", "magmaan", or a function function(data, model) returning a fitted object.
cores	Number of cores. Values above 1 use mclapply (not available on Windows, where the replications run sequentially).
seed	Optional integer seed. If NULL, a seed is drawn from the current random number generator.

Details

Estimation is done with lavaan (fit = "lavaan", the default), with magmaan (fit = "magmaan", if installed), or with any function function(data, model). With magmaan the model is constructed once per distribution with magmaan::magmaan_model() and reused for every sample, which avoids repeated parsing and setup.

Each replication draws its data from its own L'Ecuyer-CMRG random number stream, determined by seed and the position of the replication in the design. Results are therefore identical for any value of cores. Errors and warnings from fitting or extraction are caught and recorded per replication rather than stopping the run.

Value

A data frame with one row per replication and columns dist, n, rep, converged, error, warning, followed by the extracted values. The seed and package versions are stored in attr(, "info").

Extractors for magmaan

A magmaan fit f gives the parameter table with summary(f)$coefficients (columns lhs, op, rhs, free, est, se, pvalue, ...), its score and likelihood ratio tests with summary(f)$tests, and fit measures with magmaan::fit_measures(f). Note that by default fit_measures() reports robust RMSEA and CFI and a bias-corrected SRMR under the plain names; magmaan::fit_measures(f, lavaan_compat = "MLR") gives lavaan's definitions and names (rmsea, rmsea.robust, cfi.robust, ...). magmaan's default standard errors are robust (sandwich), matching lavaan with estimator = "MLR" rather than plain ML. magmaan fits a mean structure, so its parameter table also contains the intercepts.

See Also

calibrate_ig, calibrate_plsim, calibrate_vm, calibrate_vita

Examples

model <- "F =~ x1 + x2 + x3 + x4"
sigma <- matrix(0.5, 4, 4); diag(sigma) <- 1
dimnames(sigma) <- list(paste0("x", 1:4), paste0("x", 1:4))
dists <- list(
  IG = calibrate_ig(sigma, skewness = rep(2, 4), excesskurtosis = rep(7, 4)),
  VM = calibrate_vm(sigma, skewness = rep(2, 4), excesskurtosis = rep(7, 4)))

# rejection rate of the ML chi-square test
chisq <- function(fit) lavaan::fitMeasures(fit, c("chisq", "df", "pvalue"))
res <- simulate_estimate(dists, model, n = 200, reps = 20,
                         extractor = chisq, seed = 1)
aggregate(pvalue ~ dist, data = res, FUN = function(p) mean(p < 0.05))

# RMSEA and CFI, ordinary and robust (MLR)
fitind <- function(fit)
  lavaan::fitMeasures(fit, c("rmsea", "cfi", "rmsea.robust", "cfi.robust"))
res <- simulate_estimate(dists, model, n = 200, reps = 20, estimator = "MLR",
                         extractor = fitind, seed = 1)
aggregate(cbind(rmsea, rmsea.robust) ~ dist, data = res, FUN = mean)

# estimate, standard error and p-value of one loading
loading <- function(fit) {
  pe <- lavaan::parameterEstimates(fit)
  pe <- pe[pe$lhs == "F" & pe$op == "=~" & pe$rhs == "x2", ]
  c(est = pe$est, se = pe$se, pvalue = pe$pvalue)
}
res <- simulate_estimate(dists$VM, model, n = 200, reps = 20,
                         extractor = loading, seed = 1)
colMeans(res[, c("est", "se")])   # mean estimate vs. mean standard error
sd(res$est)                       # empirical standard error

# the same with magmaan
if (requireNamespace("magmaan", quietly = TRUE) &&
    utils::packageVersion("magmaan") >= "0.2.0") {
  mag_loading <- function(fit) {
    cf <- summary(fit)$coefficients
    cf <- cf[cf$lhs == "F" & cf$op == "=~" & cf$rhs == "x2", ]
    c(est = cf$est, se = cf$se, pvalue = cf$pvalue)
  }
  mag_tests <- function(fit) {
    tests <- summary(fit)$tests
    stats::setNames(tests$pvalue, paste(tests$test, tests$reference, sep = "_"))
  }
  res <- simulate_estimate(dists$VM, model, n = 200, reps = 20, fit = "magmaan",
                           extractor = function(fit) c(mag_loading(fit), mag_tests(fit)),
                           seed = 1)
  colMeans(res[, c("est", "se", "score_peba4")])
}

covsim documentation built on Oct. 29, 2026, 5:06 p.m.