Monte Carlo studies with covsim

Njål Foldnes

2026-10-29

This vignette walks through a complete simulation study in covsim: specify a population, calibrate several non-normal data generators with the same covariance matrix, skewness and kurtosis, and estimate a structural equation model on many simulated samples.

The study is done in two ways. First by hand: the calibration objects produce lists of samples, as rIG(), rPLSIM() and vita() did in earlier versions of covsim, and the researcher writes the loop that fits the model and collects the results. Then with simulate_estimate(), which runs that loop and takes care of the bookkeeping. Estimation is done with lavaan and with magmaan, which is much faster.

The question studied is a classical one: how often do chi-square tests of model fit reject a correctly specified model when the data are non-normal? The answer depends on the data generator, even when all generators share the same first four marginal moments. This is why it is useful to compare several generators in the same study.

The number of replications is kept small so the vignette builds quickly. A real study would use 1000 or more replications per condition.

The population

A two-factor model with three indicators per factor, standardized loadings 0.7 and factor correlation 0.4. The variable names of sigma must match the model syntax.

library(covsim)

model <- "F1 =~ x1 + x2 + x3
          F2 =~ x4 + x5 + x6"

L <- matrix(0, 6, 2)
L[1:3, 1] <- L[4:6, 2] <- 0.7
sigma <- L %*% matrix(c(1, 0.4, 0.4, 1), 2) %*% t(L)
diag(sigma) <- 1
dimnames(sigma) <- list(paste0("x", 1:6), paste0("x", 1:6))

skew <- rep(2, 6)
kurt <- rep(7, 6)   # excess kurtosis

Calibrating the data generators

Each calibrate_*() function returns a calibration object. Calibration is done once, and the object can then produce any number of samples.

dists <- list(
  IG    = calibrate_ig(sigma, skew, kurt),
  PLSIM = calibrate_plsim(sigma, skew, kurt, verbose = FALSE),
  VM    = calibrate_vm(sigma, skew, kurt)
)

The three generators have the same covariance matrix, skewness and excess kurtosis, but differ in their multivariate dependence. A large sample confirms the marginal moments:

moments <- function(x) c(skew = mean(scale(x)^3), exkurt = mean(scale(x)^4) - 3)
set.seed(1)
round(sapply(dists, function(d) moments(simulate(d, N = 1e5)[[1]][, 1])), 2)
#>          IG PLSIM   VM
#> skew   2.02  2.01 2.03
#> exkurt 7.47  7.14 7.22

VITA (vine copulas) can be added in the same way. Its calibration takes longer, so it is not run here:

margins <- rep(list(list(distr = "gamma", shape = 1)), 6)  # skewness 2
dists$VITA <- calibrate_vita(margins, sigma, Nmax = 1e5, verbose = FALSE)

Note that VITA works with the margins as distributions, so their skewness and kurtosis are given by the chosen family (here a gamma margin with skewness 2 and excess kurtosis 6), and the covariance matrix must be on the scale of those margins.

Part 1: the study by hand

A list of samples

simulate() draws nsim samples of size N from a calibration object and returns them as a list of data matrices with the variable names of sigma. The same seed gives the same samples:

samples <- simulate(dists$PLSIM, nsim = 50, N = 200, seed = 1)
length(samples)
#> [1] 50
head(samples[[1]], 3)
#>              x1          x2         x3        x4         x5        x6
#> [1,] -0.4101574  0.41281098 -0.2509398 1.0874435 -0.3141345 -0.393562
#> [2,]  0.0121689 -0.14186017 -1.2891077 3.1039231 -0.2188600 -0.367604
#> [3,] -0.4030072  0.09774435 -0.3928586 0.3584605  1.1263063  2.835460

The functions of earlier covsim versions still work, with the same arguments and output. They calibrate and draw in one call, using the same fast calibration as calibrate_plsim():

set.seed(1)
old_way <- rPLSIM(200, sigma, skew, kurt, reps = 50, verbose = FALSE)
length(old_way$samples)
#> [1] 50

Each call to rPLSIM() calibrates once and then draws all reps samples. When the same population is needed several times, for example for several sample sizes, it is cheaper to calibrate once with calibrate_plsim() and call simulate() for each sample size, since every new call to rPLSIM() repeats the calibration.

Fitting the model and collecting the results

The researcher now has full control: any estimator can be used and anything can be stored. Here each sample is fitted with lavaan using MLR, and the p-values of the ML test and of the scaled (Yuan-Bentler) test are kept. A fit that fails gives NA instead of stopping the loop:

fit_tests <- function(x) {
  fit <- tryCatch(lavaan::sem(model, data = as.data.frame(x), estimator = "MLR"),
                  error = function(e) NULL)
  if (is.null(fit))
    return(c(pvalue = NA, pvalue.scaled = NA))
  lavaan::fitMeasures(fit, c("pvalue", "pvalue.scaled"))
}

by_hand <- lapply(dists, function(cal) {
  samples <- simulate(cal, nsim = 50, N = 200, seed = 1)
  t(sapply(samples, fit_tests))
})
head(by_hand$VM, 3)
#>         pvalue pvalue.scaled
#> [1,] 0.7330111     0.6121612
#> [2,] 0.5085567     0.5690397
#> [3,] 0.2463290     0.2917469

Rejection rates at the 5% level:

round(sapply(by_hand, function(p) colMeans(p < 0.05, na.rm = TRUE)), 2)
#>                 IG PLSIM   VM
#> pvalue        0.16  0.24 0.20
#> pvalue.scaled 0.12  0.14 0.16

This is all that is needed for a small study. A larger study with several sample sizes, many conditions and parallel computing also needs a design table, a random number stream for each condition (so that the results do not depend on the order or the number of cores), a record of errors and warnings in each replication, and a way to combine everything into one data set. simulate_estimate() does this.

Part 2: the study with simulate_estimate()

The extractor: what to keep from each fit

simulate_estimate() draws reps samples of each size in n from each distribution, fits the model to each sample, and passes the fitted object to extractor. The extractor is a function written by the researcher, in the same way as fit_tests() above: it receives one fitted object and returns a named vector (or list) of the quantities to keep; each element becomes a column of the result. Which estimates, standard errors, test statistics or fit indices to keep is thus decided by the researcher, not by covsim.

The extractor must be written for the estimation package used. The examples below use lavaan’s accessor functions fitMeasures() and parameterEstimates(); magmaan versions follow later.

Estimation with lavaan

Arguments in ... go to lavaan::sem(). Here the model is estimated with MLR. The first extractor keeps the statistic and p-value of the normal theory ML chi-square test and of the Yuan-Bentler (MLR) scaled test:

chisq_tests <- function(fit) {
  lavaan::fitMeasures(fit, c("chisq", "pvalue", "chisq.scaled", "pvalue.scaled"))
}

res <- simulate_estimate(dists, model, n = c(200, 500), reps = 50,
                         extractor = chisq_tests, estimator = "MLR", seed = 1)
head(res)
#>   dist   n rep converged error warning     chisq    pvalue chisq.scaled
#> 1   IG 200   1      TRUE  <NA>    <NA> 13.139242 0.1071477    17.303819
#> 2   IG 200   2      TRUE  <NA>    <NA> 12.344890 0.1364663    14.964371
#> 3   IG 200   3      TRUE  <NA>    <NA>  4.089174 0.8489895     6.719689
#> 4   IG 200   4      TRUE  <NA>    <NA>  6.938831 0.5432465     7.075888
#> 5   IG 200   5      TRUE  <NA>    <NA>  6.845373 0.5534041     7.123690
#> 6   IG 200   6      TRUE  <NA>    <NA>  3.521407 0.8975219     4.980416
#>   pvalue.scaled
#> 1    0.02709642
#> 2    0.05984197
#> 3    0.56715588
#> 4    0.52846756
#> 5    0.52334827
#> 6    0.75966723

There is one row per replication. The columns converged, error and warning record what happened in each fit; a failed fit gives NA values and an error message instead of stopping the study.

table(res$converged, useNA = "ifany")
#> 
#> TRUE 
#>  300

Rejection rates at the 5% level, by distribution and sample size:

reject <- function(p) mean(p < 0.05, na.rm = TRUE)
aggregate(cbind(ML = pvalue, MLR = pvalue.scaled) ~ dist + n, data = res,
          FUN = reject)
#>    dist   n   ML  MLR
#> 1    IG 200 0.04 0.06
#> 2 PLSIM 200 0.16 0.14
#> 3    VM 200 0.16 0.16
#> 4    IG 500 0.06 0.04
#> 5 PLSIM 500 0.18 0.06
#> 6    VM 500 0.18 0.10

The rejection rates differ between the generators although they share the same covariance matrix, skewness and kurtosis. With 50 replications each rate has a standard error of about 3 percentage points, so a real study needs many more replications before such differences can be interpreted. The rates also differ from those in Part 1: simulate_estimate() gives each replication its own random number stream, so seed = 1 does not produce the same samples as simulate(cal, seed = 1).

Fit indices

To study RMSEA and CFI, ordinary and robust, extract those instead:

fit_indices <- function(fit) {
  lavaan::fitMeasures(fit, c("rmsea", "rmsea.robust", "cfi", "cfi.robust"))
}

res_fi <- simulate_estimate(dists, model, n = 200, reps = 50,
                            extractor = fit_indices, estimator = "MLR", seed = 1)
aggregate(cbind(rmsea, rmsea.robust, cfi, cfi.robust) ~ dist, data = res_fi,
          FUN = mean)
#>    dist      rmsea rmsea.robust       cfi cfi.robust
#> 1    IG 0.01698970   0.02116894 0.9955399  0.9941137
#> 2 PLSIM 0.03529958   0.02917490 0.9878785  0.9895980
#> 3    VM 0.03544332   0.03302894 0.9856679  0.9862149

Parameter estimates, standard errors and coverage

An extractor can pick out a single parameter. Here it keeps the estimate, standard error and 95% confidence limits of the loading of x2, whose population value is 1 (the indicators have unit variance and standardized loadings 0.7, and the loading of x1 is fixed to 1):

loading_x2 <- function(fit) {
  pe <- lavaan::parameterEstimates(fit)
  pe <- pe[pe$lhs == "F1" & pe$op == "=~" & pe$rhs == "x2", ]
  c(est = pe$est, se = pe$se, lower = pe$ci.lower, upper = pe$ci.upper)
}

res_ML  <- simulate_estimate(dists, model, n = 200, reps = 50,
                             extractor = loading_x2, seed = 2)
res_MLR <- simulate_estimate(dists, model, n = 200, reps = 50,
                             extractor = loading_x2, estimator = "MLR", seed = 2)

summarize <- function(r) {
  c(bias = mean(r$est) - 1,
    empirical_sd = sd(r$est),
    mean_se = mean(r$se),
    coverage = mean(r$lower < 1 & 1 < r$upper))
}
tab <- rbind(sapply(split(res_ML, res_ML$dist), summarize),
             sapply(split(res_MLR, res_MLR$dist), summarize))
rownames(tab) <- paste(rep(c("ML", "MLR"), each = 4), rownames(tab))
round(tab, 3)
#>                     IG PLSIM    VM
#> ML bias          0.058 0.072 0.012
#> ML empirical_sd  0.246 0.251 0.214
#> ML mean_se       0.154 0.157 0.147
#> ML coverage      0.820 0.820 0.820
#> MLR bias         0.058 0.072 0.012
#> MLR empirical_sd 0.246 0.251 0.214
#> MLR mean_se      0.203 0.239 0.200
#> MLR coverage     0.840 0.920 0.960

Both runs use seed = 2, so ML and MLR are computed on the same samples and differ only in the standard errors: the estimates are identical. When the mean standard error is far below the empirical standard deviation, the confidence intervals cover the true value less often than 95%.

Reproducibility and parallel computing

Each replication uses its own random number stream, determined by seed and its position in the design. The result is therefore the same for any number of cores, and the seed is stored with the result:

attr(res, "info")
#> $seed
#> [1] 1
#> 
#> $fit
#> [1] "lavaan"
#> 
#> $versions
#>  covsim  lavaan 
#> "1.2.0" "0.7.2"
res4 <- simulate_estimate(dists, model, n = c(200, 500), reps = 50,
                          extractor = chisq_tests, estimator = "MLR", seed = 1,
                          cores = 4)
all.equal(res, res4)   # TRUE (on Windows, cores is set to 1)

To keep the whole fitted objects, for analyses not planned in advance, use extractor = function(fit) list(fit = fit). This uses much more memory.

Estimation with magmaan

magmaan is a fast SEM estimator by Jonas Moss with a C++ core. It is not on CRAN; see the end of this vignette for installation. The code in this section is only run when magmaan (version 0.2.0 or later) is installed.

has_magmaan
#> [1] FALSE

With fit = "magmaan" the model is built once per distribution and reused for every sample. A magmaan fit has its own accessors, so the extractors are written differently:

Two differences from lavaan matter when results are compared. magmaan’s default standard errors are robust (sandwich), which corresponds to lavaan with estimator = "MLR", not to plain ML. And fit_measures(fit) reports robust RMSEA and CFI and a bias-corrected SRMR under the plain names; fit_measures(fit, lavaan_compat = "MLR") gives lavaan’s definitions and names.

The magmaan versions of the test and loading extractors:

mag_tests <- function(fit) {
  tests <- summary(fit)$tests
  setNames(tests$pvalue, paste(tests$test, tests$reference, sep = "_"))
}

res_mag <- simulate_estimate(dists, model, n = c(200, 500), reps = 50,
                             extractor = mag_tests, fit = "magmaan", seed = 1)
head(res_mag)
aggregate(cbind(score_peba4, lr_peba4) ~ dist + n, data = res_mag,
          FUN = reject)

Because the seed and the design are the same, res and res_mag are computed on exactly the same simulated samples, so the lavaan and magmaan tests can be compared replication by replication:

both <- merge(res[, c("dist", "n", "rep", "pvalue.scaled")],
              res_mag[, c("dist", "n", "rep", "score_peba4")])
aggregate(cbind(MLR = pvalue.scaled, PEBA4 = score_peba4) ~ dist, data = both,
          FUN = reject)
mag_loading_x2 <- function(fit) {
  cf <- summary(fit)$coefficients
  cf <- cf[cf$lhs == "F1" & cf$op == "=~" & cf$rhs == "x2", ]
  c(est = cf$est, se = cf$se, lower = cf$ci.lower, upper = cf$ci.upper)
}
mag_fit_indices <- function(fit) {
  fm <- magmaan::fit_measures(fit, lavaan_compat = "MLR")
  setNames(fm$estimate, fm$index)[c("rmsea", "rmsea.robust", "cfi", "cfi.robust")]
}

res_mag2 <- simulate_estimate(dists, model, n = 200, reps = 50, fit = "magmaan",
                              extractor = function(fit) c(mag_loading_x2(fit),
                                                          mag_fit_indices(fit)),
                              seed = 2)
round(sapply(split(res_mag2, res_mag2$dist), summarize), 3)

The extractor above combines two functions; the result has one column for each quantity. With the same seed, these numbers agree with the MLR rows of the lavaan table above.

Speed

magmaan is typically several times faster than lavaan, which matters in studies with many conditions:

t_lavaan <- system.time(
  simulate_estimate(dists$VM, model, n = 500, reps = 50,
                    extractor = chisq_tests, estimator = "MLR", seed = 1))[["elapsed"]]
t_magmaan <- system.time(
  simulate_estimate(dists$VM, model, n = 500, reps = 50,
                    extractor = mag_tests, fit = "magmaan", seed = 1))[["elapsed"]]
c(lavaan = t_lavaan, magmaan = t_magmaan, ratio = t_lavaan / t_magmaan)

Arguments in ... are passed to magmaan::magmaan(). For example, inference = FALSE skips standard errors and tests when only the estimates are needed.

Installing magmaan

magmaan needs a C++23 compiler and is installed from GitHub in two steps:

remotes::install_github("JonasMoss/magmaan", subdir = "r-package")
remotes::install_github("JonasMoss/magmaan", subdir = "r-magmaan")

Any estimator: a user-supplied fit function

fit can also be any function of the data and the model. Here lavaan is used with the WLS estimator, but the function could just as well call another package; the extractor must then match what that function returns:

fit_wls <- function(data, model) lavaan::sem(model, data = data, estimator = "WLS")
res_wls <- simulate_estimate(dists$IG, model, n = 500, reps = 20,
                             extractor = chisq_tests, fit = fit_wls, seed = 1)
reject(res_wls$pvalue)
#> [1] 0.1