# Determining sample size for Structural Equation Model with R

I want to run and experiment to test the validity of a theoretical model. I have the structural and measurement model, but do not have previous data on any of the parameters. Is there a way I can find the ideal sample size?

I use R, and have read that simsem package can be useful. However, I could not find examples that demonstrate how to proceed. I request help in pointing me to some references that can help me work on this? While I want to continue using R, if there is any other package besides simsem that can work, I am open to exploring that also. Please suggest.

• If you don't have any idea about the parameters, make a series of sensible guesses. Your bigger problem is convergence when your sample is too small. – Jeremy Miles Oct 6 '19 at 4:07

I think the best walkthrough of how to sample size plan using lavaan and simsem I have ever read was in a very short and very accessible text called Latent Variable Modeling Using R: A Step-by-Step Guide by A. Alexander Beaujean. His examples are admittedly relatively simple, but they should provide an excellent basis to start from when thinking about a defensible sample size planning strategy.

As a toy example, I can share a bit of code I developed as part of a briefly formed working group in my grad program. We were trying to develop (or at least starting to sketch out) a simulation-based statistical power tool capable of estimating power and sample size for a wide range of model and data types.

## Code Walkthrough

Start by establishing a working directory and loading the necessary libraries:

wd<-'path/to/folder/named/SimulateData Example'
setwd(wd)

out.folder<-paste0(wd, '/Output')
dir.create(out.folder)

#libraries needed for main function.
library(lavaan)
library(MASS)
library(simsem)


Then specify the population model. Note that the population model represents your prior beliefs (and they say frequentist statisticians don't incorporate priors!) about the multivariate distribution of the scores in your target population. That is a confusing mouthful, and it is admittedly difficult. Had I not been taught the mystical art of path tracing in Greg Hancock's SEM course, I might have found it much more frustrating. A decent summary of the algebraic tracing rules can be found here.

#specify model
pop.model<-'
Y~.2*X1+.5*X2
X1~~.4*X2
Y~~.63*Y
'


There is no doubt about it, specification of the population model is the most demanding aspect of this approach, both conceptually and practically. Just remember not to put too much pressure on yourself to come up with the perfect power analysis. Instead, strive for a defensible power analysis. There is no avoiding a challenge when working on models with complex latent variable and/or autocorrelation structures. Still, specifying a reasonable (and mathematically tractable) population mode can be a pain, but it is doable.

Next you specify the test.model, which is the "sample" lavaan-syntax model, just without the fixed paths.

test.model<-'
Y~X1+X2
X1~~X2
Y~~Y
'


I should quickly pause here and say that I am also making some other assumptions about my population and the the measures in this toy problem. For instance, the approach below assumes the scores are derived from a multivariate normal distribution. More subtly, the model also implies a directionality in its effects. I do not mean to imply causality, but on some level a choice was made to regress Y on X1 and X2. Why is this pattern of effects more defensible than alternative specifications? You should arguably have a reasonable answer. Okay back to the code.

Below is a quick example of how the simulateData function can generate multiple samples.

Note that by passing the vector c(10, 15, 20) to the sample.nobs argument, I am able to generate three samples that all get appended to the same data.frame called demo. The new variable group distinguishes between the three randomly generated (and in this case differently sized) data sets, which you will see if you run this code on your own machine.

demo<-simulateData(model = pop.model, sample.nobs = c(10, 15, 20))
str(demo)


And just so you can get a sense of how the new data.frame gets structured you can run:

head(demo); tail(demo)


Moving on to my toy program, which I still think if built out in an intelligent, modular way could really work as an incredibly useful power analysis tool. I just keep not having the time to really get around to it.

As you'll see, the function can accept skewness and kurtosis arguments, thus allowing some flexibility about even the need to strictly assume a multivariate normal distribution when one doesn't think that is reasonable expectation in the population (think of psychopathology symptom scores in a non-clinical population). I'll add a little annotation throughout to hopefully explain the basics.

simPwr<-function(n.max=500, n.min=20, n.sims=1000, skew.val=0,
kurtosis.val=0, pop.model, test.model, alpha=.05){

# Here is my current inelegant workaround for making this more flexible - i.e., incoporate more predictors
# Basically I fit a model using data generated from the population model just so I can get the number of regression paths
# There has got to be an easier way people...
tmp.data <- simulateData(model = pop.model, sample.nobs = 500)
tmp.fit <- sem(test.model, data = tmp.data)
tmp.par.table <- parameterEstimates(tmp.fit)
n.reg.paths <- length(tmp.par.table$$op[tmp.par.table$$op == '~'])

# Create a sequence of sample sizes from n.min to n.max by 10
n.obs<-seq(n.min, n.max, by=10)

# data.frame to retain final results
pow <- data.frame()

# You could imagine all sorts of ways to branch out these assumptions and make a model more flexible
# Right now the function is only set to monitor statistical power of "regression paths"
# You could imagine passing a character vector of lavaan-operator-defined paths you want to analyze
# You could structure the program to accept a correlation or covariance matrix as that then gets converted to lavaan syntax (may be more intuitive for users)
# Add some decent error handling and there could be something here I swear... (reflecting mainly for self-motivation purposes)

# Now loop through the the vector of different sample size
for(n in 1:length(n.obs)){
# And draw s samples of sample size n using the population model
# Also this is not most elegant version of this program but should still illustrate the process
est <- data.frame()
pval <- data.frame()

# Definitely some expensive for loops here. Code in need of optimization.
for(s in 1:n.sims){
tmp.data<-simulateData(pop.model, sample.nobs = n.obs[n], skewness = skew.val, kurtosis = kurtosis.val)
tmp.fit<-sem(test.model, data=tmp.data, fixed.x = F)
parm<-parameterestimates(tmp.fit)
parm<-parm[parm$$op=='~',] est<-rbind(est, parm$$est)
pval<-rbind(pval, parm$pvalue) } Power<-ifelse(pval<alpha, 1, 0) # Save values in a data.frame (mainly for graphing) and final analysis mu.est<-colMeans(est) mu.pow<-colMeans(Power) sim.N<-rep(n.obs[n], n.reg.paths) parm.nm <- paste(parm$$lhs[parm$$op=="~"], parm$$op[parm$$op=="~"], parm$$rhs[parm$$op=="~"]) dat.fin<-data.frame(parm.nm, mu.est, mu.pow, sim.N) colnames(dat.fin)<-c('Parameter', 'Estimate', 'Power', 'N') pow<-rbind(pow, dat.fin) } return(pow) } test<-simPwr(n.max=350, n.min = 10, n.sims = 1000, pop.model = pop.model, test.model = test.model)  So you see that I have run the mini program to obtain power estimates for samples ranging from 10 (unrealisitically low) up to 350 (a modest $$N$$ for an SEM model). Now using loess regression combined with a little upsampling, I can see the smallest sample size that yields a power of at least .80. First for Y~X1... N350<-1:350 fit.lo1<-loess(Power~N, data = test[test$$Parameter=='Y ~ X1',]) ROW<-predict(fit.lo1, newdata=N350) TEST<-cbind(ROW, N350) TEST<-TEST[TEST[,1]>.8,] TEST<-as.data.frame(TEST) TEST<-na.omit(TEST) TEST$$N350[1]  Which results in a sample size of about 150. Note this is not the most robust way to obtain this value, but it works in this case (and we'll plot the results to be sure). Now for Y~X2 fit.lo2<-loess(Power~N, data = test[test$Parameter=='Y ~ X2',])
ROW<-predict(fit.lo2, newdata=N350)
TEST<-cbind(ROW, N350)
TEST<-TEST[TEST[,1]>.8,]
TEST<-as.data.frame(TEST)
TEST<-na.omit(TEST)
TEST\$N350[1]


Which results in a sample size of about 30.

And now to plot the data just to be sure it makes sense:

library(ggplot2)
g1<-ggplot(aes(x=N, y=Power, group = Parameter, color=Parameter), data=test)
g2<-g1+stat_smooth(se=F)
g3<-g2+ylab('Power')+xlab('Sample Size')
g4<-g3+geom_hline(yintercept = .8, color='red', lty='dashed')
g5<-g4+ggtitle('Power Curves for Parameters in Model')
g5


As you can see, an artifact of using the default loess model is that one of the power curves briefly jumps above 1, which doesn't make any sense. So there are lots of opportunities to improve this workflow, but I think it offers a reasonable way to think about power analyses in SEM models.

Beaujean, A. A. (2014). Latent Variable Modeling Using R: A Step-by-Step Guide. New York: Routledge.