---
title: "Monte Carlo studies with covsim"
author: "Njål Foldnes"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Monte Carlo studies with covsim}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE
)
has_magmaan <- requireNamespace("magmaan", quietly = TRUE) &&
  utils::packageVersion("magmaan") >= "0.2.0"
```

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.

```{r}
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.

```{r}
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:

```{r}
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)
```

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

```{r, eval = FALSE}
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:

```{r}
samples <- simulate(dists$PLSIM, nsim = 50, N = 200, seed = 1)
length(samples)
head(samples[[1]], 3)
```

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()`:

```{r}
set.seed(1)
old_way <- rPLSIM(200, sigma, skew, kurt, reps = 50, verbose = FALSE)
length(old_way$samples)
```

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:

```{r}
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)
```

Rejection rates at the 5% level:

```{r}
round(sapply(by_hand, function(p) colMeans(p < 0.05, na.rm = TRUE)), 2)
```

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:

```{r}
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)
```

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.

```{r}
table(res$converged, useNA = "ifany")
```

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

```{r}
reject <- function(p) mean(p < 0.05, na.rm = TRUE)
aggregate(cbind(ML = pvalue, MLR = pvalue.scaled) ~ dist + n, data = res,
          FUN = reject)
```

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:

```{r}
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)
```

#### 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):

```{r}
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)
```

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:

```{r}
attr(res, "info")
```

```{r, eval = FALSE}
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](https://github.com/JonasMoss/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.

```{r}
has_magmaan
```

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:

- `summary(fit)$tests`: the robust score and likelihood ratio tests,
  calibrated with the PEBA4 approximation by default;
- `summary(fit)$coefficients`: the parameter table, with columns `lhs`,
  `op`, `rhs`, `est`, `se`, `pvalue`, `ci.lower` and `ci.upper`;
- `magmaan::fit_measures(fit)`: the fit indices.

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:

```{r, eval = has_magmaan}
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:

```{r, eval = has_magmaan}
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)
```

```{r, eval = has_magmaan}
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:

```{r, eval = has_magmaan}
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:

```{r, eval = FALSE}
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:

```{r}
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)
```