Structural Equation Modeling: Lecture 15 Sample Size Planning
Sunthud Pornprasertmanit
28/10/2025
Power Analysis for Parameter Esimates
In structural equation modeling (SEM), one of the most popular ways to plan sample size is through simulation-based power analysis. This method involves repeatedly generating data under a specified model and then analyzing each dataset to check the proportion of replications in which the desired parameters are statistically significant.
In this tutorial, we will use the simsem package to
conduct a simulation-based power analysis. The general steps are as
follows:
- Write a model for data generation. This model represents the population model with known population parameters.
- Write a model for data analysis. This is the model you intend to fit to the data (often the same as the generating model).
- Use the
sim()function, specifying these arguments:
nRep: Number of replications (e.g., 1000).generate: The data-generating model (the population model).model: The analysis model (the one being fitted).n: Sample size per replication.lavaanfun: The lavaan function used in the analysis (e.g.,"cfa"or"sem"). If not specified,simsemuses the genericlavaan()function by default, which is suitable unless you need special defaults fromcfa()orsem().
- Summarize and interpret the results. Use
summary()to inspect relevant output, such as statistical power for specific parameters.
Example 1
In this example, we will perform a power analysis for a confirmatory factor analysis (CFA) model with two latent factors—one measured by three indicators, and the other measured by two indicators.
The model used for data generation is written in lavaan
syntax, with parameter values specified before each parameter. In
lavaan, these parameter values are treated as
unstandardized, so we must set appropriate error variances to ensure
that both the latent variables and indicators have variances of 1.
Let’s start with the data-generating model:
library(simsem)
genmodel <- "
f1 =~ 0.7*x1 + 0.7*x2 + 0.7*x3
f2 =~ 0.8*x4 + 0.9*x5
f2 ~~ 0.3*f1
f1 ~~ 1*f1
f2 ~~ 1*f2
x1 ~~ 0.51*x1
x2 ~~ 0.51*x2
x3 ~~ 0.51*x3
x4 ~~ 0.36*x4
x5 ~~ 0.19*x5
"Let’s go through this model line by line:
f1 =~ 0.7*x1 + 0.7*x2 + 0.7*x3: Factor 1 loads on Indicators 1–3 with loadings of 0.7 each.f2 =~ 0.8*x4 + 0.9*x5: Factor 2 loads on Indicators 4–5 with loadings of 0.8 and 0.9, respectively.f2 ~~ 0.3*f1: The covariance between Factor 1 and Factor 2 is 0.3.f1 ~~ 1*f1andf2 ~~ 1*f2: Each factor has a variance of 1. Consequently, the covariance of 0.3 represents a factor correlation of 0.3.x1 ~~ 0.51*x1: The residual variance of Indicator 1 is 0.51, calculated so that its total variance equals 1. Since each indicator is standardized, \(1- \lambda^2 = 1 - 0.7^2 = 0.51\).x2 ~~ 0.51*x2andx3 ~~ 0.51*x3: The same logic applies to Indicators 2 and 3.x4 ~~ 0.36*x4: Based on \(1 - 0.8^2 = 0.36\), the residual variance of Indicator 4 is 0.36.x5 ~~ 0.19*x5: From \(1 - 0.9^2 = 0.19\), the residual variance of Indicator 5 is 0.19.
Once all relationships and variances have been specified, we can generate data using either:
simulateData()fromlavaanorgenerate()fromsimsem.
Below, we use the generate() function to draw a sample
size of 200.
set.seed(123321)
dattemp <- generate(genmodel, n=200)
head(dattemp)## x1 x2 x3 x4 x5
## 1 -0.1139278 -0.86757102 -0.53516500 -0.7590253 -0.4473281
## 2 0.1277854 -0.19721568 0.36620753 -2.0315035 -1.7736733
## 3 -0.2606923 -1.64988424 0.09594833 0.1829244 -0.1324875
## 4 -0.6337764 -0.52339101 -2.07926927 -0.3039608 -0.7709616
## 5 0.8230339 -0.00208385 2.35333017 0.7151418 0.4034206
## 6 0.2876113 -0.62497258 -1.84986097 -0.1496067 1.1965293We set a random seed (set.seed()) to ensure
reproducibility — that is, the same random sample is generated each time
the code is run.
Next, we define the analysis model, which is the model we will fit to the simulated data:
analysismodel <- "
f1 =~ x1 + x2 + x3
f2 =~ x4 + x5
"Now let’s fit this model to the generated dataset:
fittemp <- cfa(analysismodel, data=dattemp, std.lv=TRUE)
fittemp## lavaan 0.6-20 ended normally after 17 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 11
##
## Number of observations 200
##
## Model Test User Model:
##
## Test statistic 1.408
## Degrees of freedom 4
## P-value (Chi-square) 0.843Running the Simulation
Now we use the sim() function to conduct the power
analysis. In this case, we will run 1,000 replications, each using a
sample size of 200. The argument lavaanfun = "cfa" ensures
that the same defaults as the cfa() function are applied.
If additional arguments (e.g., std.lv = TRUE) are used in
the analysis model, simply include them inside the sim()
function as well.
simcfa200 <- sim(nRep=1000, model=analysismodel, generate=genmodel, n=200, lavaanfun="cfa", std.lv=TRUE)You’ll see the simulation progress in the R console (the progress output is suppressed here for readability).
Once the simulation is complete, we can summarize the results:
summary(simcfa200)## RESULT OBJECT
## Model Type
## [1] "lavaan"
## ========= Fit Indices Cutoffs ============
## Alpha
## Fit Indices 0.1 0.05 0.01 0.001 Mean SD
## chisq 8.155 10.109 13.627 23.388 4.181 3.022
## aic 2609.487 2627.570 2658.012 2693.257 2553.810 45.201
## bic 2645.768 2663.851 2694.293 2729.539 2590.092 45.201
## rmsea 0.072 0.087 0.110 0.156 0.022 0.032
## cfi 0.985 0.979 0.966 0.944 0.996 0.008
## tli 0.964 0.948 0.914 0.861 0.999 0.026
## srmr 0.037 0.041 0.052 0.064 0.023 0.010
## ========= Parameter Estimates and Standard Errors ============
## Estimate Average Estimate SD Average SE Power (Not equal 0) Std Est
## f1=~x1 0.697 0.073 0.075 1.000 0.700
## f1=~x2 0.691 0.077 0.075 1.000 0.693
## f1=~x3 0.699 0.079 0.075 1.000 0.701
## f2=~x4 0.804 0.149 0.155 0.978 0.807
## f2=~x5 0.920 0.279 0.394 0.979 0.923
## x1~~x1 0.502 0.079 0.079 1.000 0.507
## x2~~x2 0.509 0.081 0.079 1.000 0.516
## x3~~x3 0.499 0.081 0.079 1.000 0.505
## x4~~x4 0.325 0.254 0.262 0.584 0.329
## x5~~x5 0.073 1.793 3.548 0.261 0.065
## f1~~f2 0.293 0.088 0.085 0.906 0.293
## Std Est SD Std Ave SE Average Param Average Bias Coverage
## f1=~x1 0.056 0.057 0.70 -0.003 0.949
## f1=~x2 0.059 0.057 0.70 -0.009 0.942
## f1=~x3 0.060 0.057 0.70 -0.001 0.938
## f2=~x4 0.141 0.144 0.80 0.004 0.967
## f2=~x5 0.289 0.401 0.90 0.020 0.962
## x1~~x1 0.078 0.079 0.51 -0.008 0.954
## x2~~x2 0.081 0.078 0.51 -0.001 0.950
## x3~~x3 0.083 0.079 0.51 -0.011 0.947
## x4~~x4 0.258 0.266 0.36 -0.035 0.971
## x5~~x5 2.049 4.036 0.19 -0.117 0.954
## f1~~f2 0.088 0.085 0.30 -0.007 0.955
## ========= Correlation between Fit Indices ============
## chisq aic bic rmsea cfi tli srmr
## chisq 1.000 0.016 0.016 0.948 -0.935 -0.990 0.818
## aic 0.016 1.000 1.000 0.018 -0.033 -0.014 -0.009
## bic 0.016 1.000 1.000 0.018 -0.033 -0.014 -0.009
## rmsea 0.948 0.018 0.018 1.000 -0.935 -0.945 0.761
## cfi -0.935 -0.033 -0.033 -0.935 1.000 0.936 -0.723
## tli -0.990 -0.014 -0.014 -0.945 0.936 1.000 -0.822
## srmr 0.818 -0.009 -0.009 0.761 -0.723 -0.822 1.000
## ================== Replications =====================
## Number of replications = 1000
## Number of converged replications = 746
## Number of nonconverged replications:
## 1. Nonconvergent Results = 1
## 2. Nonconvergent results from multiple imputation = 0
## 3. At least one SE were negative or NA = 0
## 4. Nonpositive-definite latent or observed (residual) covariance matrix
## (e.g., Heywood case or linear dependency) = 253Understanding the Output
The output of summary(simcfa200) consists of several
sections:
- Model Type Indicates that the model was estimated using lavaan.
- Fit Indices Cutoffs Since this is a correctly specified model, this section displays Type I error cutoffs (i.e., how often a correct model is incorrectly rejected). Note that with only five indicators, the RMSEA cutoff for \(\alpha = .05\) may appear quite large (e.g., .11). This demonstrates why you should not rely solely on universal cutoff values such as those from Hu & Bentler (1999).
- Parameter Estimates and Standard Errors A detailed table summarizing results across replications, including:
- Estimate Average: The mean parameter estimate.
- Estimate SD: The standard deviation of parameter estimates.
- Average SE: The mean estimated standard error.
- Power (Not equal 0): The proportion of replications in which the parameter was statistically significant.
- Std Est: The mean standardized estimate.
- Std Est SD: The standard deviation of standardized estimates.
- Std Ave SE: The mean standard error of standardized estimates.
- Average Param: The parameter value from the generating model.
- Average Bias: Difference between Estimate Average and Average Param.
- Coverage: The proportion of replications in which the 95% confidence interval contained the true parameter.
- Correlation between Fit Indices Displays how fit indices co-vary across replications. For instance, RMSEA, CFI, and TLI are typically highly correlated (around .93–.95).
- Replications Summary Reports the number of converged replications. In this example, 746 out of 1000 models converged successfully. Most non-convergent replications were due to non–positive definite matrices (e.g., negative error variances), which is common when using small sample sizes such as 200 (especially for factors with two indicators). This means that even a “perfect” model can fail to converge about 25% of the time when the sample size is relatively small.
Finally, let’s focus on the factor correlation between
f1 and f2. The simulation indicates that the
power to detect this relationship as significant is .906 when \(n = 200\). This level of power is generally
considered more than adequate.
Extracting Simulation Results
In addition to the overall summary, users may wish to extract
specific information from the simulation results. The
summaryParam() function provides a detailed table
summarizing the parameter estimates, standard errors, and power for each
parameter across all replications.
It is important to specify improper = FALSE to exclude
any replications with improper parameter estimates, such as negative
residual variances or non–positive definite covariance matrices. These
improper solutions are typically considered invalid and should be
excluded from interpretation.
summaryParam(simcfa200, improper=FALSE)## Estimate Average Estimate SD Average SE Power (Not equal 0) Std Est
## f1=~x1 0.6991461 0.07222300 0.07492922 1.0000000 0.7008212
## f1=~x2 0.6933521 0.07649999 0.07482418 1.0000000 0.6952343
## f1=~x3 0.6986571 0.07936156 0.07487691 1.0000000 0.7001101
## f2=~x4 0.8213794 0.08731213 0.11768363 0.9986595 0.8254259
## f2=~x5 0.8731900 0.09073817 0.12259915 1.0000000 0.8753623
## x1~~x1 0.5014231 0.07828823 0.07903144 1.0000000 0.5057870
## x2~~x2 0.5081150 0.08244929 0.07879882 1.0000000 0.5132124
## x3~~x3 0.5010003 0.08032269 0.07899593 1.0000000 0.5063105
## x4~~x4 0.3090155 0.11698684 0.17155278 0.5321716 0.3135523
## x5~~x5 0.2254477 0.12326375 0.19204563 0.2828418 0.2284716
## f1~~f2 0.3126772 0.08031484 0.08287795 0.9664879 0.3126772
## Std Est SD Std Ave SE Average Param Average Bias Coverage
## f1=~x1 0.05537823 0.05645460 0.70 -0.0008538519 0.9530831
## f1=~x2 0.05866360 0.05655419 0.70 -0.0066478711 0.9477212
## f1=~x3 0.05949854 0.05644957 0.70 -0.0013428849 0.9369973
## f2=~x4 0.07160065 0.10568905 0.80 0.0213793609 0.9839142
## f2=~x5 0.07263794 0.11082188 0.90 -0.0268099729 0.9758713
## x1~~x1 0.07684298 0.07872629 0.51 -0.0085769060 0.9504021
## x2~~x2 0.08114968 0.07817908 0.51 -0.0018849852 0.9490617
## x3~~x3 0.08266550 0.07856385 0.51 -0.0089996981 0.9504021
## x4~~x4 0.11952011 0.17510244 0.36 -0.0509844807 0.9959786
## x5~~x5 0.12600061 0.19492326 0.19 0.0354477174 0.9758713
## f1~~f2 0.08031484 0.08287795 0.30 0.0126771752 0.9571046Users can also examine the replication-level summary using the
summaryConverge() function. This function reports the
number of converged and nonconverged replications, along with the
reasons for nonconvergence. Common causes include negative variances,
non–positive definite latent covariance matrices, or numerical
underflow.
summaryConverge(simcfa200)## $Converged
## num.converged num.nonconverged
## 746 254
##
## $`Nonconvergent Reasons`
## count
## Nonconvergent 1
## Improper SE 0
## Nonpositive definite matrix 253
## Optimal estimates were not guaranteed 0Power at Different Sample Sizes
Let us now focus specifically on the factor correlation between
f1 and f2 and examine how its power changes
with different sample sizes.
We begin by reducing the sample size to 120 and running the simulation again.
simcfa120 <- sim(nRep=1000, model=analysismodel, generate=genmodel, n=120, lavaanfun="cfa", std.lv=TRUE)After completing the simulation, we can extract the power for the
factor correlation parameter using the summaryParam()
function:
summaryParam(simcfa120, improper=FALSE)["f1~~f2", "Power (Not equal 0)"]## [1] 0.8654781The result indicates that the power of the factor correlation remains relatively high, at approximately .86, even with a sample size of 120. This suggests that, in terms of statistical significance alone, 120 participants might be sufficient to detect the factor correlation.
However, statistical power is not the only concern in SEM. A smaller sample size can lead to a higher rate of nonconvergent solutions—replications in which the estimation process fails to produce admissible parameter estimates.
summaryConverge(simcfa120)## $Converged
## num.converged num.nonconverged
## 617 383
##
## $`Nonconvergent Reasons`
## count
## Nonconvergent 5
## Improper SE 0
## Nonpositive definite matrix 378
## Optimal estimates were not guaranteed 0At n = 120, the nonconvergence rate is approximately 38%, which is unacceptably high. Such a high failure rate indicates that, even though the power is adequate, the estimation procedure frequently encounters problematic solutions. In practice, this level of instability would make the analysis unreliable.
Increasing the Sample Size
To address this issue, let’s increase the sample size to 400 and observe the convergence rate.
simcfa400 <- sim(nRep=1000, model=analysismodel, generate=genmodel, n=400, lavaanfun="cfa", std.lv=TRUE)summaryConverge(simcfa400)## $Converged
## num.converged num.nonconverged
## 861 139
##
## $`Nonconvergent Reasons`
## count
## Nonconvergent 0
## Improper SE 0
## Nonpositive definite matrix 139
## Optimal estimates were not guaranteed 0With n = 400, the convergence rate improves substantially to 86%, suggesting that a majority of replications yield admissible solutions.
Next, let’s further increase the sample size to 800 to check whether convergence continues to improve.
simcfa800 <- sim(nRep=1000, model=analysismodel, generate=genmodel, n=800, lavaanfun="cfa", std.lv=TRUE)summaryConverge(simcfa800)## $Converged
## num.converged num.nonconverged
## 951 49
##
## $`Nonconvergent Reasons`
## count
## Nonconvergent 0
## Improper SE 0
## Nonpositive definite matrix 49
## Optimal estimates were not guaranteed 0At n = 800, the convergence rate rises to approximately 95%, indicating that nearly all replications now converge successfully.
Interpreting the Results
Although larger sample sizes improve convergence, they also require substantially data collection resources. In this example, the power for detecting the factor correlation was already high (approximately .90) with a sample size of 200. Therefore, increasing the sample size to 800 primarily improves convergence rather than statistical power.
However, it is important to consider model structure as well as sample size. The high nonconvergence rate observed at smaller sample sizes stems largely from the presence of a two-indicator factor, which provides limited information to estimate the factor’s variance and covariance accurately. To mitigate this issue, researchers are strongly encouraged to avoid two-indicator factors whenever possible and to identify or construct additional indicators for each latent variable. Doing so generally improves both convergence rates and parameter stability, reducing the need for excessively large sample sizes.
Example 2
This example demonstrates how to evaluate the power of indirect
effects in a structural equation model (SEM). Here, the effect of
f1 on f3 is mediated by f2. To
ensure that all regression coefficients are standardized, the residual
variances of the endogenous factors must be specified carefully. This is
particularly important when a dependent variable is predicted by more
than one predictor, because in such cases, the simple equation \(1 - \beta^2\) for residual variance no
longer holds.
Step 1: Determining Residual Variances
Suppose that the effect of f1 on f2 is .3,
and the effects of f1 and f2 on
f3 are .1 and −.2, respectively. First, let’s compute the
residual variance of f2.
\(1 = Var(F_2) = Var(0.3F_1 + \zeta_2) = Var(0.3F_1) + Var(\zeta_2) = 0.3^2Var(F_1) + Var(\zeta_2)\)
Because \(Var(F_1) = 1\), we have
\(Var(\zeta_2) = 1 - 0.3^2 = 0.91\)
Next, let’s compute the residual variance of f3. Since
f3 is predicted by both f1 and
f2, we need to include their covariance in the
calculation.
\(\begin{align*} 1 &= Var(F_3) \\ &= Var(0.1F_1 - 0.2F_2 + \zeta_3) \\ &= Var(0.1F_1 - 0.2F_2) + Var(\zeta_3) \\ &= Var(0.1F_1) + Var(-0.2F_2) + 2Cov(0.1F_1, -0.2F_2) + Var(\zeta_3) \\ &= 0.1^2 Var(F_1) + (-0.2)^2 Var(F_2) + 2(0.1)(-0.2)Cov(F_1,F_2) + Var(\zeta_3) \end{align*}\)
Substituting \(Var(F_1) = Var(F_2) = 1\) and \(Cov(F_1, F_2) = 0.3\), we obtain:
\(Var(\zeta_3) = 1 - 0.01 - 0.04 - (-0.012) = 0.962\)
If this calculation seems intimidating, don’t worry—later we will discuss easier ways to obtain residual variances automatically.
Step 2: Specifying the Population Model
Now that we know all the required parameters, let’s specify the
population (data-generating) model in simsem. This includes
factor loadings, residual indicator variances, regression coefficients,
and residual variances for the latent factors.
library(simsem)
genmed <- "
f1 =~ 0.6*x1 + 0.6*x2 + 0.6*x3 + 0.6*x4
f2 =~ 0.6*x5 + 0.6*x6 + 0.6*x7 + 0.6*x8
f3 =~ 0.6*x9 + 0.6*x10 + 0.6*x11
f2 ~ 0.3*f1
f3 ~ 0.1*f1 + -0.2*f2
f1 ~~ 1*f1
f2 ~~ 0.91*f2
f3 ~~ 0.962*f3
x1 ~~ 0.64*x1
x2 ~~ 0.64*x2
x3 ~~ 0.64*x3
x4 ~~ 0.64*x4
x5 ~~ 0.64*x5
x6 ~~ 0.64*x6
x7 ~~ 0.64*x7
x8 ~~ 0.64*x8
x9 ~~ 0.64*x9
x10 ~~ 0.64*x10
x11 ~~ 0.64*x11
"
dattemp <- generate(genmed, n=200)
head(dattemp)## x1 x2 x3 x4 x5 x6
## 1 1.1046490 0.30705863 0.93402291 -0.3023924 1.2350816 1.6083761
## 2 0.3188236 1.06875810 1.91989519 1.0914871 0.6654968 0.5884684
## 3 -0.6720708 0.39961185 -0.16427390 -0.1702366 -0.5058581 -0.6115593
## 4 1.2419701 0.88531192 -0.06311092 -0.1391267 1.0080673 -0.9842115
## 5 0.8222886 -0.06788384 0.41423811 0.9998163 -1.0222258 -0.7009612
## 6 2.5157736 1.59182957 -0.65234573 -0.8307726 0.1600336 1.8107035
## x7 x8 x9 x10 x11
## 1 1.7213499 0.56714083 -2.2724811 -1.86025845 -2.0076222
## 2 0.9693583 1.24877529 -0.9132921 -1.14671678 -0.2680800
## 3 -0.5752259 -0.21434695 0.3314968 -0.06498312 0.3970441
## 4 0.6452555 0.42635180 1.7420058 0.69104422 0.7620133
## 5 0.7541515 -0.63784763 -0.9270282 1.29447431 -0.8397843
## 6 -0.4061052 -0.06070242 -2.1900025 -1.28537311 0.0678565It is always a good idea to generate an example dataset to verify
that the simulated data look as expected. The resulting dataset should
not contain latent factors (f1, f2,
f3) as columns, since only indicators are simulated.
Step 3: Specifying and Fitting the Analysis Model
We now define the analysis model for testing the indirect effect and fit it to the simulated data.
medmodel <- "
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7 + x8
f3 =~ x9 + x10 + x11
f2 ~ a*f1
f3 ~ c*f1 + b*f2
ind := a*b
"
fittemp <- sem(medmodel, dattemp)
fittemp## lavaan 0.6-20 ended normally after 31 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 25
##
## Number of observations 200
##
## Model Test User Model:
##
## Test statistic 31.334
## Degrees of freedom 41
## P-value (Chi-square) 0.862Always analyze a small simulated dataset first to confirm that the model syntax is correct and that parameter estimates are as expected before running the full simulation.
Step 4: Running the Simulation
Next, we perform a Monte Carlo simulation with 1,000 replications and a sample size of 200.
simsem200 <- sim(nRep=1000, model=medmodel, n=200, generate=genmed, lavaanfun="sem")We can inspect the simulation results using the
summary() function.
summary(simsem200)## RESULT OBJECT
## Model Type
## [1] "lavaan"
## ========= Fit Indices Cutoffs ============
## Alpha
## Fit Indices 0.1 0.05 0.01 0.001 Mean SD
## chisq 54.317 58.641 66.321 76.543 42.348 9.250
## aic 6016.243 6037.980 6079.261 6109.826 5932.479 65.473
## bic 6098.701 6120.438 6161.719 6192.284 6014.936 65.473
## rmsea 0.040 0.046 0.056 0.066 0.015 0.017
## cfi 0.959 0.948 0.925 0.910 0.987 0.019
## tli 0.945 0.930 0.900 0.879 0.996 0.039
## srmr 0.054 0.057 0.062 0.066 0.047 0.006
## ========= Parameter Estimates and Standard Errors ============
## Estimate Average Estimate SD Average SE Power (Not equal 0) Std Est
## f1=~x2 1.021 0.197 0.188 1.000 0.599
## f1=~x3 1.019 0.188 0.188 1.000 0.598
## f1=~x4 1.017 0.199 0.188 1.000 0.596
## f2=~x6 1.005 0.190 0.185 1.000 0.595
## f2=~x7 1.012 0.201 0.186 1.000 0.596
## f2=~x8 1.017 0.192 0.186 1.000 0.600
## f3=~x10 1.032 0.236 0.232 1.000 0.604
## f3=~x11 1.025 0.253 0.231 1.000 0.601
## a 0.305 0.118 0.115 0.815 0.298
## c 0.109 0.129 0.121 0.137 0.105
## b -0.201 0.131 0.124 0.359 -0.198
## x1~~x1 0.632 0.090 0.085 1.000 0.639
## x2~~x2 0.631 0.088 0.085 1.000 0.637
## x3~~x3 0.631 0.087 0.085 1.000 0.637
## x4~~x4 0.634 0.087 0.085 1.000 0.640
## x5~~x5 0.629 0.087 0.085 1.000 0.634
## x6~~x6 0.633 0.085 0.084 1.000 0.642
## x7~~x7 0.637 0.085 0.085 1.000 0.640
## x8~~x8 0.632 0.086 0.085 1.000 0.636
## x9~~x9 0.628 0.105 0.102 0.996 0.632
## x10~~x10 0.624 0.101 0.102 0.998 0.628
## x11~~x11 0.624 0.106 0.101 0.992 0.632
## f1~~f1 0.360 0.094 0.094 1.000 1.000
## f2~~f2 0.331 0.094 0.087 1.000 0.901
## f3~~f3 0.345 0.102 0.104 0.999 0.938
## ind -0.062 0.048 0.045 0.109 -0.060
## Std Est SD Std Ave SE
## f1=~x2 0.071 0.065
## f1=~x3 0.070 0.065
## f1=~x4 0.070 0.065
## f2=~x6 0.067 0.065
## f2=~x7 0.067 0.065
## f2=~x8 0.066 0.065
## f3=~x10 0.080 0.079
## f3=~x11 0.083 0.079
## a 0.099 0.095
## c 0.119 0.113
## b 0.119 0.112
## x1~~x1 0.079 0.077
## x2~~x2 0.084 0.077
## x3~~x3 0.083 0.077
## x4~~x4 0.083 0.077
## x5~~x5 0.082 0.077
## x6~~x6 0.079 0.077
## x7~~x7 0.079 0.077
## x8~~x8 0.079 0.077
## x9~~x9 0.095 0.095
## x10~~x10 0.097 0.096
## x11~~x11 0.101 0.095
## f1~~f1 0.000 0.000
## f2~~f2 0.060 0.056
## f3~~f3 0.048 0.048
## ind 0.044 0.042
## ========= Correlation between Fit Indices ============
## chisq aic bic rmsea cfi tli srmr
## chisq 1.000 -0.021 -0.021 0.930 -0.899 -0.990 0.855
## aic -0.021 1.000 1.000 -0.022 0.030 0.023 -0.004
## bic -0.021 1.000 1.000 -0.022 0.030 0.023 -0.004
## rmsea 0.930 -0.022 -0.022 1.000 -0.941 -0.915 0.780
## cfi -0.899 0.030 0.030 -0.941 1.000 0.889 -0.748
## tli -0.990 0.023 0.023 -0.915 0.889 1.000 -0.854
## srmr 0.855 -0.004 -0.004 0.780 -0.748 -0.854 1.000
## ================== Replications =====================
## Number of replications = 1000
## Number of converged replications = 1000
## Number of nonconverged replications:
## 1. Nonconvergent Results = 0
## 2. Nonconvergent results from multiple imputation = 0
## 3. At least one SE were negative or NA = 0
## 4. Nonpositive-definite latent or observed (residual) covariance matrix
## (e.g., Heywood case or linear dependency) = 0
## NOTE: The data generation model is not the same as the analysis model. See the summary of the population underlying data generation by the summaryPopulation function.The results show a convergence rate of 100%, which is excellent.
However, the power of the indirect effect (a*b) is only
.109—that is, only about 109 out of 1,000 replications yielded a
statistically significant indirect effect.
Note that this power estimate is based on the Wald test, not the bootstrap test. In practice, the bootstrap test generally provides lower power, but for sample-size planning purposes, the Wald-based estimate is often sufficient as a first approximation. If needed, researchers can always increase the actual sample size later to account for differences in testing methods.
To examine parameter estimates in greater detail, we can use
summaryParam().
summaryParam(simsem200, improper=FALSE)## Estimate Average Estimate SD Average SE Power (Not equal 0)
## f1=~x2 1.02116745 0.19664597 0.18822186 1.000
## f1=~x3 1.01934032 0.18756145 0.18777942 1.000
## f1=~x4 1.01706446 0.19927468 0.18782068 1.000
## f2=~x6 1.00507090 0.18964101 0.18481650 1.000
## f2=~x7 1.01222245 0.20090293 0.18617942 1.000
## f2=~x8 1.01660112 0.19183432 0.18627152 1.000
## f3=~x10 1.03178332 0.23597507 0.23185628 1.000
## f3=~x11 1.02495064 0.25271979 0.23107382 1.000
## a 0.30502978 0.11808994 0.11480007 0.815
## c 0.10914811 0.12875494 0.12087348 0.137
## b -0.20140902 0.13052288 0.12363822 0.359
## x1~~x1 0.63196149 0.08960664 0.08464820 1.000
## x2~~x2 0.63054411 0.08839002 0.08483803 1.000
## x3~~x3 0.63096058 0.08745519 0.08472041 1.000
## x4~~x4 0.63352680 0.08665634 0.08479024 1.000
## x5~~x5 0.62864019 0.08731376 0.08459574 1.000
## x6~~x6 0.63311478 0.08537712 0.08412928 1.000
## x7~~x7 0.63661059 0.08462874 0.08485429 1.000
## x8~~x8 0.63183166 0.08592624 0.08467656 1.000
## x9~~x9 0.62773431 0.10486424 0.10160691 0.996
## x10~~x10 0.62383316 0.10140649 0.10205279 0.998
## x11~~x11 0.62391778 0.10630735 0.10133892 0.992
## f1~~f1 0.36008838 0.09375476 0.09400615 1.000
## f2~~f2 0.33073259 0.09409080 0.08727893 1.000
## f3~~f3 0.34464110 0.10212251 0.10436825 0.999
## ind -0.06190609 0.04837481 0.04544080 0.109
## Std Est Std Est SD Std Ave SE
## f1=~x2 0.59870210 0.0707515478776193224020 0.06519651
## f1=~x3 0.59849662 0.0696942773552343547427 0.06516415
## f1=~x4 0.59602725 0.0702214797955686836906 0.06529239
## f2=~x6 0.59486372 0.0669674226461103166974 0.06505880
## f2=~x7 0.59597321 0.0672674154921543004937 0.06507062
## f2=~x8 0.59956462 0.0662644131736266905897 0.06490721
## f3=~x10 0.60434367 0.0804676696139831337806 0.07919957
## f3=~x11 0.60097425 0.0832823748200491320581 0.07918687
## a 0.29780403 0.0991873656033078338057 0.09504398
## c 0.10505619 0.1189607509053121903175 0.11323486
## b -0.19845945 0.1191803701305017837431 0.11222265
## x1~~x1 0.63873320 0.0789330775727356270810 0.07741571
## x2~~x2 0.63655502 0.0843295547124761624147 0.07738356
## x3~~x3 0.63694936 0.0827510468881289124043 0.07731245
## x4~~x4 0.63982539 0.0829794175404538619789 0.07715952
## x5~~x5 0.63408752 0.0821508896238050817518 0.07733606
## x6~~x6 0.64165700 0.0794089214012806132725 0.07676816
## x7~~x7 0.64029555 0.0787414105538081915414 0.07696559
## x8~~x8 0.63613569 0.0791026060602835617841 0.07723236
## x9~~x9 0.63227756 0.0950223199370531512731 0.09508932
## x10~~x10 0.62830016 0.0971673339888071418979 0.09550690
## x11~~x11 0.63190093 0.1008868437390423311584 0.09508585
## f1~~f1 1.00000000 0.0000000000000001459311 0.00000000
## f2~~f2 0.90148447 0.0600665072763092444674 0.05588842
## f3~~f3 0.93763636 0.0484203700247180096605 0.04799093
## ind -0.05994133 0.0441309693613317552408 0.04212399Step 5: Varying Sample Size
Clearly, we need a larger sample size to achieve adequate power. In
simsem, the argument n can be specified as a
vector, allowing the simulation to run across multiple sample sizes in a
single call. Each replication corresponds to one value of
n, and the total number of replications equals the length
of the sample-size vector. Here, we vary n from 200 to
1200, increasing by 1.
simsemvary <- sim(model=medmodel, n=200:1200, generate=genmed, lavaanfun="sem")Step 6: Plotting and Estimating Required Sample Size
We can visualize how the power of the indirect effect changes with
sample size using the plotPower() function. The argument
powerParam specifies which parameter’s power is to be
plotted. You can check the exact parameter names by inspecting the row
names in summaryParam() output.
plotPower(simsemvary, powerParam="ind")
To determine the estimated sample size required to achieve .80 power,
we first extract the predicted power values across sample sizes using
getPower(), and then use findPower().
pow <- getPower(simsemvary)
findPower(pow, "N", 0.80)## f1=~x2 f1=~x3 f1=~x4 f2=~x6 f2=~x7 f2=~x8 f3=~x10 f3=~x11
## Inf Inf Inf Inf Inf Inf Inf Inf
## a c b x1~~x1 x2~~x2 x3~~x3 x4~~x4 x5~~x5
## 205 NA 551 Inf Inf Inf Inf Inf
## x6~~x6 x7~~x7 x8~~x8 x9~~x9 x10~~x10 x11~~x11 f1~~f1 f2~~f2
## Inf Inf Inf Inf Inf Inf Inf Inf
## f3~~f3 ind
## Inf 644The result indicates that a sample size of approximately 644 is needed to achieve a power of .80 for the indirect effect.
Step 7: Confirming the Result
Finally, let’s rerun the simulation with a rounded sample size of 650 to verify that the desired power level is achieved.
simsem650 <- sim(nRep=1000, model=medmodel, n=650, generate=genmed, lavaanfun="sem")summaryParam(simsem650, improper=FALSE)["ind", "Power (Not equal 0)"]## [1] 0.824The resulting power for the indirect effect is approximately .824, which meets the conventional .80 criterion. Again, keep in mind that this estimate is based on the Wald test (delta method/Sobel test); in practice, a bootstrap test will yield slightly lower power. Therefore, it would be prudent to plan for an actual sample size of around 700 to ensure sufficient power in real data analysis.
Example 2 using simstandard
Let’s revisit how to calculate residual variances. When there are two
or more predictors—as in regression models or models with
cross-loadings—the simple residual variance formula \(1 - \beta^2\) no longer applies. In such
cases, a more reliable and flexible approach is to generate population
data using the simstandard package. Then, the
sim() function from simsem can randomly draw
samples from that population to perform the simulation study.
In this example, we use the same mediation model as before, but we do
not specify any factor or indicator variances manually. Instead, we use
the sim_standardized() function to generate population data
with one million cases (n = 1,000,000). In this dataset,
only the observed variables are saved—latent factors and latent
indicator residuals are not included.
library(simstandard)
set.seed(123321)
genmodel <- "
f1 =~ 0.6*x1 + 0.6*x2 + 0.6*x3 + 0.6*x4
f2 =~ 0.6*x5 + 0.6*x6 + 0.6*x7 + 0.6*x8
f3 =~ 0.6*x9 + 0.6*x10 + 0.6*x11
f2 ~ 0.3*f1
f3 ~ 0.1*f1 + -0.2*f2
"
popdat <- sim_standardized(genmodel, n = 1000000, observed=TRUE,
latent=FALSE, errors=FALSE)
head(popdat)## # A tibble: 6 × 11
## x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.98 -0.447 0.645 -0.374 0.552 0.532 -0.338 -0.184 0.765 -0.0842
## 2 -2.55 -1.27 0.120 0.563 -1.43 -0.870 -0.595 0.0710 -0.322 -0.616
## 3 -0.732 0.0824 0.186 0.319 -0.117 -0.576 -0.0226 -0.733 1.32 0.121
## 4 -1.23 -1.98 -0.211 -0.0946 1.15 0.899 2.11 0.834 0.848 0.748
## 5 0.433 0.183 0.0466 -0.342 1.55 0.216 0.146 -0.473 1.80 -0.181
## 6 0.642 -0.167 -1.76 -0.236 0.139 -0.884 -0.289 -0.694 0.330 -1.24
## # ℹ 1 more variable: x11 <dbl>Next, we define the analysis model, just as in previous examples. Before running a simulation, it’s always a good practice to fit the model to the population data and inspect the parameter estimates. This step ensures that the population model behaves as expected.
medmodel <- "
f1 =~ x1 + x2 + x3 + x4
f2 =~ x5 + x6 + x7 + x8
f3 =~ x9 + x10 + x11
f2 ~ a*f1
f3 ~ c*f1 + b*f2
ind := a*b
"
paramout <- sem(medmodel, popdat)
parameterEstimates(paramout)## lhs op rhs label est se z pvalue ci.lower ci.upper
## 1 f1 =~ x1 1.000 0.000 NA NA 1.000 1.000
## 2 f1 =~ x2 0.999 0.003 392.048 0 0.994 1.004
## 3 f1 =~ x3 0.999 0.003 392.082 0 0.994 1.004
## 4 f1 =~ x4 0.997 0.003 391.711 0 0.992 1.002
## 5 f2 =~ x5 1.000 0.000 NA NA 1.000 1.000
## 6 f2 =~ x6 0.999 0.003 393.544 0 0.994 1.004
## 7 f2 =~ x7 0.998 0.003 393.273 0 0.993 1.003
## 8 f2 =~ x8 0.999 0.003 393.628 0 0.994 1.004
## 9 f3 =~ x9 1.000 0.000 NA NA 1.000 1.000
## 10 f3 =~ x10 0.997 0.003 322.516 0 0.991 1.003
## 11 f3 =~ x11 0.998 0.003 322.521 0 0.992 1.004
## 12 f2 ~ f1 a 0.297 0.002 189.237 0 0.294 0.300
## 13 f3 ~ f1 c 0.099 0.002 60.351 0 0.096 0.103
## 14 f3 ~ f2 b -0.201 0.002 -118.107 0 -0.204 -0.198
## 15 x1 ~~ x1 0.640 0.001 528.132 0 0.637 0.642
## 16 x2 ~~ x2 0.639 0.001 528.344 0 0.637 0.641
## 17 x3 ~~ x3 0.639 0.001 528.231 0 0.636 0.641
## 18 x4 ~~ x4 0.641 0.001 529.439 0 0.638 0.643
## 19 x5 ~~ x5 0.638 0.001 529.668 0 0.636 0.641
## 20 x6 ~~ x6 0.639 0.001 530.189 0 0.637 0.641
## 21 x7 ~~ x7 0.641 0.001 531.042 0 0.638 0.643
## 22 x8 ~~ x8 0.638 0.001 529.923 0 0.636 0.641
## 23 x9 ~~ x9 0.638 0.001 450.714 0 0.635 0.641
## 24 x10 ~~ x10 0.640 0.001 453.155 0 0.637 0.643
## 25 x11 ~~ x11 0.642 0.001 453.249 0 0.639 0.644
## 26 f1 ~~ f1 0.361 0.001 266.821 0 0.359 0.364
## 27 f2 ~~ f2 0.329 0.001 262.577 0 0.327 0.332
## 28 f3 ~~ f3 0.347 0.001 232.346 0 0.344 0.350
## 29 ind := a*b ind -0.060 0.001 -101.236 0 -0.061 -0.059By examining these estimates, you can confirm that all parameter values (factor loadings, regression paths, and indirect effects) align with those specified in your conceptual model.
Now that the population dataset is ready, we can use the
sim() function to run the simulation. Instead of specifying
a generate model as before, we use the rawData argument to
feed in the pre-generated population dataset.
simsem200b <- sim(nRep=1000, model=medmodel, n=650, rawData=popdat, lavaanfun="sem")Although we specify 1,000 replications, you might notice that simsem
actually runs 1,001 replications. The extra replication is used
internally to analyze the target model on the population data itself,
allowing simsem to store the “true” parameter values for
comparison.
summary(simsem200b)## RESULT OBJECT
## Model Type
## [1] "lavaan"
## ========= Fit Indices Cutoffs ============
## Alpha
## Fit Indices 0.1 0.05 0.01 0.001 Mean SD
## chisq 54.203 56.871 64.156 70.685 41.633 9.126
## aic 19400.249 19443.993 19509.593 19586.589 19250.027 119.211
## bic 19512.173 19555.917 19621.517 19698.513 19361.951 119.211
## rmsea 0.022 0.024 0.029 0.033 0.008 0.009
## cfi 0.988 0.985 0.980 0.974 0.996 0.005
## tli 0.984 0.980 0.973 0.965 0.999 0.012
## srmr 0.030 0.032 0.034 0.037 0.026 0.003
## ========= Parameter Estimates and Standard Errors ============
## Estimate Average Estimate SD Average SE Power (Not equal 0) Std Est
## f1=~x2 1.009 0.101 0.101 1.000 0.603
## f1=~x3 1.002 0.100 0.101 1.000 0.599
## f1=~x4 1.005 0.104 0.101 1.000 0.599
## f2=~x6 1.007 0.104 0.101 1.000 0.601
## f2=~x7 1.004 0.102 0.101 1.000 0.600
## f2=~x8 1.006 0.102 0.101 1.000 0.601
## f3=~x10 1.011 0.128 0.124 1.000 0.601
## f3=~x11 1.010 0.128 0.124 1.000 0.601
## a 0.300 0.062 0.062 1.000 0.299
## c 0.103 0.066 0.065 0.339 0.102
## b -0.204 0.068 0.067 0.883 -0.203
## x1~~x1 0.637 0.048 0.047 1.000 0.638
## x2~~x2 0.635 0.046 0.047 1.000 0.636
## x3~~x3 0.638 0.047 0.047 1.000 0.640
## x4~~x4 0.639 0.048 0.047 1.000 0.639
## x5~~x5 0.636 0.047 0.047 1.000 0.638
## x6~~x6 0.637 0.047 0.047 1.000 0.638
## x7~~x7 0.638 0.046 0.047 1.000 0.639
## x8~~x8 0.636 0.046 0.047 1.000 0.638
## x9~~x9 0.635 0.058 0.056 1.000 0.638
## x10~~x10 0.636 0.056 0.056 1.000 0.637
## x11~~x11 0.636 0.058 0.056 1.000 0.637
## f1~~f1 0.362 0.053 0.053 1.000 1.000
## f2~~f2 0.329 0.051 0.049 1.000 0.908
## f3~~f3 0.345 0.062 0.058 1.000 0.954
## ind -0.061 0.025 0.024 0.823 -0.061
## Std Est SD Std Ave SE Average Param Average Bias Coverage
## f1=~x2 0.036 0.036 0.999 0.010 0.949
## f1=~x3 0.037 0.036 0.999 0.003 0.954
## f1=~x4 0.038 0.036 0.997 0.007 0.955
## f2=~x6 0.035 0.036 0.999 0.007 0.944
## f2=~x7 0.035 0.036 0.998 0.006 0.943
## f2=~x8 0.036 0.036 0.999 0.007 0.954
## f3=~x10 0.044 0.044 0.997 0.014 0.945
## f3=~x11 0.045 0.044 0.998 0.012 0.947
## a 0.055 0.053 0.297 0.003 0.959
## c 0.064 0.063 0.099 0.004 0.957
## b 0.064 0.063 -0.201 -0.003 0.945
## x1~~x1 0.044 0.043 0.640 -0.002 0.937
## x2~~x2 0.043 0.043 0.639 -0.004 0.956
## x3~~x3 0.045 0.043 0.639 -0.001 0.943
## x4~~x4 0.045 0.043 0.641 -0.001 0.936
## x5~~x5 0.045 0.043 0.638 -0.002 0.948
## x6~~x6 0.043 0.043 0.639 -0.002 0.955
## x7~~x7 0.042 0.043 0.641 -0.003 0.955
## x8~~x8 0.043 0.043 0.638 -0.002 0.950
## x9~~x9 0.055 0.052 0.638 -0.003 0.942
## x10~~x10 0.053 0.052 0.640 -0.004 0.941
## x11~~x11 0.054 0.052 0.642 -0.006 0.942
## f1~~f1 0.000 0.000 0.361 0.001 0.939
## f2~~f2 0.033 0.032 0.329 0.000 0.931
## f3~~f3 0.025 0.025 0.347 -0.002 0.924
## ind 0.023 0.023 -0.060 -0.001 0.932
## ========= Correlation between Fit Indices ============
## chisq aic bic rmsea cfi tli srmr
## chisq 1.000 0.002 0.002 0.927 -0.895 -0.998 0.871
## aic 0.002 1.000 1.000 0.017 -0.004 -0.002 -0.002
## bic 0.002 1.000 1.000 0.017 -0.004 -0.002 -0.002
## rmsea 0.927 0.017 0.017 1.000 -0.954 -0.923 0.790
## cfi -0.895 -0.004 -0.004 -0.954 1.000 0.891 -0.762
## tli -0.998 -0.002 -0.002 -0.923 0.891 1.000 -0.870
## srmr 0.871 -0.002 -0.002 0.790 -0.762 -0.870 1.000
## ================== Replications =====================
## Number of replications = 1000
## Number of converged replications = 1000
## Number of nonconverged replications:
## 1. Nonconvergent Results = 0
## 2. Nonconvergent results from multiple imputation = 0
## 3. At least one SE were negative or NA = 0
## 4. Nonpositive-definite latent or observed (residual) covariance matrix
## (e.g., Heywood case or linear dependency) = 0The output from this simulation should look very similar to the
results obtained earlier when specifying the data-generating model in
genmodel. The key difference here is that you don’t need to
calculate residual variances manually—simstandard takes
care of that process automatically when generating population data.
With this approach, you can avoid tedious residual variance calculations altogether. This method not only saves time but also reduces the chance of mistakes when dealing with multiple predictors or correlated latent variables. It is therefore highly recommended whenever your population model involves complex regression paths or cross-loadings.
Example 3 Cohen’s d
In this example, suppose we have two factors measuring the same construct but rated by different sources—for instance, two groups of raters evaluating the same targets. Let \(F_1\) represent ratings from standardized raters, and \(F_2\) represent ratings from newly-trained raters.
An assessment company would like to determine how many trials are required to ensure that both rater groups produce statistically equivalent ratings. The target is to achieve high statistical power to detect a serious difference first (defined as Cohen’s d = 0.5). Later in this tutorial, we will also demonstrate an approach using the width of the confidence interval as the evaluation criterion.
Step 1: Define the Population Model
Assume the standardized raters have four standardized items with loadings of 0.8, 0.7, 0.6, and 0.5, respectively. The corresponding residual variances are therefore 0.36, 0.51, 0.64, and 0.75 (computed as \(1 - \lambda^2\)). The factor variance for standardized raters (\(F_1\)) is fixed at 1.
For the newly-trained raters, assume they share the same unstandardized loadings and residual variances, but their factor variance is larger at 1.2, indicating slightly more variability among them. The correlation between \(F_1\) and \(F_2\) is set at 0.6. Hence, their factor covariance is:
\(0.6\times \sqrt{1} \times \sqrt{1.2} = .657\).
To represent shared measurement context, assume error correlations of 0.3 between the same indicators across raters. Thus, the error covariances between corresponding indicators are:
\(\begin{align*} 0.3\times \sqrt{0.36 \times 0.36} &= 0.108 \\ 0.3\times \sqrt{0.51 \times 0.51} &= 0.153 \\ 0.3\times \sqrt{0.64 \times 0.64} &= 0.192 \\ 0.3\times \sqrt{0.75 \times 0.75} &= 0.225 \end{align*}\)
Step 2: Specify Means and Intercepts
Suppose the indicator means of the standardized raters are 3.3, 2.5, 4.2, and 3.7. We fix the factor mean of the standardized raters (\(F_1\)) at 0, making the indicator intercepts equal to those means. For simplicity, assume equal indicator intercepts across both rater groups.
We wish to detect a factor mean difference of 0.5 on the standardized scale. Given that the variance of \(F_1\) is 1, a Cohen’s d of 0.5 corresponds to a factor mean of 0.5 for the newly-trained raters (\(F_2\)).
Step 3: Generate Population Data
The complete population model is defined below:
gencohend <- "
f1 =~ 0.8*x1 + 0.7*x2 + 0.6*x3 + 0.5*x4
f2 =~ 0.8*x5 + 0.7*x6 + 0.6*x7 + 0.5*x8
f2 ~~ 0.657*f1
f1 ~~ 1*f1
f2 ~~ 1.2*f2
x1 ~~ 0.108*x5
x2 ~~ 0.153*x6
x3 ~~ 0.192*x7
x4 ~~ 0.225*x8
x1 ~~ 0.36*x1
x2 ~~ 0.51*x2
x3 ~~ 0.64*x3
x4 ~~ 0.75*x4
x5 ~~ 0.36*x5
x6 ~~ 0.51*x6
x7 ~~ 0.64*x7
x8 ~~ 0.75*x8
f1 ~ 0*1
f2 ~ 0.5*1
x1 ~ 3.3*1
x2 ~ 2.5*1
x3 ~ 4.2*1
x4 ~ 3.7*1
x5 ~ 3.3*1
x6 ~ 2.5*1
x7 ~ 4.2*1
x8 ~ 3.7*1
"
dattemp <- generate(gencohend, n=200)
head(dattemp)## x1 x2 x3 x4 x5 x6 x7 x8
## 1 3.727076 2.055761 5.134956 3.421463 2.912067 3.069674 4.219464 4.046194
## 2 4.804308 2.875282 5.224742 4.473247 3.964362 2.322858 4.167434 3.393640
## 3 5.052443 3.079578 5.031502 3.384489 3.791693 2.997792 4.948670 3.067496
## 4 4.062140 1.846778 4.854868 2.831205 4.709710 1.374933 5.830632 3.977962
## 5 4.450930 3.385976 3.362921 2.014782 2.337866 1.069573 2.950040 1.442557
## 6 3.915005 4.477268 4.611444 2.198267 4.257015 4.092424 4.458248 3.296052Note that The simstandard package cannot be used in this
example because we must specify unstandardized
parameters for the newly-trained raters’ factor. Be extremely
careful when defining your data-generating model—any mistake here will
distort your simulated power or other estimated properties.
Step 4: Define the Analysis Model
The analysis model represents a scalar invariance model with
correlated residuals. Cohen’s d is computed as the difference
between the factor means (M2 - M1) divided by the square
root of the variance of \(F_1\)
(sqrt(V1)).
The marker-variable approach is used for model identification, meaning the first indicator’s loading is fixed to 1 for each factor. Keep in mind that the indicator intercepts in this model differ from those in the population model due to this parameterization.
cohendmodel <- "
f1 =~ 1*x1 + L2*x2 + L3*x3 + L4*x4
f2 =~ 1*x5 + L2*x6 + L3*x7 + L4*x8
f1 ~~ NA*f2
f1 ~~ NA*f1 + V1*f1
f2 ~~ NA*f2 + V2*f2
x1 ~~ x5
x2 ~~ x6
x3 ~~ x7
x4 ~~ x8
x1 ~ 0*1
x2 ~ I2*1
x3 ~ I3*1
x4 ~ I4*1
x5 ~ 0*1
x6 ~ I2*1
x7 ~ I3*1
x8 ~ I4*1
f1 ~ M1*1
f2 ~ M2*1
eff := (M2 - M1)/sqrt(V1)
"
fittemp <- cfa(cohendmodel, dattemp)
fittemp## lavaan 0.6-20 ended normally after 54 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 29
## Number of equality constraints 6
##
## Number of observations 200
##
## Model Test User Model:
##
## Test statistic 21.464
## Degrees of freedom 21
## P-value (Chi-square) 0.431Step 5: Run the Simulation
We first run a simulation with 200 participants per group:
simcohend200 <- sim(nRep=1000, model=cohendmodel, generate=gencohend, n=200, lavaanfun="cfa")Then, we summarize the results:
summary(simcohend200)## RESULT OBJECT
## Model Type
## [1] "lavaan"
## ========= Fit Indices Cutoffs ============
## Alpha
## Fit Indices 0.1 0.05 0.01 0.001 Mean SD
## chisq 30.260 33.763 39.499 46.372 21.835 6.453
## aic 4178.912 4203.246 4231.326 4286.706 4108.625 56.892
## bic 4254.773 4279.108 4307.187 4362.567 4184.487 56.892
## rmsea 0.047 0.055 0.066 0.078 0.017 0.020
## cfi 0.982 0.976 0.964 0.955 0.994 0.008
## tli 0.975 0.968 0.952 0.940 0.998 0.017
## srmr 0.051 0.054 0.060 0.065 0.041 0.008
## ========= Parameter Estimates and Standard Errors ============
## Estimate Average Estimate SD Average SE Power (Not equal 0)
## L2 <- (f1=~x2) 0.876 0.080 0.072 1.000
## L3 <- (f1=~x3) 0.748 0.072 0.069 1.000
## L4 <- (f1=~x4) 0.628 0.072 0.068 1.000
## L2 <- (f2=~x6) 0.876 0.080 0.072 1.000
## L3 <- (f2=~x7) 0.748 0.072 0.069 1.000
## L4 <- (f2=~x8) 0.628 0.072 0.068 1.000
## f1~~f2 0.417 0.081 0.077 1.000
## V1 0.638 0.099 0.096 1.000
## V2 0.766 0.117 0.112 1.000
## x1~~x5 0.105 0.050 0.047 0.606
## x2~~x6 0.153 0.050 0.049 0.893
## x3~~x7 0.192 0.054 0.054 0.968
## x4~~x8 0.219 0.059 0.059 0.977
## I2 <- (x2~1) -0.392 0.285 0.259 0.328
## I3 <- (x3~1) 1.733 0.258 0.247 1.000
## I4 <- (x4~1) 1.629 0.255 0.244 1.000
## I2 <- (x6~1) -0.392 0.285 0.259 0.328
## I3 <- (x7~1) 1.733 0.258 0.247 1.000
## I4 <- (x8~1) 1.629 0.255 0.244 1.000
## M1 3.300 0.070 0.068 1.000
## M2 3.699 0.073 0.073 1.000
## x1~~x1 0.355 0.067 0.065 1.000
## x2~~x2 0.508 0.069 0.068 1.000
## x3~~x3 0.635 0.076 0.074 1.000
## x4~~x4 0.744 0.083 0.081 1.000
## x5~~x5 0.355 0.070 0.069 0.998
## x6~~x6 0.508 0.071 0.069 1.000
## x7~~x7 0.640 0.074 0.075 1.000
## x8~~x8 0.744 0.083 0.081 1.000
## eff 0.504 0.088 0.087 1.000
## [ L2 ] - [ L2 ] 0.000 0.000 NA 1.000
## [ L3 ] - [ L3 ] 0.000 0.000 NA 1.000
## [ L4 ] - [ L4 ] 0.000 0.000 NA 1.000
## [ I2 ] - [ I2 ] 0.000 0.000 NA 1.000
## [ I3 ] - [ I3 ] 0.000 0.000 NA 1.000
## [ I4 ] - [ I4 ] 0.000 0.000 NA 1.000
## Std Est Std Est SD Std Ave SE
## L2 <- (f1=~x2) 0.697 0.043 0.042
## L3 <- (f1=~x3) 0.597 0.047 0.045
## L4 <- (f1=~x4) 0.500 0.051 0.048
## L2 <- (f2=~x6) 0.729 0.043 0.040
## L3 <- (f2=~x7) 0.630 0.045 0.044
## L4 <- (f2=~x8) 0.534 0.051 0.049
## f1~~f2 0.595 0.065 0.062
## V1 1.000 0.000 0.000
## V2 1.000 0.000 0.000
## x1~~x5 0.290 0.119 0.112
## x2~~x6 0.300 0.085 0.082
## x3~~x7 0.300 0.073 0.073
## x4~~x8 0.294 0.069 0.069
## I2 <- (x2~1) -0.391 0.283 0.256
## I3 <- (x3~1) 1.751 0.300 0.286
## I4 <- (x4~1) 1.642 0.291 0.279
## I2 <- (x6~1) -0.372 0.268 0.244
## I3 <- (x7~1) 1.688 0.292 0.279
## I4 <- (x8~1) 1.602 0.287 0.276
## M1 4.170 0.349 0.334
## M2 4.264 0.337 0.322
## x1~~x1 0.359 0.065 0.063
## x2~~x2 0.512 0.060 0.058
## x3~~x3 0.641 0.056 0.054
## x4~~x4 0.747 0.051 0.048
## x5~~x5 0.319 0.064 0.061
## x6~~x6 0.467 0.062 0.058
## x7~~x7 0.601 0.057 0.056
## x8~~x8 0.712 0.055 0.052
## eff 0.094 0.338 0.332
## [ L2 ] - [ L2 ] -0.032 0.042 0.041
## [ L3 ] - [ L3 ] -0.033 0.042 0.041
## [ L4 ] - [ L4 ] -0.034 0.040 0.040
## [ I2 ] - [ I2 ] -0.019 0.031 0.029
## [ I3 ] - [ I3 ] 0.063 0.100 0.096
## [ I4 ] - [ I4 ] 0.040 0.095 0.094
## ========= Correlation between Fit Indices ============
## chisq aic bic rmsea cfi tli srmr
## chisq 1.000 -0.005 -0.005 0.930 -0.898 -0.993 0.782
## aic -0.005 1.000 1.000 -0.002 -0.004 0.008 0.029
## bic -0.005 1.000 1.000 -0.002 -0.004 0.008 0.029
## rmsea 0.930 -0.002 -0.002 1.000 -0.938 -0.922 0.712
## cfi -0.898 -0.004 -0.004 -0.938 1.000 0.897 -0.679
## tli -0.993 0.008 0.008 -0.922 0.897 1.000 -0.785
## srmr 0.782 0.029 0.029 0.712 -0.679 -0.785 1.000
## ================== Replications =====================
## Number of replications = 1000
## Number of converged replications = 1000
## Number of nonconverged replications:
## 1. Nonconvergent Results = 0
## 2. Nonconvergent results from multiple imputation = 0
## 3. At least one SE were negative or NA = 0
## 4. Nonpositive-definite latent or observed (residual) covariance matrix
## (e.g., Heywood case or linear dependency) = 0
## NOTE: The data generation model is not the same as the analysis model. See the summary of the population underlying data generation by the summaryPopulation function.The power to detect the effect size (Cohen’s d = 0.5) is 1.000, indicating that 200 trials are more than sufficient. The convergence rate is also 100%. We can therefore reduce the sample size to find the smallest number that still yields adequate power.
Step 6: Vary Sample Size
Next, we run a simulation varying sample sizes from 20 to 200:
simcohendvary <- sim(model=cohendmodel, generate=gencohend, n=20:200, lavaanfun="cfa")We can visualize the relationship between sample size and power:
plotPower(simcohendvary, powerParam="eff")
To determine the sample size needed for 80% power, we use
getPower() and findPower():
pow <- getPower(simcohendvary)
findPower(pow, "N", 0.80)## L2 <- (f1=~x2) L3 <- (f1=~x3) L4 <- (f1=~x4) L2 <- (f2=~x6) L3 <- (f2=~x7)
## Inf Inf 20 Inf Inf
## L4 <- (f2=~x8) f1~~f2 V1 V2 x1~~x5
## 20 39 22 22 NA
## x2~~x6 x3~~x7 x4~~x8 I2 <- (x2~1) I3 <- (x3~1)
## 139 140 126 NA 53
## I4 <- (x4~1) I2 <- (x6~1) I3 <- (x7~1) I4 <- (x8~1) M1
## 40 NA 53 40 Inf
## M2 x1~~x1 x2~~x2 x3~~x3 x4~~x4
## Inf 34 29 Inf Inf
## x5~~x5 x6~~x6 x7~~x7 x8~~x8 eff
## 41 29 22 Inf 46
## [ L2 ] - [ L2 ] [ L3 ] - [ L3 ] [ L4 ] - [ L4 ] [ I2 ] - [ I2 ] [ I3 ] - [ I3 ]
## Inf Inf Inf Inf Inf
## [ I4 ] - [ I4 ]
## InfThe estimated sample size is 46.
Step 7: Verify with 50 participants
Finally, let’s confirm that 50 participants achieve sufficient power:
simcohend50 <- sim(nRep=1000, model=cohendmodel, generate=gencohend, n=50, lavaanfun="cfa")summaryParam(simcohend50, improper=FALSE)["eff", "Power (Not equal 0)"]## [1] 0.8787565The resulting power is 0.878, indicating that 50 participants (or trials) provide adequate power to detect a Cohen’s d = 0.5 difference between the standardized and newly-trained raters.
Next, we check the model convergence rate to ensure that the estimation is stable.
summaryConverge(simcohend50)## $Converged
## num.converged num.nonconverged
## 965 35
##
## $`Nonconvergent Reasons`
## count
## Nonconvergent 0
## Improper SE 0
## Nonpositive definite matrix 35
## Optimal estimates were not guaranteed 0The convergence rate is 96.5%, which is acceptable and supports retaining 50 trials as the recommended sample size for achieving both sufficient power and model stability.
Power Analysis in Model Fit Evaluation in General
Let’s now turn to model fit evaluation — examining whether fit indices can effectively reject a misspecified model.
Step 1: Define the Population Model
We begin with a two-factor model, where each factor has standardized loadings of 0.7 and the two factors are correlated at 0.9.
library(simsem)
twofactormodel <- "
f1 =~ 0.7*x1 + 0.7*x2 + 0.7*x3
f2 =~ 0.7*x4 + 0.7*x5 + 0.7*x6
f2 ~~ 0.9*f1
f1 ~~ 1*f1
f2 ~~ 1*f2
x1 ~~ 0.51*x1
x2 ~~ 0.51*x2
x3 ~~ 0.51*x3
x4 ~~ 0.51*x4
x5 ~~ 0.51*x5
x6 ~~ 0.51*x6
"
dattemp <- generate(twofactormodel, n=200)
head(dattemp)## x1 x2 x3 x4 x5 x6
## 1 -0.7177618 0.007489944 0.004931526 1.74349396 1.03006442 -0.3459891
## 2 0.1969952 0.539139271 1.458209314 0.20094467 1.27915887 1.1237937
## 3 -0.4989032 0.970419698 0.257101636 0.09355938 1.20593765 0.5308583
## 4 0.5049516 0.331151059 0.745602836 -0.59066632 0.03417567 -0.6913791
## 5 0.3496283 -0.082953929 0.156229878 -1.58733121 -0.93516006 0.2713161
## 6 -1.0231527 -1.147991315 -1.948349524 -1.21957341 -0.53685906 0.0590555Step 2: Define the Misspecified Analysis Model
Next, we analyze the data using a one-factor CFA model, which incorrectly assumes that all six items load on a single latent variable.
onefactormodel <- "
f1 =~ x1 + x2 + x3 + x4 + x5 + x6
"
fittemp <- cfa(onefactormodel, data=dattemp)
fittemp## lavaan 0.6-20 ended normally after 19 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 12
##
## Number of observations 200
##
## Model Test User Model:
##
## Test statistic 15.296
## Degrees of freedom 9
## P-value (Chi-square) 0.083Step 3: Run the Simulation
We now simulate 1,000 replications, where the population data are generated from the two-factor model but analyzed using the one-factor model, with a sample size of 200.
simtwotoone200 <- sim(nRep=1000, model=onefactormodel, generate=twofactormodel, n=200, lavaanfun="cfa")Step 4: Summarize the Results
summary(simtwotoone200)## RESULT OBJECT
## Model Type
## [1] "lavaan"
## ========= Fit Indices Cutoffs ============
## Alpha
## Fit Indices 0.1 0.05 0.01 0.001 Mean SD
## chisq 23.112 26.791 32.960 38.952 14.688 6.560
## aic 3098.723 3114.683 3146.533 3186.958 3034.011 49.769
## bic 3138.303 3154.263 3186.112 3226.537 3073.591 49.769
## rmsea 0.089 0.099 0.115 0.129 0.049 0.033
## cfi 0.960 0.953 0.937 0.919 0.984 0.016
## tli 0.934 0.922 0.896 0.864 0.975 0.029
## srmr 0.042 0.045 0.051 0.055 0.032 0.008
## ========= Parameter Estimates and Standard Errors ============
## Estimate Average Estimate SD Average SE Power (Not equal 0) Std Est
## f1=~x2 1.005 0.125 0.125 1 0.676
## f1=~x3 1.004 0.125 0.126 1 0.674
## f1=~x4 1.011 0.137 0.126 1 0.678
## f1=~x5 1.005 0.137 0.126 1 0.675
## f1=~x6 1.008 0.138 0.126 1 0.676
## x1~~x1 0.533 0.069 0.064 1 0.539
## x2~~x2 0.533 0.065 0.064 1 0.541
## x3~~x3 0.537 0.065 0.064 1 0.543
## x4~~x4 0.533 0.064 0.064 1 0.539
## x5~~x5 0.535 0.066 0.064 1 0.542
## x6~~x6 0.536 0.066 0.064 1 0.541
## f1~~f1 0.460 0.095 0.091 1 1.000
## Std Est SD Std Ave SE
## f1=~x2 0.049 0.047
## f1=~x3 0.048 0.047
## f1=~x4 0.048 0.046
## f1=~x5 0.047 0.047
## f1=~x6 0.049 0.047
## x1~~x1 0.066 0.063
## x2~~x2 0.066 0.063
## x3~~x3 0.064 0.063
## x4~~x4 0.064 0.063
## x5~~x5 0.063 0.063
## x6~~x6 0.065 0.063
## f1~~f1 0.000 0.000
## ========= Correlation between Fit Indices ============
## chisq aic bic rmsea cfi tli srmr
## chisq 1.000 0.038 0.038 0.966 -0.968 -0.984 0.949
## aic 0.038 1.000 1.000 0.043 -0.028 -0.032 0.042
## bic 0.038 1.000 1.000 0.043 -0.028 -0.032 0.042
## rmsea 0.966 0.043 0.043 1.000 -0.929 -0.955 0.936
## cfi -0.968 -0.028 -0.028 -0.929 1.000 0.984 -0.942
## tli -0.984 -0.032 -0.032 -0.955 0.984 1.000 -0.972
## srmr 0.949 0.042 0.042 0.936 -0.942 -0.972 1.000
## ================== Replications =====================
## Number of replications = 1000
## Number of converged replications = 1000
## Number of nonconverged replications:
## 1. Nonconvergent Results = 0
## 2. Nonconvergent results from multiple imputation = 0
## 3. At least one SE were negative or NA = 0
## 4. Nonpositive-definite latent or observed (residual) covariance matrix
## (e.g., Heywood case or linear dependency) = 0
## NOTE: The data generation model is not the same as the analysis model. See the summary of the population underlying data generation by the summaryPopulation function.All replications converge successfully. In the “Fit Indices Cutoff” section, we observe the average fit indices: RMSEA = .049, CFI = .984, TLI = .975, SRMR = .032. Surprisingly, these values look quite good even though the model is misspecified, meaning that the power to reject the model is likely low.
Step 5: Apply Hu & Bentler (1999) Cutoffs
hubentler99 <- c(RMSEA=0.05, CFI=0.95, TLI=0.95, SRMR=0.06)Now we can visualize the power of each fit index to detect the model misspecification.
plotPowerFit(simtwotoone200, cutoff=hubentler99, alpha=0.05,
usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
Interpretation:
- RMSEA: Many replications show RMSEA > 0.05, but some are 0.
- CFI: Most replications have CFI > 0.95, indicating low power.
- TLI: Similar to CFI but slightly more sensitive.
- SRMR: Nearly all replications have SRMR < 0.06, providing almost no power.
Step 6: Quantify Power of Fit Indices
getPowerFit(simtwotoone200, cutoff=hubentler99, alpha=0.05,
usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))## cfi tli rmsea srmr
## 0.037 0.182 0.523 0.000The highest power belongs to RMSEA (\(\approx\) 0.523), while the others show power below 0.20 — a clear indication that these cutoffs perform poorly for this misspecified model.
Step 7: Create Dynamic Fit Cutoffs
Dynamic fit indices (DFI) adapt the cutoff values to specific models and sample sizes. We first generate empirical population data (matching the model-implied parameters).
dattemp2 <- generate(twofactormodel, n=200, empirical=TRUE)To confirm that this empirical data reproduces the population parameters:
twofactormodelforanalysis <- "
f1 =~ x1 + x2 + x3
f2 =~ x4 + x5 + x6
"
fittemptrue <- cfa(twofactormodelforanalysis, data=dattemp2)
summary(fittemptrue, std=TRUE)## lavaan 0.6-20 ended normally after 17 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 13
##
## Number of observations 200
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 8
## P-value (Chi-square) 1.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## f1 =~
## x1 1.000 0.700 0.700
## x2 1.000 0.121 8.270 0.000 0.700 0.700
## x3 1.000 0.121 8.270 0.000 0.700 0.700
## f2 =~
## x4 1.000 0.700 0.700
## x5 1.000 0.121 8.270 0.000 0.700 0.700
## x6 1.000 0.121 8.270 0.000 0.700 0.700
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## f1 ~~
## f2 0.441 0.071 6.208 0.000 0.900 0.900
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.510 0.066 7.703 0.000 0.510 0.510
## .x2 0.510 0.066 7.703 0.000 0.510 0.510
## .x3 0.510 0.066 7.703 0.000 0.510 0.510
## .x4 0.510 0.066 7.703 0.000 0.510 0.510
## .x5 0.510 0.066 7.703 0.000 0.510 0.510
## .x6 0.510 0.066 7.703 0.000 0.510 0.510
## f1 0.490 0.096 5.114 0.000 1.000 1.000
## f2 0.490 0.096 5.114 0.000 1.000 1.000The standardized estimates perfectly match the data-generating values.
Step 8: Fit the Misspecified Model and Run DFI
We now analyze the empirical data with the one-factor model.
fittemp2 <- cfa(onefactormodel, data=dattemp2)
summary(fittemp2, std=TRUE)## lavaan 0.6-20 ended normally after 19 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 12
##
## Number of observations 200
##
## Model Test User Model:
##
## Test statistic 5.392
## Degrees of freedom 9
## P-value (Chi-square) 0.799
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## f1 =~
## x1 1.000 0.679 0.679
## x2 1.000 0.123 8.127 0.000 0.679 0.679
## x3 1.000 0.123 8.127 0.000 0.679 0.679
## x4 1.000 0.123 8.127 0.000 0.679 0.679
## x5 1.000 0.123 8.127 0.000 0.679 0.679
## x6 1.000 0.123 8.127 0.000 0.679 0.679
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.539 0.065 8.305 0.000 0.539 0.539
## .x2 0.539 0.065 8.305 0.000 0.539 0.539
## .x3 0.539 0.065 8.305 0.000 0.539 0.539
## .x4 0.539 0.065 8.305 0.000 0.539 0.539
## .x5 0.539 0.065 8.305 0.000 0.539 0.539
## .x6 0.539 0.065 8.305 0.000 0.539 0.539
## f1 0.461 0.092 5.026 0.000 1.000 1.000Then, run the Dynamic Fit Index procedure:
library(dynamic)
outdynamic <- cfaOne(fittemp2, plot=TRUE)
outdynamic## Your DFI cutoffs:
## SRMR RMSEA CFI
## Level-0 0.034 0.066 0.979
## Specificity 95% 95% 95%
##
## Level-1 0.034 0.066 0.979
## Sensitivity 59% 66% 62%
##
## Level-2 0.035 0.072 0.978
## Sensitivity 92% 95% 95%
##
## Level-3 0.043 0.112 0.951
## Sensitivity 95% 95% 95%
##
## Empirical fit indices:
## Chi-Square df p-value SRMR RMSEA CFI
## 5.392 9 0.799 0.02 0 1
##
## The distributions for each level are in the Plots tab
## [[1]]
##
## [[2]]
##
## [[3]]
Extract the Level 1–3 cutoff values:
level1 <- outdynamic$cutoffs["Level-1",]
level2 <- outdynamic$cutoffs["Level-2",]
level3 <- outdynamic$cutoffs["Level-3",]Step 9: Evaluate Power Using DFI Cutoffs
getPowerFit(simtwotoone200, cutoff=level1, alpha=0.05,
usedFit=c("RMSEA", "CFI", "SRMR"))## cfi rmsea srmr
## 0.335 0.335 0.373getPowerFit(simtwotoone200, cutoff=level2, alpha=0.05,
usedFit=c("RMSEA", "CFI", "SRMR"))## cfi rmsea srmr
## 0.315 0.263 0.328getPowerFit(simtwotoone200, cutoff=level3, alpha=0.05,
usedFit=c("RMSEA", "CFI", "SRMR"))## cfi rmsea srmr
## 0.043 0.015 0.091As expected, Level-1 cutoffs yield the highest power (\(\approx\) 0.335–0.421), though still quite low. The dynamic approach is more sensitive than fixed cutoffs but still struggles to reject the misspecified model.
Step 10: Increase the Sample Size
Let’s check whether a larger sample size (n = 800) increases power.
simtwotoone800 <- sim(nRep=1000, model=onefactormodel, generate=twofactormodel, n=800, lavaanfun="cfa")Generate empirical data and obtain new DFI cutoffs:
dattemp3 <- generate(twofactormodel, n=800, empirical=TRUE)
fittemp3 <- cfa(onefactormodel, data=dattemp3)
outdynamic3 <- cfaOne(fittemp3)
outdynamic3## Your DFI cutoffs:
## SRMR RMSEA CFI
## Level-0 0.017 0.032 0.995
## Specificity 95% 95% 95%
##
## Level-1 0.023 0.059 0.985
## Sensitivity 95% 95% 95%
##
## Level-2 0.036 0.103 0.956
## Sensitivity 95% 95% 95%
##
## Level-3 0.049 0.142 0.923
## Sensitivity 95% 95% 95%
##
## Empirical fit indices:
## Chi-Square df p-value SRMR RMSEA CFI
## 21.568 9 0.01 0.02 0.042 0.992Then evaluate power again:
level1_3 <- outdynamic3$cutoffs["Level-1",]
level2_3 <- outdynamic3$cutoffs["Level-2",]
level3_3 <- outdynamic3$cutoffs["Level-3",]
getPowerFit(simtwotoone800, cutoff=level1_3, alpha=0.05,
usedFit=c("RMSEA", "CFI", "SRMR"))## cfi rmsea srmr
## 0.422 0.345 0.530getPowerFit(simtwotoone800, cutoff=level2_3, alpha=0.05,
usedFit=c("RMSEA", "CFI", "SRMR"))## cfi rmsea srmr
## 0.000 0.000 0.002getPowerFit(simtwotoone800, cutoff=level3_3, alpha=0.05,
usedFit=c("RMSEA", "CFI", "SRMR"))## cfi rmsea srmr
## 0 0 0With a larger sample, Level-1 power slightly improves (\(\approx\) 0.345–0.530). However, Level-2 and Level-3 cutoffs yield zero power, meaning no replications were rejected.
Step 11: Visualize Sampling Distributions
To understand why power decreases, let’s visualize the sampling variability of fit indices:
plotPowerFit(simtwotoone800, cutoff=hubentler99, alpha=0.05,
usedFit=c("RMSEA", "CFI", "TLI", "SRMR"))
The x-axis range becomes narrower at larger samples, showing reduced variability. As the standard errors of fit indices decrease, the probability of obtaining extreme (bad) fit values declines, thus reducing the power to reject the model.
Step 12: Interpretation and Key Takeaways
The relationship between sample size and power to reject misspecified models is not straightforward:
- Larger n decreases sampling variability in fit indices.
- If the misspecified model’s population fit is better than the cutoff, power decreases with larger n.
- If it’s worse than the cutoff, power increases with larger n.
Therefore, always evaluate the power of your fit indices (either traditional or dynamic). If your indices show insufficient power to reject a misspecified model, consider using alternative methods or additional evidence to assess model adequacy.
Power Analysis in Test of Close Fit and Test of Not Close Fit
According to Browne & Cudeck (1993), researchers can statistically test whether the obtained RMSEA value is greater than or less than a given threshold (commonly 0.05). Two types of tests are defined:
- Test of Close Fit: \(H_0: \varepsilon \le 0.05\) vs. \(H_A: \varepsilon > 0.05\). Researchers aim to retain the null hypothesis, concluding that the observed fit indices could come from a population model with RMSEA = 0.05.
- Test of Not Close Fit: \(H_0: \varepsilon \ge 0.05\) vs. \(H_A: \varepsilon < 0.05\). Researchers aim to reject the null hypothesis, claiming that the population model fits more closely than RMSEA = 0.05.
It is important to emphasize that both RMSEA = 0.05 and RMSEA = 0.08 are arbitrary thresholds historically used to distinguish “good” from “poor” model fit. As demonstrated in Example 1, even a perfectly fitting model can yield many replications with RMSEA values exceeding 0.05 (> 10%). Nonetheless, because these tests continue to be discussed in the literature, this section illustrates how to evaluate their statistical power.
Step 1: Examine the Sampling Distribution of RMSEA
Suppose the population RMSEA is 0.05, with a sample size of 200 and
20 degrees of freedom. We can visualize its sampling distribution using
the plotRMSEAdist() function from the semTools
package.
library(semTools)
plotRMSEAdist(0.05, n = 200, df = 20)
The distribution shows the expected variability of RMSEA estimates around the population value.
Step 2: Define the Critical Region for the Test of Close Fit
For the test of close fit, the rejection region lies on the right tail (i.e., larger RMSEA values are considered evidence against close fit).
plotRMSEAdist(0.05, n = 200, df = 20, ptile = .95)
Step 3: Define the Critical Region for the Test of Not Close Fit
Conversely, for the test of not close fit, the rejection region lies on the left tail of the distribution (i.e., smaller RMSEA values are considered evidence of close fit). To make the plot clearer, let’s increase the sample size to 1000, because at n = 200, many simulated RMSEA values are close to zero.
plotRMSEAdist(0.05, n = 1000, df = 20, ptile = .05)
Step 4: Compute the Power of the Test of Close Fit
To evaluate the power of these tests, we use the
findRMSEApower() function in semTools. The
method follows the approach proposed by MacCallum
et al. (1996).
Let the null hypothesis specify a population RMSEA of 0.05, and the true population RMSEA be 0.08. This corresponds to a test of close fit, where we wish to reject \(H_0: \varepsilon \le 0.05\).
findRMSEApower(rmsea0 = 0.05, rmseaA = 0.08, df = 20, n = 200)## [1] 0.454368The resulting power is 0.45, indicating that with n = 200 and df = 20, there is only a 45% chance of correctly rejecting the null hypothesis when the population RMSEA truly equals 0.08.
Step 5: Compute the Power of the Test of Not Close Fit
Now let’s reverse the logic to test not close fit. Suppose the null hypothesis specifies RMSEA = 0.05, while the true population RMSEA is 0.04. We wish to reject \(H_0: \varepsilon \ge 0.05\), indicating evidence of a close-fitting model.
findRMSEApower(rmsea0 = 0.05, rmseaA = 0.04, df = 25, n = 500)## [1] 0.2265141Here, the power is 0.23, which is quite low. In other words, even with 500 participants and 25 degrees of freedom, there is only a 23% chance of correctly identifying that the population RMSEA is truly less than 0.05.
Let’s increase the sample size to 2000:
findRMSEApower(rmsea0 = 0.05, rmseaA = 0.04, df = 25, n = 2000)## [1] 0.6784104Now, the power increases to 0.68 (68%), showing that larger samples
help but still do not guarantee adequate power. Please note that in the
findRMSEApower() function,
- If
rmseaA>rmsea0, the test of close fit is performed. - If
rmseaA<rmsea0, the test of not close fit is performed.
This distinction determines which tail of the distribution defines the rejection region.
Step 6: Determine Required Sample Size
You can also determine the required sample size to achieve a desired
power level (e.g., 0.80) using findRMSEAsamplesize().
findRMSEAsamplesize(rmsea0 = 0.05, rmseaA = 0.04, df = 25, power = 0.8)## [1] 2699The result indicates that approximately 2,699 cases are needed to achieve 80% power for detecting a true RMSEA of 0.04 when testing against 0.05 with 25 degrees of freedom.
Interpretation and Caution
As these examples demonstrate, the test of close fit and test of not close fit are highly sensitive to sample size and degrees of freedom. However, because RMSEA thresholds such as 0.05 or 0.08 are arbitrary and model-dependent, relying solely on these tests can be misleading.
Therefore, while this section illustrates how such tests work and how to evaluate their power, I do not recommend using fixed RMSEA cutoffs or their associated tests for practical model evaluation. Researchers should instead consider multiple fit indices, model complexity, theoretical grounding, and simulation-based evidence when assessing model fit.
Accuracy in Parameter Estimation in Parameters
In recent years, researchers have been increasingly encouraged to report confidence intervals (CIs) for parameter estimates. The width of a CI reflects the accuracy of parameter estimation—narrower intervals indicate greater precision. However, achieving a narrow interval often requires a larger sample size, which may not always be feasible or justifiable given limited resources. Thus, researchers must balance desired precision with practical considerations such as cost and effort.
Let us revisit the Cohen’s d example (Example 3 above). The data-generating model, analysis model, and simulation setup can be specified as follows. Please refer back to Example 3 for details on the setup.
gencohend <- "
f1 =~ 0.8*x1 + 0.7*x2 + 0.6*x3 + 0.5*x4
f2 =~ 0.8*x5 + 0.7*x6 + 0.6*x7 + 0.5*x8
f2 ~~ 0.657*f1
f1 ~~ 1*f1
f2 ~~ 1.2*f2
x1 ~~ 0.108*x5
x2 ~~ 0.153*x6
x3 ~~ 0.192*x7
x4 ~~ 0.225*x8
x1 ~~ 0.36*x1
x2 ~~ 0.51*x2
x3 ~~ 0.64*x3
x4 ~~ 0.75*x4
x5 ~~ 0.36*x5
x6 ~~ 0.51*x6
x7 ~~ 0.64*x7
x8 ~~ 0.75*x8
f1 ~ 0*1
f2 ~ 0.5*1
x1 ~ 3.3*1
x2 ~ 2.5*1
x3 ~ 4.2*1
x4 ~ 3.7*1
x5 ~ 3.3*1
x6 ~ 2.5*1
x7 ~ 4.2*1
x8 ~ 3.7*1
"
cohendmodel <- "
f1 =~ 1*x1 + L2*x2 + L3*x3 + L4*x4
f2 =~ 1*x5 + L2*x6 + L3*x7 + L4*x8
f1 ~~ NA*f2
f1 ~~ NA*f1 + V1*f1
f2 ~~ NA*f2 + V2*f2
x1 ~~ x5
x2 ~~ x6
x3 ~~ x7
x4 ~~ x8
x1 ~ 0*1
x2 ~ I2*1
x3 ~ I3*1
x4 ~ I4*1
x5 ~ 0*1
x6 ~ I2*1
x7 ~ I3*1
x8 ~ I4*1
f1 ~ M1*1
f2 ~ M2*1
eff := (M2 - M1)/sqrt(V1)
"
simcohend200 <- sim(nRep=1000, model=cohendmodel, generate=gencohend, n=200, lavaanfun="cfa")Suppose an assessment company wishes to estimate Cohen’s d with a 95%
confidence interval whose total width is less than 0.20 (corresponding
to a margin of error of \(\pm 0.10\)).
The getCIwidth() function can be used to evaluate this
width:
getCIwidth(simcohend200)## L2 <- (f1=~x2) L3 <- (f1=~x3) L4 <- (f1=~x4) L2 <- (f2=~x6) L3 <- (f2=~x7)
## 50% 0.2815853 0.269113 0.264643 0.2815853 0.269113
## L4 <- (f2=~x8) f1~~f2 V1 V2 x1~~x5 x2~~x6 x3~~x7
## 50% 0.264643 0.3029075 0.3770516 0.4377194 0.1841124 0.1929652 0.2116934
## x4~~x8 I2 <- (x2~1) I3 <- (x3~1) I4 <- (x4~1) I2 <- (x6~1) I3 <- (x7~1)
## 50% 0.2320527 1.005879 0.9624364 0.9503902 1.005879 0.9624364
## I4 <- (x8~1) M1 M2 x1~~x1 x2~~x2 x3~~x3 x4~~x4
## 50% 0.9503902 0.267946 0.2852692 0.2558156 0.2642428 0.2880625 0.31552
## x5~~x5 x6~~x6 x7~~x7 x8~~x8 eff [ L2 ] - [ L2 ]
## 50% 0.2678335 0.2721939 0.2926336 0.3179406 0.3373998 0
## [ L3 ] - [ L3 ] [ L4 ] - [ L4 ] [ I2 ] - [ I2 ] [ I3 ] - [ I3 ]
## 50% 0 0 0 0
## [ I4 ] - [ I4 ]
## 50% 0This function returns the median (50th percentile) of the CI width
obtained across all replications. For the parameter eff,
the median 95% CI width is 0.337. Researchers may also be interested in
the 80th percentile of the CI width distribution, which represents the
width achieved in 80% of the replications:
getCIwidth(simcohend200, assurance = 0.8)## L2 <- (f1=~x2) L3 <- (f1=~x3) L4 <- (f1=~x4) L2 <- (f2=~x6) L3 <- (f2=~x7)
## 80% 0.3079244 0.2908378 0.2858425 0.3079244 0.2908378
## L4 <- (f2=~x8) f1~~f2 V1 V2 x1~~x5 x2~~x6 x3~~x7
## 80% 0.2858425 0.3326791 0.4135065 0.4789868 0.1957453 0.2059657 0.2258218
## x4~~x8 I2 <- (x2~1) I3 <- (x3~1) I4 <- (x4~1) I2 <- (x6~1) I3 <- (x7~1)
## 80% 0.248309 1.095954 1.041218 1.022381 1.095954 1.041218
## I4 <- (x8~1) M1 M2 x1~~x1 x2~~x2 x3~~x3 x4~~x4
## 80% 1.022381 0.2795544 0.2978972 0.2748137 0.2840672 0.3124003 0.3432419
## x5~~x5 x6~~x6 x7~~x7 x8~~x8 eff [ L2 ] - [ L2 ]
## 80% 0.2876932 0.2920166 0.3153916 0.3436684 0.3615181 0
## [ L3 ] - [ L3 ] [ L4 ] - [ L4 ] [ I2 ] - [ I2 ] [ I3 ] - [ I3 ]
## 80% 0 0 0 0
## [ I4 ] - [ I4 ]
## 80% 0The degree of assurance (Lai & Kelley, 2011) refers to the proportion of replications whose CI width is below a specified threshold. In this case, an assurance level of 0.80 means that 80% of the replications yield CI widths smaller than a given value. Here, 80% of the replications have widths below 0.362. However, our goal is to ensure that the width is below 0.20 with 80% assurance.
To achieve this, we can conduct simulations with varying sample sizes, ranging from 200 to 600:
simcohendvary <- sim(model=cohendmodel, generate=gencohend, n=200:600, lavaanfun="cfa")We can then plot the CI width (with 80% assurance) as a function of sample size:
plotCIwidth(simcohendvary, "eff", assurance = 0.80)
abline(h = 0.2, col="blue")
A horizontal blue line representing the target width (0.20) is added
to the plot for reference. From the plot, it appears that the desired
precision (width \(\le\) 0.20) is
achieved when the sample size is approximately 570. We can confirm this
by extracting the estimated CI widths for sample sizes of 570 and 600
using getCIwidth():
getCIwidth(simcohendvary, assurance = 0.80, nVal=570)["eff"]## eff
## 1 0.2010979getCIwidth(simcohendvary, assurance = 0.80, nVal=600)["eff"]## eff
## 1 0.1904233The estimated widths with 80% assurance are 0.201 for n = 570 and 0.190 for n = 600. Note that these are simulation-based estimates and may vary slightly across replications. To obtain a more precise estimate, we can run a focused simulation with n = 600 and recheck the 80th percentile of the CI width:
simcohend600 <- sim(nRep=1000, model=cohendmodel, generate=gencohend, n=600, lavaanfun="cfa")getCIwidth(simcohend600, assurance = 0.8)## L2 <- (f1=~x2) L3 <- (f1=~x3) L4 <- (f1=~x4) L2 <- (f2=~x6) L3 <- (f2=~x7)
## 80% 0.1719212 0.1631439 0.1599227 0.1719212 0.1631439
## L4 <- (f2=~x8) f1~~f2 V1 V2 x1~~x5 x2~~x6 x3~~x7
## 80% 0.1599227 0.1863229 0.2309709 0.267515 0.1110514 0.1159172 0.127961
## x4~~x8 I2 <- (x2~1) I3 <- (x3~1) I4 <- (x4~1) I2 <- (x6~1) I3 <- (x7~1)
## 80% 0.1408465 0.6130051 0.5843383 0.5731134 0.6130051 0.5843383
## I4 <- (x8~1) M1 M2 x1~~x1 x2~~x2 x3~~x3 x4~~x4
## 80% 0.5731134 0.1591607 0.1695457 0.1547499 0.1600343 0.176596 0.1930036
## x5~~x5 x6~~x6 x7~~x7 x8~~x8 eff [ L2 ] - [ L2 ]
## 80% 0.1619992 0.1643044 0.1770307 0.193613 0.20281 0
## [ L3 ] - [ L3 ] [ L4 ] - [ L4 ] [ I2 ] - [ I2 ] [ I3 ] - [ I3 ]
## 80% 0 0 0 0
## [ I4 ] - [ I4 ]
## 80% 0The resulting CI width (with 80% assurance) is 0.203. Therefore, the assessment company should plan to collect data from approximately 600 trials to ensure that the confidence interval for Cohen’s d is no wider than 0.203 with an 80% guarantee.
It is worth noting that the earlier estimate from the varying sample size simulation predicted a width of 0.190, which slightly overestimated the precision (i.e., underestimated the true width). In some situations, the predictive model could instead underestimate the accuracy. For this reason, it is generally advisable to rerun the full simulation—as demonstrated in this example—to obtain more reliable and precise results. In this case, maintaining a sample size of 600 trials is a prudent choice, as it yields a CI width very close to the desired target of 0.20.
Accuracy in RMSEA Estimation
Researchers may also wish to determine the sample size required to achieve a desired width of the RMSEA confidence interval (Kelley & Lai, 2011). This approach focuses on the precision of RMSEA estimation rather than on model fit testing.
First, let’s compute the confidence interval of RMSEA for a given
observed value (assuming multivariate normality). We can use the
ci.rmsea() function from the MBESS
package:
library(MBESS)
ci.rmsea(rmsea=0.05, df=20, N=500, conf.level=.90)## $Lower.Conf.Limit
## [1] 0.03043259
##
## $RMSEA
## [1] 0.05
##
## $Upper.Conf.Limit
## [1] 0.06962194For df = 20 and n = 500, the 90% confidence interval for RMSEA ranges from 0.03 to 0.07, yielding a total width of 0.04.
Suppose we want to find the sample size required so that the width of
the 90% confidence interval is no greater than 0.03, with the same
degrees of freedom (df = 20) and an expected RMSEA of 0.05. We
can use the ss.aipe.rmsea() function for this purpose:
ss.aipe.rmsea(RMSEA=0.05, df=20, width=0.03, conf.level=.90)## Necessary sample size so that the expected width of the 90% confidence interval
## is no greater than 0.03, given a population RMSEA of 0.05, is:## [1] 760The output shows that 760 participants are needed to achieve a 90% confidence interval width of less than 0.03 for RMSEA.
Note that a 90% confidence level is typically used (rather than 95%) because it corresponds to the one-tailed nature of the test of close fit and test of not close fit procedures.