Package 'simsem'

Description

Data analysis using the model specification (linkS4class{SimSem}) or the mx model object (MxModel). Data will be multiply imputed if the miss argument is specified.

Usage

analyze(model, data, package="lavaan", miss=NULL, aux=NULL, group = NULL, 
	mxMixture = FALSE, ...)

Arguments

Value

The lavaan object containing the output

Author(s)

Patrick Miller (University of Notre Dame; [email protected]), Sunthud Pornprasertmanit ([email protected])

See Also

Note that users can use functions provided by lavaan package (lavaan, cfa, sem, or growth) if they wish to analyze data by lavaan directly.

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
LY <- bind(loading, 0.7)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VY <- bind(rep(NA,6),2)

CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType = "CFA")

dat <- generate(CFA.Model,200)
out <- analyze(CFA.Model,dat)

Description

This function will provide averages of model fit statistics and indices for nested models. It will also provide average differences of fit indices and power for likelihood ratio tests of nested models.

Arguments

Value

A data frame that provides the statistics described above from all parameters. For using with linkS4class{SimResult}, the result is a list with two or three elements:

• summary: Average of fit indices across all replications

• diff: Average of the differences in fit indices across all replications

• varyParam: The statistical power of chi-square difference test given values of varying parameters (such as sample size or percent missing)

Author(s)

Alexander M. Schoemann (East Carolina University; [email protected]), Sunthud Pornprasertmanit ([email protected])

See Also

SimResult for the object input

Examples

## Not run: 
loading1 <- matrix(0, 6, 1)
loading1[1:6, 1] <- NA
loading2 <- loading1
loading2[6,1] <- 0
LY1 <- bind(loading1, 0.7)
LY2 <- bind(loading2, 0.7)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model1 <- model(LY = LY1, RPS = RPS, RTE = RTE, modelType="CFA")
CFA.Model2 <- model(LY = LY2, RPS = RPS, RTE = RTE, modelType="CFA")

# We make the examples running only 5 replications to save time.
# In reality, more replications are needed.
# Need to make sure that both simResult calls have the same seed!
Output1 <- sim(5, n=500, model=CFA.Model1, generate=CFA.Model1, seed=123567)
Output2 <- sim(5, n=500, model=CFA.Model2, generate=CFA.Model1, seed=123567)
anova(Output1, Output2)

# The example when the sample size is varying
Output1b <- sim(NULL, n=seq(50, 500, 50), model=CFA.Model1, generate=CFA.Model1, seed=123567)
Output2b <- sim(NULL, n=seq(50, 500, 50), model=CFA.Model2, generate=CFA.Model1, seed=123567)
anova(Output1b, Output2b)

## End(Not run)

Description

Create SimMatrix or SimVector object that specifies

1. Pattern of fixed/freed parameters for analysis

2. Population parameter values for data generation

3. Any model misspecification (true population parameter is different than the one specified) for these parameters.

Each matrix in the Lisrel-style notation is specified in this way (e.g. LY, PS, and TE) and is used to create a model analysis template and a data generation template for simulation through the model function.

Usage

bind(free = NULL, popParam = NULL, misspec = NULL, symmetric = FALSE)
binds(free = NULL, popParam = NULL, misspec = NULL, symmetric = TRUE)

Arguments

Details

Bind is the first step in the bind -> model -> sim workflow of simsem, and this document outlines the user interface or language used to describe these simulations. This interface, while complex, enables a wide array of simulation specifications for structural equation models by building on LISREL-style parameter specifications.

In simulations supported by simsem, a given parameter may be either fixed or freed for analysis, but may optionally also have a population value or distribution for data generation, or a value or distribution of misspecification. The purpose of bind is to stack these multiple meanings of a parameter into an object recognized by simsem, a SimMatrix. Each matrix in the Lisrel notation (e.g. LY, PS, TE, BE) becomes a SimMatrix, and is passed to the function model, which builds the data generation template and an analysis template (a lavaan parameter table), collectively forming a SimSem object, which can be passed to the function sim for simulation.

Note that any (dim)names attributes will be set to NULL for any vectors or matrices passed to free, popParam, or misspec in order to prevent errors elsewhere in the workflow. To set custom variable names, please use any of the indLab, facLab, covLab, or groupLab arguments to model().

Value

SimMatrix or SimVector object that used for model specification for analysis and data generation in simsem.

Author(s)

Patrick Miller (University of Notre Dame; [email protected]), Sunthud Pornprasertmanit ([email protected])

See Also

• model To combine simMatrix objects into a complete data analysis and data generation template, which is a SimSem object

• generate To generate data using the simsem template.

• analyze To analyze real or generated data using the simsem template.

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
loadingValues <- matrix(0, 6, 2)
loadingValues[1:3, 1] <- 0.7
loadingValues[4:6, 2] <- 0.7
LY <- bind(loading, loadingValues)
summary(LY)

# Set both factor correlations to .05
latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

# Misspecify all error covarainces
error.cor <- matrix(0, 6, 6)
diag(error.cor) <- NA
RTE <- binds(error.cor,1,"runif(1,-.05,.05)")

Description

Create a data distribution object. There are two ways to specify nonnormal data-generation model. To create nonnormal data by the copula method, margins and ... arguments are required. To create data by Vale and Maurelli's method, skewness and/or kurtosis arguments are required.

Usage

bindDist(margins = NULL, ..., p = NULL, keepScale = TRUE, reverse = FALSE, 
	copula = NULL, skewness = NULL, kurtosis = NULL)

Arguments

Value

SimDataDist that saves analysis result from simulate data.

Author(s)

Sunthud Pornprasertmanit ([email protected])

References

Mair, P., Satorra, A., & Bentler, P. M. (2012). Generating nonnormal multivariate data using copulas: Applications to SEM. Multivariate Behavioral Research, 47, 547-565.

Vale, C. D. & Maurelli, V. A. (1983) Simulating multivariate nonormal distributions. Psychometrika, 48, 465-471.

See Also

• SimResult for the type of resulting object

Examples

# Create data based on Vale and Maurelli's method by specifying skewness and kurtosis
dist <- bindDist(skewness = c(0, -2, 2), kurtosis = c(0, 8, 4))

## Not run: 
library(copula)

# Create three-dimensional distribution by gaussian copula with 
# the following marginal distributions
#   1. t-distribution with df = 2
# 	2. chi-square distribution with df = 3
#	3. normal distribution with mean = 0 and sd = 1

# Setting the attribute of each marginal distribution
d1 <- list(df=2)
d2 <- list(df=3)
d3 <- list(mean=0, sd=1)

# Create a data distribution object by setting the names of each distribution
# and their arguments
dist <- bindDist(c("t", "chisq", "norm"), d1, d2, d3)

# Create data by using Gumbel Copula as the multivariate distribution
dist <- bindDist(c("t", "chisq", "norm"), d1, d2, d3, copula = gumbelCopula(2, dim = 3))

# Reverse the direction of chi-square distribution from positively skew to negatively skew
dist <- bindDist(c("t", "chisq", "norm"), d1, d2, d3, copula = gumbelCopula(2, dim = 3),
	reverse = c(FALSE, TRUE, FALSE))

## End(Not run)

Description

Extract parameter estimates from a simulation result. This function is similar to the inspect method with what = "coef".

Arguments

Value

Parameter estimates of each replication

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

SimResult for the object input

Examples

## Not run: 
loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
LY <- bind(loading, 0.7)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VY <- bind(rep(NA,6),2)

CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType = "CFA")

# In reality, more than 5 replications are needed.
Output <- sim(5, CFA.Model, n=200)
coef(Output)
coef(Output, improper = TRUE)

## End(Not run)

Description

Combine result objects into a single result object

Usage

combineSim(...)

Arguments

Value

A combined result object

Author(s)

Terry Jorgensen (University of Kansas; [email protected]), Sunthud Pornprasertmanit ([email protected])

See Also

Result object (SimResult)

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
LY <- bind(loading, 0.7)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VY <- bind(rep(NA,6),2)

CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType = "CFA")
Output1 <- sim(5, CFA.Model, n=200, seed=123321)
Output2 <- sim(4, CFA.Model, n=200, seed=324567)
Output3 <- sim(3, CFA.Model, n=200, seed=789987)
Output <- combineSim(Output1, Output2, Output3)
summary(Output)

Description

A function to find the coverage rate of confidence intervals in a model when one or more of the simulations parameters vary randomly across replications.

Usage

continuousCoverage(simResult, coverValue = NULL, contN = TRUE, contMCAR = FALSE, 
    contMAR = FALSE, contParam = NULL, coverParam = NULL, pred = NULL)

Arguments

Details

In this function, the coverage (which can be 0 or 1) is regressed on randomly varying simulation parameters (e.g., sample size, percentage of missing data, or model parameters) using logistic regression. For a set of independent variables values, the predicted probability from the logistic regression equation is the predicted coverage rate.

Value

Data frame containing columns representing values of the randomly varying simulation parameters, and coverage rates for model parameters of interest.

Author(s)

Sunthud Pornprasertmanit ([email protected]), Alexander M. Schoemann (East Carolina University; [email protected])

See Also

• SimResult to see how to create a simResult object with randomly varying parameters.

Examples

## Not run: 
# Specify Sample Size by n
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, 0.7)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# Specify both continuous sample size and percent missing completely at random. 
# Note that more fine-grained values of n and pmMCAR is needed, e.g., n=seq(50, 500, 1) 
# and pmMCAR=seq(0, 0.2, 0.01)
Output <- sim(NULL, CFA.Model, n=seq(100, 200, 20), pmMCAR=c(0, 0.1, 0.2))
summary(Output)

# Find the coverage rates of all combinations of different sample size and percent MCAR missing
Ccover <- continuousCoverage(Output, contN = TRUE, contMCAR = TRUE)
Ccover

# Find the coverage rates of parameter estimates when sample size is 200 
# and percent MCAR missing is 0.3
Ccover2 <- continuousCoverage(Output, coverValue=0, contN = TRUE, contMCAR = TRUE, 
     pred=list(N = 200, pmMCAR = 0.3))
Ccover2

## End(Not run)

Description

A function to find the power of parameters in a model when one or more of the simulations parameters vary randomly across replications.

Usage

continuousPower(simResult, contN = TRUE, contMCAR = FALSE, contMAR = FALSE, 
	contParam = NULL, alpha = .05, powerParam = NULL, pred = NULL)

Arguments

Details

A common use of simulations is to conduct power analyses, especially when using SEM (Muthen & Muthen, 2002). Here, researchers select values for each parameter and a sample size and run a simulation to determine power in those conditions (the proportion of generated datasets in which a particular parameter of interest is significantly different from zero). To evaluate power at multiple sample sizes, one simulation for each sample size must be run. By continuously varying sample size across replications, only a single simulation is needed. In this simulation, the sample size for each replication varies randomly across plausible sample sizes (e.g., sample sizes between 200 and 500). For each replication, the sample size and significance of each parameter (0 = not significant, 1 = significant) are recorded. When the simulation is complete, parameter significance is regressed on sample size using logistic regression. For a given sample size, the predicted probability from the logistic regression equation is the power to detect an effect at that sample size. This approach can be extended to other randomly varying simulation parameters such as the percentage of missing data, and model parameters.

Value

Data frame containing columns representing values of the randomly varying simulation parameters, and power for model parameters of interest.

Author(s)

Alexander M. Schoemann (East Carolina University; [email protected]), Sunthud Pornprasertmanit ([email protected])

References

Muthen, L. K., & Muthen, B. O. (2002). How to use a Monte Carlo study to decide on sample size and determine power. Structural Equation Modeling, 4, 599-620.

See Also

• SimResult to see how to create a simResult object with randomly varying parameters.

Examples

## Not run: 
# Specify Sample Size by n
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, 0.7)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")
dat <- generate(CFA.Model, 50)
out <- analyze(CFA.Model, dat)

# Specify both continuous sample size and percent missing completely at random. 
# Note that more fine-grained values of n and pmMCAR is needed, e.g., n=seq(50, 500, 1) 
# and pmMCAR=seq(0, 0.2, 0.01)
Output <- sim(NULL, CFA.Model, n=seq(100, 200, 20), pmMCAR=c(0, 0.1, 0.2))
summary(Output)

# Find the power of all combinations of different sample size and percent MCAR missing
Cpow <- continuousPower(Output, contN = TRUE, contMCAR = TRUE)
Cpow

# Find the power of parameter estimates when sample size is 200 and percent MCAR missing is 0.3
Cpow2 <- continuousPower(Output, contN = TRUE, contMCAR = TRUE, pred=list(N = 200, pmMCAR = 0.3))
Cpow2

## End(Not run)

Description

This function can be used to create data from a set of parameters created from draw, called a codeparamSet. This function is used internally to create data, and is available publicly for accessibility and debugging.

Usage

createData(paramSet, n, indDist=NULL, sequential=FALSE, facDist=NULL, 
errorDist=NULL, saveLatentVar = FALSE, indLab=NULL, facLab = NULL, 
modelBoot=FALSE, realData=NULL, covData=NULL, empirical = FALSE)

Arguments

Details

This function will use the modified mvrnorm function (from the MASS package) by Paul E. Johnson to create data from model implied covariance matrix if the data distribution object (SimDataDist) is not specified. The modified function is just a small modification from the original mvrnorm function such that the data generated with the sample sizes of n and n + k (where k > 0) will be replicable in the first n rows.

It the data distribution object is specified, either the copula model or the Vale and Maurelli's method is used. For the copula approach, if the copula argument is not specified in the data distribution object, the naive Gaussian copula is used. The correlation matrix is direct applied to the multivariate Gaussian copula. The correlation matrix will be equivalent to the Spearman's correlation (rank correlation) of the resulting data. If the copula argument is specified, such as ellipCopula, normalCopula, or archmCopula, the data-transformation method from Mair, Satorra, and Bentler (2012) is used. In brief, the data ($X$) are created from the multivariate copula. The covariance from the generated data is used as the starting point ($S$). Then, the target data ($Y$) with the target covariance as model-implied covariance matrix ($\Sigma_0$) can be created:

$Y = XS^{-1/2}\Sigma^{1/2}_0.$

See bindDist for further details. For the Vale and Maurelli's (1983) method, the code is brought from the lavaan package.

For the model-based bootstrap, the transformation proposed by Yung & Bentler (1996) is used. This procedure is the expansion from the Bollen and Stine (1992) bootstrap including a mean structure. The model-implied mean vector and covariance matrix with trivial misspecification will be used in the model-based bootstrap if misspec is specified. See page 133 of Bollen and Stine (1992) for a reference.

Internally, parameters are first drawn, and data is then created from these parameters. Both of these steps are available via the draw and createData functions respectively.

Value

A data.frame containing simulated data from the data generation template. A variable "group" is appended indicating group membership.

Author(s)

Sunthud Pornprasertmanit ([email protected]), Patrick Miller (University of Notre Dame; [email protected]). The original code of mvrnorm function is based on the MASS package slightly modified by Paul E. Johnson. The code for data-transformation in multivariate copula is based on Mair et al. (2012) article. The code for Vale and Maurelli (1983) is slightly modified from the function provided in the lavaan package.

References

Bollen, K. A., & Stine, R. A. (1992). Bootstrapping goodness-of-fit measures in structural equation models. Sociological Methods and Research, 21, 205-229.

Mair, P., Satorra, A., & Bentler, P. M. (2012). Generating nonnormal multivariate data using copulas: Applications to SEM. Multivariate Behavioral Research, 47, 547-565.

Vale, C. D. & Maurelli, V. A. (1983) Simulating multivariate nonormal distributions. Psychometrika, 48, 465-471.

Yung, Y.-F., & Bentler, P. M. (1996). Bootstrapping techniques in analysis of mean and covariance structures. In G. A. Marcoulides & R. E. Schumacker (Eds.), Advanced structural equation modeling: Issues and techniques (pp. 195-226). Mahwah, NJ: Erlbaum.

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
LY <- bind(loading, 0.7)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VY <- bind(rep(NA,6),2)

CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType = "CFA")

# Draw a parameter set for data generation.
param <- draw(CFA.Model)

# Generate data from the first group in the paramList.
dat <- createData(param[[1]], n = 200)

Draw parameters from a SimSem object.

Description

This function draws parameters from a SimSem template, for debugging or other use. Used internally to create data. Data can be created in one step from a SimSem object using generate.

Usage

draw(model, maxDraw=50, misfitBounds=NULL, averageNumMisspec=FALSE, 
optMisfit = NULL, optDraws = 50, misfitType = "f0", createOrder = c(1, 2, 3),
covData = NULL)

Arguments

Value

Nested list of drawn parameters in the form [[Group]][[param,misspec,misOnly]][[SimMatrix]]. E.g. The LY parameter matrix of the first group would be indexed as obj[[1]]$param$LY. The values in $param are the raw parameter values with no misspecification. The values in $misspec are raw parameter values + misspecification. The values in $misOnly are only the misspecification values.

Author(s)

Sunthud Pornprasertmanit ([email protected]), Patrick Miller (University of Notre Dame; [email protected])

See Also

createData To generate random data using a set of parameters from draw

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
LY <- bind(loading, 0.7)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VY <- bind(rep(NA,6),2)

CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType = "CFA")

# Draw a parameter set for data generation.
param <- draw(CFA.Model)

# Example of Multiple Group Model with Weak Invariance

loading.in <- matrix(0, 6, 2)
loading.in[1:3, 1] <- c("load1", "load2", "load3")
loading.in[4:6, 2] <- c("load4", "load5", "load6")
mis <- matrix(0,6,2)
mis[loading.in == "0"] <- "runif(1, -0.1, 0.1)"
LY.in <- bind(loading.in, "runif(1, 0.7, 0.8)", mis)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VTE <- bind(rep(NA, 6), 0.51)

VPS1 <- bind(rep(1, 2))

VPS2 <- bind(rep(NA, 2), c(1.1, 1.2))

# Inequality constraint
script <- "
sth := load1 + load2 + load3
load4 == (load5 + load6) / 2
load4 > 0
load5 > 0
sth2 := load1 - load2
"

# Model Template
weak <- model(LY = LY.in, RPS = RPS, VPS=list(VPS1, VPS2), RTE = RTE, VTE=VTE, ngroups=2, 
     modelType = "CFA", con=script)

# Constraint --> Misspecification --> Fill Parameters
draw(weak, createOrder=c(1, 2, 3))

# Constraint --> Fill Parameters --> Misspecification 
draw(weak, createOrder=c(1, 3, 2))

# Misspecification --> Constraint --> Fill Parameters
draw(weak, createOrder=c(2, 1, 3))

# Misspecification --> Fill Parameters --> Constraint
draw(weak, createOrder=c(2, 3, 1))

# Fill Parameters --> Constraint --> Misspecification
draw(weak, createOrder=c(3, 1, 2))

# Fill Parameters --> Misspecification --> Constraint
draw(weak, createOrder=c(3, 2, 1))

Description

Creates a data analysis template (lavaan parameter table) for simulations with structural equation models based on Y-side LISREL design matrices. Each corresponds to a LISREL matrix, but must be a matrix or a vector. In addition to the usual Y-side matrices in LISREL, both PS and TE can be specified using correlations (RPS, RTE) and scaled by a vector of residual variances (VTE, VPS) or total variances (VY, VE). Multiple groups are supported by passing lists of matrices or vectors to arguments, or by specifying the number of groups.

Usage

estmodel(LY = NULL, PS = NULL, RPS = NULL, TE = NULL, RTE = NULL, BE = NULL, 
	VTE = NULL, VY = NULL, VPS = NULL, VE = NULL, TY = NULL, AL = NULL, 
	MY = NULL, ME = NULL, KA = NULL, GA = NULL, modelType, indLab = NULL, 
	facLab = NULL, covLab = NULL, groupLab = "group", ngroups = 1, con = NULL)
estmodel.cfa(LY = NULL,PS = NULL,RPS = NULL, TE = NULL,RTE = NULL, VTE = NULL, 
	VY = NULL, VPS = NULL, VE = NULL, TY = NULL, AL = NULL, MY = NULL, ME = NULL, 
	KA = NULL, GA = NULL, indLab = NULL, facLab = NULL, covLab = NULL, 
	groupLab = "group", ngroups = 1, con = NULL)
estmodel.path(PS = NULL, RPS = NULL, BE = NULL, VPS = NULL, VE = NULL, AL = NULL, 
	ME = NULL, KA = NULL, GA = NULL, indLab = NULL, facLab = NULL, covLab = NULL, 
	groupLab = "group", ngroups = 1,con = NULL)
estmodel.sem(LY = NULL,PS = NULL,RPS = NULL, TE = NULL,RTE = NULL, BE = NULL, 
	VTE = NULL, VY = NULL, VPS = NULL, VE = NULL, TY = NULL, AL = NULL, MY = NULL, 
	ME = NULL, KA = NULL, GA = NULL, indLab = NULL, facLab = NULL, covLab = NULL, 
	groupLab = "group", ngroups = 1, con = NULL)

Arguments

Details

This function contains default settings:

For modelType="CFA", LY is required. As the default, the on-diagonal elements of PS are fixed as 1 and the off-diagonal elements of PS are freely estimated. The off-diagonal elements of TE are freely estimated and the off-diagonal elements of TE are fixed to 0. The AL elements are fixed to 0. The TY elements are freely estimated.

For modelType="Path", BE is required. As the default, the on-diagonal elements of PS are freely estimated, the off-diagonal elements between exogenous variables (covariance between exogenous variables) are freely estimated, and the other off-diagonal elements are fixed to 0. The AL elements are freely estimated.

For modelType="SEM", LY and BE are required. As the default, the on-diagonal elements of PS are fixed to 1, the off-diagonal elements between exogenous factors (covariance between exogenous factors) are freely estimated, and the other off-diagonal elements are fixed to 0. The off-diagonal elements of TE are freely estimated and the off-diagonal elements of TE are fixed to 0. The AL elements are fixed to 0. The TY elements are freely estimated.

The estmodel.cfa, estmodel.path, and estmodel.sem are the shortcuts for the estmodel function when modelType are "CFA", "Path", and "SEM", respectively.

Value

SimSem object that contains the data generation template (@dgen) and analysis template (@pt).

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• model To build data generation and data analysis template for simulation.

• sim For simulations using the SimSem template.

• generate To generate data using the SimSem template.

• analyze To analyze real or generated data using the SimSem template.

• draw To draw parameters using the SimSem template.

Examples

loading <- matrix(0, 12, 4)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
loading[7:9, 3] <- NA
loading[10:12, 4] <- NA

CFA.Model <- estmodel(LY = loading, modelType = "CFA")

path <- matrix(0, 4, 4)
path[3, 1:2] <- NA
path[4, 3] <- NA
Path.Model <- estmodel(BE = path, modelType = "Path")

SEM.Model <- estmodel(BE = path, LY = loading, modelType="SEM")

# Shortcut
CFA.Model <- estmodel.cfa(LY = loading)
Path.Model <- estmodel.path(BE = path)
SEM.Model <- estmodel.sem(BE = path, LY = loading)

Description

This function can be used to export data created from a set of parameters created from draw, called a codeparamSet. This function can export data to be analyzed with either Mplus or LISREL.

Usage

exportData(nRep, model, n, program = "Mplus", fileStem = "sim", miss = NULL, 
	missCode = -999, datafun=NULL, pmMCAR = NULL, pmMAR = NULL, facDist = NULL, 
	indDist = NULL, errorDist = NULL, sequential = FALSE, modelBoot = FALSE, 
	realData = NULL, maxDraw = 50, misfitType = "f0", misfitBounds = NULL, 
	averageNumMisspec = NULL, optMisfit=NULL, optDraws = 50, seed = 123321, 
	silent = FALSE, multicore = FALSE, numProc = NULL,  params = FALSE)

Arguments

Value

Text files saved to the current working directory. If program = "Mplus" one file is output for each replication, and an extra file is output with the names of all saved data sets (this file can be used with the MONTECARLO command in Mplus). If program = "LISREL" one file is output with each replication stacked on top of the next (this file can be used with the RP command in LISREL). If program = TRUE, a list of parameter values for each replication is returned.

Author(s)

Alexander M. Schoemann (East Carolina University; [email protected])

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
LY <- bind(loading, 0.7)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VY <- bind(rep(NA,6),2)

CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType = "CFA")

## Export 20 replications to an external data file (not run).
#exportData(20, CFA.Model, 200)

Description

Find a value of independent variable that provides a given value of coverage rate. If there are more than one varying parameters, this function will find the value of the target varying parameters given the values of the other varying parameters.

Usage

findCoverage(coverTable, iv, target)

Arguments

Value

There are five possible types of values provided:

• Value The varying parameter value that provides the coverage rate just under the specified coverage rate (the adjacent value of varying parameter provides over power than the specified power value).

• Minimum value The minimum value has already provided the low coverage rate (way under the specified coverage rate). The value of varying parameters that provides exact coverage rate may be lower than the minimum value. The example of varying parameter that can provides the minimum value is sample size.

• Maximum value The maximum value has already provided the low coverage rate (way under the specified coverage rate). The value of varying parameters that provides exact desired power may be higher than the maximum value. The example of varying parameter that can provides the maximum value is percent missing.

• NA There is no value in the domain of varying parameters that provides the coverage rate lower than the desired coverage rate.

• Inf The coverage rate of all values in the varying parameters is 0 (specifically more than 0.0001) and any values of the varying parameters can be picked and still provide enough power.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• getCoverage to find the coverage rate of parameter estimates

• continuousCoverage to find the coverage rate of parameter estimates for the result object (linkS4class{SimResult}) with varying parameters.

Examples

## Not run: 
# Specify Sample Size by n
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, 0.4)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# Specify both sample size and percent missing completely at random. Note that more fine-grained 
# values of n and pmMCAR is needed, e.g., n=seq(50, 500, 1) and pmMCAR=seq(0, 0.2, 0.01)
Output <- sim(NULL, model=CFA.Model, n=seq(100, 200, 20), pmMCAR=c(0, 0.1, 0.2))

# Find the power of all possible combination of N and pmMCAR
cover <- getCoverage(Output, coverValue = 0)

# Find the sample size that provides the power of 0.8
findCoverage(cover, "N", 0.20)

## End(Not run)

Description

Find factor intercept from regression coefficient matrix and factor total means for latent variable models. In the path analysis model, this function will find indicator intercept from regression coefficient and indicator total means.

Usage

findFactorIntercept(beta, factorMean = NULL, gamma = NULL, covmean = NULL)

Arguments

Value

A vector of factor (indicator) intercepts

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndMean to find indicator (measurement) total means

• findIndResidualVar to find indicator (measurement) residual variances

• findIndTotalVar to find indicator (measurement) total variances

• findFactorMean to find factor means

• findFactorResidualVar to find factor residual variances

• findFactorTotalVar to find factor total variances

• findFactorTotalCov to find factor covariances

Examples

path <- matrix(0, 9, 9)
path[4, 1] <- path[7, 4] <- 0.6
path[5, 2] <- path[8, 5] <- 0.6
path[6, 3] <- path[9, 6] <- 0.6
path[5, 1] <- path[8, 4] <- 0.4
path[6, 2] <- path[9, 5] <- 0.4
factorMean <- c(5, 2, 3, 0, 0, 0, 0, 0, 0)
findFactorIntercept(path, factorMean)

Description

Find factor total means from regression coefficient matrix and factor intercepts for latent variable models. In the path analysis model, this function will find indicator total means from regression coefficient and indicator intercept.

Usage

findFactorMean(beta, alpha = NULL, gamma = NULL, covmean = NULL)

Arguments

Value

A vector of factor (indicator) total means

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndMean to find indicator (measurement) total means

• findIndResidualVar to find indicator (measurement) residual variances

• findIndTotalVar to find indicator (measurement) total variances

• findFactorIntercept to find factor intercepts

• findFactorResidualVar to find factor residual variances

• findFactorTotalVar to find factor total variances

• findFactorTotalCov to find factor covariances

Examples

path <- matrix(0, 9, 9)
path[4, 1] <- path[7, 4] <- 0.6
path[5, 2] <- path[8, 5] <- 0.6
path[6, 3] <- path[9, 6] <- 0.6
path[5, 1] <- path[8, 4] <- 0.4
path[6, 2] <- path[9, 5] <- 0.4
intcept <- c(5, 2, 3, 0, 0, 0, 0, 0, 0)
findFactorMean(path, intcept)

Description

Find factor residual variances from regression coefficient matrix, factor (residual) correlation matrix, and total factor variances for latent variable models. In the path analysis model, this function will find indicator residual variances from regression coefficient, indicator (residual) correlation matrix, and total indicator variances.

Usage

findFactorResidualVar(beta, corPsi, totalVarPsi = NULL, gamma = NULL, covcov = NULL)

Arguments

Value

A vector of factor (indicator) residual variances

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndMean to find indicator (measurement) total means

• findIndResidualVar to find indicator (measurement) residual variances

• findIndTotalVar to find indicator (measurement) total variances

• findFactorIntercept to find factor intercepts

• findFactorMean to find factor means

• findFactorTotalVar to find factor total variances

• findFactorTotalCov to find factor covariances

Examples

path <- matrix(0, 9, 9)
path[4, 1] <- path[7, 4] <- 0.6
path[5, 2] <- path[8, 5] <- 0.6
path[6, 3] <- path[9, 6] <- 0.6
path[5, 1] <- path[8, 4] <- 0.4
path[6, 2] <- path[9, 5] <- 0.4
facCor <- diag(9)
facCor[1, 2] <- facCor[2, 1] <- 0.4
facCor[1, 3] <- facCor[3, 1] <- 0.4
facCor[2, 3] <- facCor[3, 2] <- 0.4
totalVar <- rep(1, 9)
findFactorResidualVar(path, facCor, totalVar)

Description

Find factor total covariances from regression coefficient matrix, factor residual covariance matrix. The residual covaraince matrix might be derived from factor residual correlation, total variance, and error variance. This function can be applied for path analysis model as well.

Usage

findFactorTotalCov(beta, psi = NULL, corPsi = NULL, totalVarPsi = NULL, 
    errorVarPsi = NULL, gamma = NULL, covcov = NULL)

Arguments

Value

A matrix of factor (model-implied) total covariance

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndMean to find indicator (measurement) total means

• findIndResidualVar to find indicator (measurement) residual variances

• findIndTotalVar to find indicator (measurement) total variances

• findFactorIntercept to find factor intercepts

• findFactorMean to find factor means

• findFactorResidualVar to find factor residual variances

• findFactorTotalVar to find factor total variances

Examples

path <- matrix(0, 9, 9)
path[4, 1] <- path[7, 4] <- 0.6
path[5, 2] <- path[8, 5] <- 0.6
path[6, 3] <- path[9, 6] <- 0.6
path[5, 1] <- path[8, 4] <- 0.4
path[6, 2] <- path[9, 5] <- 0.4
facCor <- diag(9)
facCor[1, 2] <- facCor[2, 1] <- 0.4
facCor[1, 3] <- facCor[3, 1] <- 0.4
facCor[2, 3] <- facCor[3, 2] <- 0.4
residualVar <- c(1, 1, 1, 0.64, 0.288, 0.288, 0.64, 0.29568, 0.21888)
findFactorTotalCov(path, corPsi=facCor, errorVarPsi=residualVar)

Description

Find factor total variances from regression coefficient matrix, factor (residual) correlation matrix, and factor residual variances for latent variable models. In the path analysis model, this function will find indicator total variances from regression coefficient, indicator (residual) correlation matrix, and indicator residual variances.

Usage

findFactorTotalVar(beta, corPsi, residualVarPsi, gamma = NULL, covcov = NULL)

Arguments

Value

A vector of factor (indicator) total variances

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndMean to find indicator (measurement) total means

• findIndResidualVar to find indicator (measurement) residual variances

• findIndTotalVar to find indicator (measurement) total variances

• findFactorIntercept to find factor intercepts

• findFactorMean to find factor means

• findFactorResidualVar to find factor residual variances

• findFactorTotalCov to find factor covariances

Examples

path <- matrix(0, 9, 9)
path[4, 1] <- path[7, 4] <- 0.6
path[5, 2] <- path[8, 5] <- 0.6
path[6, 3] <- path[9, 6] <- 0.6
path[5, 1] <- path[8, 4] <- 0.4
path[6, 2] <- path[9, 5] <- 0.4
facCor <- diag(9)
facCor[1, 2] <- facCor[2, 1] <- 0.4
facCor[1, 3] <- facCor[3, 1] <- 0.4
facCor[2, 3] <- facCor[3, 2] <- 0.4
residualVar <- c(1, 1, 1, 0.64, 0.288, 0.288, 0.64, 0.29568, 0.21888)
findFactorTotalVar(path, facCor, residualVar)

Description

Find indicator (measurement) intercepts from a factor loading matrix, total factor mean, and indicator mean.

Usage

findIndIntercept(lambda, factorMean = NULL, indicatorMean = NULL, 
	kappa = NULL, covmean = NULL)

Arguments

Value

A vector of indicator (measurement) intercepts.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndMean to find indicator (measurement) total means

• findIndResidualVar to find indicator (measurement) residual variances

• findIndTotalVar to find indicator (measurement) total variances

• findFactorIntercept to find factor intercepts

• findFactorMean to find factor means

• findFactorResidualVar to find factor residual variances

• findFactorTotalVar to find factor total variances

• findFactorTotalCov to find factor covariances

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- c(0.6, 0.7, 0.8)
loading[4:6, 2] <- c(0.6, 0.7, 0.8)
facMean <- c(0.5, 0.2)
indMean <- rep(1, 6)
findIndIntercept(loading, facMean, indMean)

Description

Find indicator total means from a factor loading matrix, total factor means, and indicator (measurement) intercepts.

Usage

findIndMean(lambda, factorMean = NULL, tau = NULL, kappa = NULL, 
	covmean = NULL)

Arguments

Value

A vector of indicator total means.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndResidualVar to find indicator (measurement) residual variances

• findIndTotalVar to find indicator (measurement) total variances

• findFactorIntercept to find factor intercepts

• findFactorMean to find factor means

• findFactorResidualVar to find factor residual variances

• findFactorTotalVar to find factor total variances

• findFactorTotalCov to find factor covariances

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- c(0.6, 0.7, 0.8)
loading[4:6, 2] <- c(0.6, 0.7, 0.8)
facMean <- c(0.5, 0.2)
intcept <- rep(0, 6)
findIndMean(loading, facMean, intcept)

Description

Find indicator (measurement) residual variances from a factor loading matrix, total factor covariance matrix, and total indicator variances.

Usage

findIndResidualVar(lambda, totalFactorCov, totalVarTheta = NULL, 
	kappa = NULL, covcov = NULL)

Arguments

Value

A vector of indicator residual variances.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndMean to find indicator (measurement) total means

• findIndTotalVar to find indicator (measurement) total variances

• findFactorIntercept to find factor intercepts

• findFactorMean to find factor means

• findFactorResidualVar to find factor residual variances

• findFactorTotalVar to find factor total variances

• findFactorTotalCov to find factor covariances

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- c(0.6, 0.7, 0.8)
loading[4:6, 2] <- c(0.6, 0.7, 0.8)
facCov <- matrix(c(1, 0.5, 0.5, 1), 2, 2)
totalVar <- rep(1, 6)
findIndResidualVar(loading, facCov, totalVar)

Description

Find indicator total variances from a factor loading matrix, total factor covariance matrix, and indicator (measurement) residual variances.

Usage

findIndTotalVar(lambda, totalFactorCov, residualVarTheta, kappa = NULL, 
	covcov = NULL)

Arguments

Value

A vector of indicator total variances.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findIndIntercept to find indicator (measurement) intercepts

• findIndMean to find indicator (measurement) total means

• findIndResidualVar to find indicator (measurement) residual variances

• findFactorIntercept to find factor intercepts

• findFactorMean to find factor means

• findFactorResidualVar to find factor residual variances

• findFactorTotalVar to find factor total variances

• findFactorTotalCov to find factor covariances

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- c(0.6, 0.7, 0.8)
loading[4:6, 2] <- c(0.6, 0.7, 0.8)
facCov <- matrix(c(1, 0.5, 0.5, 1), 2, 2)
resVar <- c(0.64, 0.51, 0.36, 0.64, 0.51, 0.36)
findIndTotalVar(loading, facCov, resVar)

Description

Find the appropriate position for freely estimated correlation (or covariance) given a regression coefficient matrix. The appropriate position is the pair of variables that are not causally related.

Usage

findPossibleFactorCor(beta)

Arguments

Value

The symmetric matrix containing the appropriate position for freely estimated correlation.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findRecursiveSet to group variables regarding the position in mediation chain.

Examples

path <- matrix(0, 9, 9)
path[4, 1] <- path[7, 4] <- NA
path[5, 2] <- path[8, 5] <- NA
path[6, 3] <- path[9, 6] <- NA
path[5, 1] <- path[8, 4] <- NA
path[6, 2] <- path[9, 5] <- NA
findPossibleFactorCor(path)

Description

Find a value of independent variable that provides a given value of power. If there are more than one varying parameters, this function will find the value of the target varying parameters given the values of the other varying parameters.

Usage

findPower(powerTable, iv, power)

Arguments

Value

There are five possible types of values provided:

• Value The varying parameter value that provides the power just over the specified power value (the adjacent value of varying parameter provides lower power than the specified power value).

• Minimum value The minimum value has already provided enough power (way over the specified power value). The value of varying parameters that provides exact desired power may be lower than the minimum value. The example of varying parameter that can provides the minimum value is sample size.

• Maximum value The maximum value has already provided enough power (way over the specified power value). The value of varying parameters that provides exact desired power may be higher than the maximum value. The example of varying parameter that can provides the maximum value is percent missing.

• NA There is no value in the domain of varying parameters that provides the power greater than the desired power.

• Inf The power of all values in the varying parameters is 1 (specifically more than 0.9999) and any values of the varying parameters can be picked and still provide enough power.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• getPower to find the power of parameter estimates

• continuousPower to find the power of parameter estimates for the result object (linkS4class{SimResult}) with varying parameters.

Examples

## Not run: 
# Specify Sample Size by n
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, 0.4)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# Specify both sample size and percent missing completely at random. Note that more fine-grained 
# values of n and pmMCAR is needed, e.g., n=seq(50, 500, 1) and pmMCAR=seq(0, 0.2, 0.01)
Output <- sim(NULL, model=CFA.Model, n=seq(100, 200, 20), pmMCAR=c(0, 0.1, 0.2))

# Find the power of all possible combination of N and pmMCAR
pow <- getPower(Output)

# Find the sample size that provides the power of 0.8
findPower(pow, "N", 0.80)

## End(Not run)

Description

In mediation analysis, variables affects other variables as a chain. This function will group variables regarding the chain of mediation analysis.

Usage

findRecursiveSet(beta)

Arguments

Value

The list of set of variables in the mediation chain. The variables in position 1 will be the independent variables. The variables in the last variables will be the end of the chain.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• findPossibleFactorCor to find the possible position for latent correlation given a regression coefficient matrix

Examples

path <- matrix(0, 9, 9)
path[4, 1] <- path[7, 4] <- NA
path[5, 2] <- path[8, 5] <- NA
path[6, 3] <- path[9, 6] <- NA
path[5, 1] <- path[8, 4] <- NA
path[6, 2] <- path[9, 5] <- NA
findRecursiveSet(path)

Description

This function can be used to generate random data based on the 1. SimSem objects created with the model function, 2. lavaan script or parameter tables, or 3. an MxModel object from the OpenMx package. Some notable features include fine control of misspecification and misspecification optimization (for SimSem only), as well as the ability to generate non-normal data. When using simsem for simulations, this function is used internally to generate data in the function sim, and can be helpful for debugging, or in creating data for use with other analysis programs.

Usage

generate(model, n, maxDraw=50, misfitBounds=NULL, misfitType="f0",
	averageNumMisspec=FALSE, optMisfit=NULL, optDraws=50, 
	createOrder = c(1, 2, 3), indDist=NULL, sequential=FALSE,	
	facDist=NULL, errorDist=NULL, saveLatentVar = FALSE, indLab=NULL, 
	modelBoot=FALSE, realData=NULL, covData=NULL, params=FALSE, group = NULL, 
	empirical = FALSE, ...)

Arguments

Details

If the lavaan script or the MxModel are provided, the model-implied covariance matrix will be computed and internally use createData function to generate data. The data-generation method is based on whether the indDist argument is specified. For the lavaan script, the code for data generation is modified from the simulateData function.

If the SimSem object is specified, it will check whether there are any random parameters or trivial misspecification in the model. If so, real or misspecified parameters are drawn via the draw function. Next, there are two methods to generate data. First, the function will calculate the model-implied covariance matrix (including model misspecification) and generate data similar to the lavaan script or the MxModel object. The second method is referred to as the sequential method, which can be used by specifying the sequential argument as TRUE. This function will create data based on the chain of equations in structural equation modeling such that independent variables and errors are generated and added as dependent variables and the dependent variables will be treated as independent variables in the next equation. For example, in the model with factor A and B are independent variables, factor C are dependent variables, factors A and B are generated first. Then, residual in factor C are created and added with factors A and B. This current step has all factor scores. Then, measurement errors are created and added with factor scores to create indicator scores. During each step, independent variables and errors can be nonnormal by setting facDist or errorDist arguments. The data generation in each step is based on the createData function.

For the model-based bootstrap (providing the realData argument), the transformation proposed by Yung & Bentler (1996) is used. This procedure is the expansion from the Bollen and Stine (1992) bootstrap including a mean structure. The model-implied mean vector and covariance matrix with trivial misspecification will be used in the model-based bootstrap if misspec is specified. See page 133 of Bollen and Stine (1992) for a reference.

Value

A data.frame containing simulated data from the data generation template. A variable "group" is appended indicating group membership.

Author(s)

Sunthud Pornprasertmanit ([email protected]), Patrick Miller (University of Notre Dame; [email protected]), the data generation code for lavaan script is modifed from the simulateData function in lavaan written by Yves Rosseel

References

Bollen, K. A., & Stine, R. A. (1992). Bootstrapping goodness-of-fit measures in structural equation models. Sociological Methods and Research, 21, 205-229.

Yung, Y.-F., & Bentler, P. M. (1996). Bootstrapping techniques in analysis of mean and covariance structures. In G. A. Marcoulides & R. E. Schumacker (Eds.), Advanced structural equation modeling: Issues and techniques (pp. 195-226). Mahwah, NJ: Erlbaum.

See Also

• draw To draw parameters using the SimSem template.

• createData To generate random data using a set of parameters from draw

Examples

loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
LY <- bind(loading, 0.7)

latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)

RTE <- binds(diag(6))

VY <- bind(rep(NA,6),2)

CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType = "CFA")

dat <- generate(CFA.Model, 200)

# Get the latent variable scores

dat2 <- generate(CFA.Model, 20, sequential = TRUE, saveLatentVar = TRUE)
dat2
attr(dat2, "latentVar")

Description

Find the median of confidence interval width or a confidence interval value given a degree of assurance (Lai & Kelley, 2011)

Usage

getCIwidth(object, assurance = 0.50, nVal = NULL, pmMCARval = NULL, 
	pmMARval = NULL, df = 0)

Arguments

Value

The median of confidence interval width or a confidence interval given a degree of assurance

Author(s)

Sunthud Pornprasertmanit ([email protected])

References

Lai, K., & Kelley, K. (2011). Accuracy in parameter estimation for targeted effects in structural equation modeling: Sample size planning for narrow confidence intervals. Psychological Methods, 16, 127-148.

See Also

SimResult for a detail of simResult

Examples

## Not run: 
loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
loadingValues <- matrix(0, 6, 2)
loadingValues[1:3, 1] <- 0.7
loadingValues[4:6, 2] <- 0.7
LY <- bind(loading, loadingValues)
latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)
error.cor <- matrix(0, 6, 6)
diag(error.cor) <- 1
RTE <- binds(error.cor)
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# We make the examples running only 5 replications to save time.
# In reality, more replications are needed.
Output <- sim(5, n = 200, model=CFA.Model)

# Get the cutoff (critical value) when alpha is 0.05
getCIwidth(Output, assurance=0.80)

# Finding the cutoff when the sample size is varied. Note that more fine-grained 
# values of n is needed, e.g., n=seq(50, 500, 1)
Output2 <- sim(NULL, model=CFA.Model, n=seq(50, 100, 10))

# Get the fit index cutoff when sample size is 75.
getCIwidth(Output2, assurance=0.80, nVal = 75)

## End(Not run)

Description

A function to find the coverage rate of confidence intervals in a model when none, one, or more of the simulations parameters vary randomly across replications.

Usage

getCoverage(simResult, coverValue = NULL, contParam = NULL, coverParam = NULL, 
    nVal = NULL, pmMCARval = NULL, pmMARval = NULL, paramVal = NULL)

Arguments

Details

In this function, the coverage (which can be 0 or 1) is regressed on randomly varying simulation parameters (e.g., sample size, percentage of missing data, or model parameters) using logistic regression. For a set of independent variables values, the predicted probability from the logistic regression equation is the predicted coverage rate.

Value

Data frame containing columns representing values of the randomly varying simulation parameters, and coverage rates for model parameters of interest.

Author(s)

Sunthud Pornprasertmanit ([email protected]), Alexander M. Schoemann (East Carolina University; [email protected])

See Also

• SimResult to see how to create a simResult object with randomly varying parameters.

Examples

## Not run: 
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, 0.7)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# Specify both sample size and percent missing completely at random. Note that more fine-grained 
# values of n and pmMCAR is needed, e.g., n=seq(50, 500, 1) and pmMCAR=seq(0, 0.2, 0.01)
Output <- sim(NULL, model=CFA.Model, n=seq(100, 200, 20), pmMCAR=c(0, 0.1, 0.2))
summary(Output)

# Get the coverage rates of all possible combinations of n and pmMCAR
getCoverage(Output)

# Get the coverage rates of the combinations of n of 100 and 200 and pmMCAR of 0, 0.1, and 0.2
getCoverage(Output, coverValue = 0, nVal=c(100, 200), pmMCARval=c(0, 0.1, 0.2))

## End(Not run)

Description

Extract fit indices information from the SimResult and get the cutoffs of fit indices given a priori alpha level

Usage

getCutoff(object, alpha, revDirec = FALSE, usedFit = NULL, nVal = NULL, 
	pmMCARval = NULL, pmMARval = NULL, df = 0)

Arguments

Value

One-tailed cutoffs of several fit indices with a priori alpha level

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

SimResult for a detail of simResult

Examples

## Not run: 
loading <- matrix(0, 6, 2)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
loadingValues <- matrix(0, 6, 2)
loadingValues[1:3, 1] <- 0.7
loadingValues[4:6, 2] <- 0.7
LY <- bind(loading, loadingValues)
latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, 0.5)
error.cor <- matrix(0, 6, 6)
diag(error.cor) <- 1
RTE <- binds(error.cor)
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# We make the examples running only 5 replications to save time.
# In reality, more replications are needed.
Output <- sim(5, n = 200, model=CFA.Model)

# Get the cutoff (critical value) when alpha is 0.05
getCutoff(Output, 0.05)

# Finding the cutoff when the sample size is varied. Note that more fine-grained 
# values of n is needed, e.g., n=seq(50, 500, 1)
Output2 <- sim(NULL, model=CFA.Model, n=seq(50, 100, 10))

# Get the fit index cutoff when sample size is 75.
getCutoff(Output2, 0.05, nVal = 75)

## End(Not run)

Description

Extract fit indices information from the simulation of parent and nested models and getCutoff of fit indices given a priori alpha level

Usage

getCutoffNested(nested, parent, alpha = 0.05, usedFit = NULL, nVal = NULL, 
	pmMCARval = NULL, pmMARval = NULL, df = 0)

Arguments

Value

One-tailed cutoffs of several fit indices with a priori alpha level

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

SimResult for a detail of simResult getCutoff for a detail of finding cutoffs for absolute fit

Examples

## Not run: 
# Nested Model
loading.null <- matrix(0, 6, 1)
loading.null[1:6, 1] <- NA
LY.NULL <- bind(loading.null, 0.7)
RPS.NULL <- binds(diag(1))

error.cor.mis <- matrix("rnorm(1, 0, 0.1)", 6, 6)
diag(error.cor.mis) <- 1
RTE <- binds(diag(6), misspec=error.cor.mis)
CFA.Model.NULL <- model(LY = LY.NULL, RPS = RPS.NULL, RTE = RTE, modelType="CFA")

# Parent Model
loading.alt <- matrix(0, 6, 2)
loading.alt[1:3, 1] <- NA
loading.alt[4:6, 2] <- NA
LY.ALT <- bind(loading.alt, 0.7)
latent.cor.alt <- matrix(NA, 2, 2)
diag(latent.cor.alt) <- 1
RPS.ALT <- binds(latent.cor.alt, "runif(1, 0.7, 0.9)")
CFA.Model.ALT <- model(LY = LY.ALT, RPS = RPS.ALT, RTE = RTE, modelType="CFA")

# The actual number of replications should be greater than 10.
Output.NULL.NULL <- sim(10, n=500, model=CFA.Model.NULL, generate=CFA.Model.NULL)
Output.NULL.ALT <- sim(10, n=500, model=CFA.Model.ALT, generate=CFA.Model.NULL)

# Find the fix index cutoff from the sampling distribution of the difference
# in fit index of nested models where the alpha is 0.05.
getCutoffNested(Output.NULL.NULL, Output.NULL.ALT, alpha=0.05)

## End(Not run)

Description

Extract fit indices information from the simulation of two models fitting on the datasets created from both models and getCutoff of fit indices given a priori alpha level

Usage

getCutoffNonNested(dat1Mod1, dat1Mod2, dat2Mod1=NULL, dat2Mod2=NULL, 
alpha=.05, usedFit=NULL, onetailed=FALSE, nVal = NULL, pmMCARval = NULL, 
pmMARval = NULL, df = 0)

Arguments

Value

One- or two-tailed cutoffs of several fit indices with a priori alpha level. The cutoff is based on the fit indices from Model 1 subtracted by the fit indices from Model 2.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

SimResult for a detail of simResult getCutoff for a detail of finding cutoffs for absolute fit getCutoffNested for a detail of finding cutoffs for nested model comparison plotCutoffNonNested Plot cutoffs for non-nested model comparison

Examples

## Not run: 
# Model A: Factor 1 with items 1-3 and Factor 2 with items 4-8
loading.A <- matrix(0, 8, 2)
loading.A[1:3, 1] <- NA
loading.A[4:8, 2] <- NA
LY.A <- bind(loading.A, 0.7)
latent.cor <- matrix(NA, 2, 2)
diag(latent.cor) <- 1
RPS <- binds(latent.cor, "runif(1, 0.7, 0.9)")
RTE <- binds(diag(8))
CFA.Model.A <- model(LY = LY.A, RPS = RPS, RTE = RTE, modelType="CFA")

# Model B: Factor 1 with items 1-4 and Factor 2 with items 5-8
loading.B <- matrix(0, 8, 2)
loading.B[1:4, 1] <- NA
loading.B[5:8, 2] <- NA
LY.B <- bind(loading.B, 0.7)
CFA.Model.B <- model(LY = LY.B, RPS = RPS, RTE = RTE, modelType="CFA")

# The actual number of replications should be greater than 10.
Output.A.A <- sim(10, n=500, model=CFA.Model.A, generate=CFA.Model.A)
Output.A.B <- sim(10, n=500, model=CFA.Model.B, generate=CFA.Model.A)
Output.B.A <- sim(10, n=500, model=CFA.Model.A, generate=CFA.Model.B)
Output.B.B <- sim(10, n=500, model=CFA.Model.B, generate=CFA.Model.B)

# Find the cutoffs from the sampling distribution to reject model A (model 1)
# and to reject model B (model 2)
getCutoffNonNested(Output.A.A, Output.A.B, Output.B.A, Output.B.B)

# Find the cutoffs from the sampling distribution to reject model A (model 1)
getCutoffNonNested(Output.A.A, Output.A.B)

# Find the cutoffs from the sampling distribution to reject model B (model 1)
getCutoffNonNested(Output.B.B, Output.B.A)

## End(Not run)

Description

Get extra outputs from a simulation result object (SimResult). Users can ask this package to extra output from the lavaan object in each iteration by setting the outfun argument (in the sim function). See the example below.

Usage

getExtraOutput(object, improper = TRUE, nonconverged = FALSE)

Arguments

Value

A list of extra outputs

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• sim A function to run a Monte Carlo simulation

Examples

## Not run: 
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, 0.7)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# Write a function to extract the modification index from lavaan object
outfun <- function(out) {
	result <- inspect(out, "mi")
}

# We will use only 5 replications to save time.
# In reality, more replications are needed.
Output <- sim(5, n=200, model=CFA.Model, outfun=outfun)

# Get the modification index of each replication
getExtraOutput(Output)

## End(Not run)

Description

This function will extract the data generation population model underlying a result object (linkS4class{SimResult}).

Usage

getPopulation(object, std = FALSE, improper = TRUE, nonconverged = FALSE)

Arguments

Value

A data frame contained the population of each replication

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• SimResult for result object

Examples

## Not run: 
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, "runif(1, 0.4, 0.9)")
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# We will use only 10 replications to save time.
# In reality, more replications are needed.
Output <- sim(10, n=200, model=CFA.Model)

# Get the population parameters
getPopulation(Output)

## End(Not run)

Description

A function to find the power of parameters in a model when none, one, or more of the simulations parameters vary randomly across replications.

Usage

getPower(simResult, alpha = 0.05, contParam = NULL, powerParam = NULL, 
	nVal = NULL, pmMCARval = NULL, pmMARval = NULL, paramVal = NULL)

Arguments

Details

A common use of simulations is to conduct power analyses, especially when using SEM (Muthen & Muthen, 2002). Here, researchers could select values for each parameter and a sample size and run a simulation to determine power in those conditions (the proportion of generated datasets in which a particular parameter of interest is significantly different from zero). To evaluate power at multiple sample sizes, one simulation for each sample size must be run. This function not only calculate power for each sample size but also calculate power for multiple sample sizes varying continuously. By continuously varying sample size across replications, only a single simulation is needed. In this simulation, the sample size for each replication varies randomly across plausible sample sizes (e.g., sample sizes between 200 and 500). For each replication, the sample size and significance of each parameter (0 = not significant, 1 = significant) are recorded. When the simulation is complete, parameter significance is regressed on sample size using logistic regression. For a given sample size, the predicted probability from the logistic regression equation is the power to detect an effect at that sample size. This approach can be extended to other randomly varying simulation parameters such as the percentage of missing data, and model parameters.

Value

Data frame containing columns representing values of the randomly varying simulation parameters, and power for model parameters of interest.

Author(s)

Alexander M. Schoemann (East Carolina University; [email protected]), Sunthud Pornprasertmanit ([email protected])

References

Muthen, L. K., & Muthen, B. O. (2002). How to use a Monte Carlo study to decide on sample size and determine power. Structural Equation Modeling, 4, 599-620.

See Also

• SimResult to see how to create a simResult object with randomly varying parameters.

Examples

## Not run: 
loading <- matrix(0, 6, 1)
loading[1:6, 1] <- NA
LY <- bind(loading, 0.7)
RPS <- binds(diag(1))
RTE <- binds(diag(6))
CFA.Model <- model(LY = LY, RPS = RPS, RTE = RTE, modelType="CFA")

# Specify both sample size and percent missing completely at random. Note that more fine-grained 
# values of n and pmMCAR is needed, e.g., n=seq(50, 500, 1) and pmMCAR=seq(0, 0.2, 0.01)
Output <- sim(NULL, model=CFA.Model, n=seq(100, 200, 20), pmMCAR=c(0, 0.1, 0.2))
summary(Output)

# Get the power of all possible combinations of n and pmMCAR
getPower(Output)

# Get the power of the combinations of n of 100 and 200 and pmMCAR of 0, 0.1, and 0.2
getPower(Output, nVal=c(100, 200), pmMCARval=c(0, 0.1, 0.2))

## End(Not run)

Description

Find the proportion of fit indices that indicate worse fit than a specified cutoffs. The cutoffs may be calculated from getCutoff of the null model.

Usage

getPowerFit(altObject, cutoff = NULL, nullObject = NULL, revDirec = FALSE, 
usedFit = NULL, alpha = 0.05, nVal = NULL, pmMCARval = NULL, pmMARval = NULL, 
condCutoff = TRUE, df = 0)

Arguments

Value

List of power given different fit indices. The TraditionalChi means the proportion of replications that are rejected by the traditional chi-square test.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

• getCutoff to find the cutoffs from null model.

• SimResult to see how to create simResult