Package {covsim}

Calibrate the IG (Independent Generator) model

Description

Performs the calibration step of the IG algorithm: computes the transformation matrix A and fits a Pearson distribution to each independent generator such that samples drawn from the calibrated model match the target covariance, skewness, and excess kurtosis. Calibration is decoupled from simulation so the same calibration object can be reused to draw many samples without recomputation.

Usage

calibrate_ig(sigma.target, skewness, excesskurtosis, typeA = "triang")

Arguments

Value

An object of class c("covsim_calib_ig", "covsim_calib") that can be passed to simulate to draw samples.

Author(s)

Njål Foldnes (njal.foldnes@gmail.com)

References

Foldnes, N. and Olson, U. H. (2016). A simple simulation technique for nonnormal data with prespecified skewness, kurtosis, and covariance matrix. Multivariate behavioral research, 51(2-3), 207-219

See Also

rIG, simulate

Examples

sigma.target <- diag(3)
cal <- calibrate_ig(sigma.target, skewness = rep(1, 3), excesskurtosis = rep(3, 3))
samples <- stats::simulate(cal, nsim = 5, N = 1000)

Calibrate the PLSIM (Piecewise Linear) model

Description

Performs the calibration step of the PLSIM algorithm: fits piecewise linear marginal transforms with the target skewness and kurtosis and searches for intermediate normal correlations that yield the target covariance after transformation. The (expensive) pairwise calibration is decoupled from simulation so the resulting calibration object can be reused to draw many samples without recomputation.

Usage

calibrate_plsim(
  sigma.target,
  skewness,
  excesskurtosis,
  numsegments = 4,
  gammalist = NULL,
  monot = FALSE,
  verbose = TRUE
)

Arguments

Value

An object of class c("covsim_calib_plsim", "covsim_calib") that can be passed to simulate to draw samples.

Author(s)

Njål Foldnes (njal.foldnes@gmail.com)

References

Foldnes, N. and Grønneberg S. (2021). Non-normal data simulation using piecewise linear transforms.

See Also

rPLSIM, simulate

Examples

set.seed(1)
sigma.target <- cov(MASS::mvrnorm(5, rep(0, 3), diag(3)))
cal <- calibrate_plsim(sigma.target, skewness = rep(1, 3), excesskurtosis = rep(4, 3))
samples <- stats::simulate(cal, nsim = 2, N = 1000)

Calibrate a regular vine (VITA)

Description

Wraps vita and returns a calibration object that can be passed to simulate. The returned object embeds the calibrated vine_dist so it can also be used directly with the rvinecopulib simulator if desired.

Usage

calibrate_vita(
  margins,
  sigma.target,
  vc = NULL,
  family_set = c("clayton", "gauss", "joe", "gumbel", "frank"),
  Nmax = 10^6,
  numrootpoints = 10,
  conflevel = 0.995,
  numpoints = 4,
  verbose = TRUE,
  cores = parallel::detectCores()
)

Arguments

Value

If a feasible solution was found, an object of class c("covsim_calib_vita", "covsim_calib") that can be passed to simulate; otherwise NULL.

See Also

vita, simulate

Examples

set.seed(1)
sigma.target <- cov(MASS::mvrnorm(10, mu = rep(0, 3), Sigma = diag(1, 3)))
marginsnorm <- lapply(X = sqrt(diag(sigma.target)),
                     function(X) list(distr = "norm", sd = X))
cal <- calibrate_vita(marginsnorm, sigma.target = sigma.target,
                      Nmax = 10^5, cores = 1)
samples <- stats::simulate(cal, nsim = 1, N = 1000)

Calibrate the VM (Vale-Maurelli) model

Description

Performs the calibration step of the Vale-Maurelli algorithm: for each variable the Fleishman polynomial a + bZ + cZ^2 + dZ^3 (with a = -c) matching the target skewness and excess kurtosis is found, and for each pair of variables the intermediate normal correlation that yields the target correlation after transformation is solved from the Vale-Maurelli cubic. Calibration is decoupled from simulation so the same calibration object can be reused to draw many samples without recomputation.

Usage

calibrate_vm(sigma.target, skewness, excesskurtosis)

Arguments

Details

When the Fleishman equations have several solutions, a monotone increasing polynomial is preferred. For many common skewness/kurtosis combinations (e.g. skewness 1 and excess kurtosis 1.5) only polynomials that turn down in the far tails exist; these are used, and flagged in monotone.

Value

An object of class c("covsim_calib_vm", "covsim_calib") that can be passed to simulate to draw samples. The element coef holds the Fleishman coefficients (one row per variable), monotone whether each polynomial is monotone increasing, and z.corr the intermediate correlation matrix.

References

Fleishman, A. I. (1978). A method for simulating non-normal distributions. Psychometrika, 43(4), 521-532.

Vale, C. D. and Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.

See Also

rVM, simulate

Examples

sigma.target <- matrix(c(1, 0.5, 0.5, 1), 2)
cal <- calibrate_vm(sigma.target, skewness = c(1, 2), excesskurtosis = c(2, 7))
samples <- stats::simulate(cal, nsim = 5, N = 1000)

Simulation of non-normal data

Description

Backward-compatible wrapper that combines calibrate_ig and simulate in one call. New code should prefer the calibrate/simulate split so calibration can be reused.

Usage

rIG(N, sigma.target, skewness, excesskurtosis, reps = 1, typeA = "triang")

Arguments

Value

A list of simulated samples.

Author(s)

Njål Foldnes (njal.foldnes@gmail.com)

References

Foldnes, N. and Olson, U. H. (2016). A simple simulation technique for nonnormal data with prespecified skewness, kurtosis, and covariance matrix. Multivariate behavioral research, 51(2-3), 207-219

Examples

set.seed(1234)
model <- '
 # measurement model
   ind60 =~ x1 + x2 + x3
   dem60 =~ y1 + y2 + y3 + y4
   dem65 =~ y5 + y6 + y7 + y8
 # regressions
   dem60 ~ ind60
   dem65 ~ ind60 + dem60
 # residual correlations
   y1 ~~ y5
   y2 ~~ y4 + y6
   y3 ~~ y7
   y4 ~~ y8
   y6 ~~ y8'
fit  <- lavaan::sem(model, data=lavaan::PoliticalDemocracy)
population.sigma <- lavaan::lavInspect(fit, "sigma.hat")
population.skew  <- c(0, 0, 0, 0, 1, 1, 1, 1, 2,2,2 )
population.excesskurt <- c( 1 , 1, 1, 1, 3, 3, 3, 3, 15, 15, 15)
my.samples <- rIG(N=10^3, sigma=population.sigma,
        skewness=population.skew,
        excesskurt=population.excesskurt,
        reps=5)

Simulation of non-normal data

Description

Backward-compatible wrapper that combines calibrate_plsim and simulate in one call. New code should prefer the calibrate/simulate split so the (expensive) pairwise calibration can be reused.

Usage

rPLSIM(
  N,
  sigma.target,
  skewness,
  excesskurtosis,
  reps = 1,
  numsegments = 4,
  gammalist = NULL,
  monot = FALSE,
  verbose = TRUE
)

Arguments

Value

A list with two elements. First element: the list of simulated samples. Second element: the fitted piecewise linear functions and the intermediate correlations matrix.

Author(s)

Njål Foldnes (njal.foldnes@gmail.com)

References

Foldnes, N. and Grønneberg S. (2021). Non-normal data simulation using piecewise linear transforms.Under review.

Examples

set.seed(1)
sigma.target  <- cov(MASS::mvrnorm(5, rep(0,3), diag(3)))
res  <- covsim::rPLSIM(10^5, sigma.target, skewness=rep(1,3), excesskurtosis=rep(4,3))
my.sample  <- res[[1]][[1]]

Simulation of non-normal data with the Vale-Maurelli method

Description

Convenience wrapper that combines calibrate_vm and simulate in one call. Prefer the calibrate/simulate split when drawing many samples from the same target.

Usage

rVM(N, sigma.target, skewness, excesskurtosis, reps = 1)

Arguments

Value

A list of simulated samples.

References

Vale, C. D. and Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.

Examples

sigma.target <- matrix(c(1, 0.5, 0.5, 1), 2)
my.samples <- rVM(N = 1000, sigma.target, skewness = c(1, 2),
                  excesskurtosis = c(2, 7), reps = 5)

Simulate from a calibrated IG model

Description

S3 method for simulate that draws samples from a covsim_calib_ig object produced by calibrate_ig.

Usage

## S3 method for class 'covsim_calib_ig'
simulate(object, nsim = 1, seed = NULL, N, ...)

Arguments

Value

A list of nsim matrices, each with N rows.

Simulate from a calibrated PLSIM model

Description

S3 method for simulate that draws samples from a covsim_calib_plsim object produced by calibrate_plsim.

Usage

## S3 method for class 'covsim_calib_plsim'
simulate(object, nsim = 1, seed = NULL, N, ...)

Arguments

Value

A list of nsim matrices, each with N rows and the column names of sigma.target.

Simulate from a calibrated VITA vine

Description

S3 method for simulate that draws samples from a covsim_calib_vita object produced by calibrate_vita.

Usage

## S3 method for class 'covsim_calib_vita'
simulate(object, nsim = 1, seed = NULL, N, cores = 1, ...)

Arguments

Value

A list of nsim data frames (or matrices), each with N rows.

Simulate from a calibrated VM model

Description

S3 method for simulate that draws samples from a covsim_calib_vm object produced by calibrate_vm.

Usage

## S3 method for class 'covsim_calib_vm'
simulate(object, nsim = 1, seed = NULL, N, ...)

Arguments

Value

A list of nsim matrices, each with N rows.

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

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")])
}

Calibrate a regular vine

Description

vita implements the VITA (VIne-To-Anything) algorithm. Covariance matrix and margins are specified, and vita calibrates the pair-copulas in each node of the tree to match the target covariance.

Usage

vita(
  margins,
  sigma.target,
  vc = NULL,
  family_set = c("clayton", "gauss", "joe", "gumbel", "frank"),
  Nmax = 10^6,
  numrootpoints = 10,
  conflevel = 0.995,
  numpoints = 4,
  verbose = TRUE,
  cores = parallel::detectCores()
)

Arguments

Value

If a feasible solution was found, a vine to be used for simulation

References

Grønneberg, S., Foldnes, N., & Marcoulides, K. M. (2021). covsim: An r package for simulating non-normal data for structural equation models using copulas. Journal of Statistical Software. doi:10.18637/jss.v102.i03

Examples

set.seed(1)# define a target covariance. 3 dimensions.
sigma.target <- cov(MASS::mvrnorm(10, mu=rep(0,3), Sigma=diag(1, 3)))

#normal margins that match the covariances:
marginsnorm <- lapply(X=sqrt(diag(sigma.target)),function(X) list(distr="norm", sd=X) )

#calibrate with a default D-vine, with rather low precision (default Nmax is 10^6)
# if cores=1 is removed, all cores are used, with a speed gain
calibrated.vine <- vita(marginsnorm, sigma.target =sigma.target, Nmax=10^5, cores=1)
#check
#round(cov(rvinecopulib::rvine(10^5, calibrated.vine))-sigma.target, 3)

#margins are normal but dependence structure is not
#pairs(rvinecopulib::rvine(500, calibrated.vine))