MGCFA Power simulations lavaan

[Camille Williams]

Hi, 

I am trying to do a power simulation to evaluate the number of participants I need to have in each of my two groups to observe metric and scalar invariance between groups at 80% power or a power sensitivity simulation to identify the minimum effect size for which I can observe a group difference in scalar or metric invariance with 80% power. 

 

If there is metric and/or scalar invariance between my groups, I want to calculate the sample required to observe a group difference in each intercept or factor loading. 

Any help and papers welcome! 

Best,

Camille

[Camille Williams]

I know that some papers recommend 200 participants while others say that 50 is enough if certain requirements are met (https://www.tandfonline.com/doi/abs/10.1207/s15327574ijt0502_4). 

However, I want to know the effect size such a sample can identify given my model.

[Terrence Jorgensen]

Any help and papers welcome! 

Monte Carlo simulations are used for power analysis in SEM:

https://doi.org/10.1207/S15328007SEM0904_8

https://doi.org/10.1207/S15328007SEM0802_7

https://doi.org/10.1177/0165025413515169

However, I want to know the effect size such a sample can identify given my model.

You can manipulate effect size the same way you manipulate N in a Monte Carlo study.  The simsem package facilitates such a study using lavaan, although the lavaanList() functionality might be sufficient for your power analysis.  There are vignette's to help with simsem, some specifically about power.

https://github.com/simsem/simsem/wiki/Vignette

Terrence D. Jorgensen

Assistant Professor, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[Camille Williams]

Thank you for these papers! 

To clarify, I am interested in group differences in intercepts and loadings.  

I want to examine the effect size of the group difference in factor loadings and intercepts that I can identify given my sample size and model. 

[Camille Williams]

Hi Terrence, 

Unfortunately, I can't seem to figure it out... I am interested in group differences in intercepts and loadings. I want to examine the effect size of the group difference in factor loadings and intercepts that I can identify given my sample size and model.

A) I don't understand how to use the LavaanList functionality with multiple groups based on the package info (https://cran.r-project.org/web/packages/lavaan/lavaan.pdf). 

B) I am currently trying to understand how to manipulate the simsem package: https://github.com/simsem/simsem/wiki/Example-6:-Multiple-Groups to establish the power of detecting a difference in intercepts or slopes between my groups based on my model.

My model is one factor explained by 22 observed variables, with several residual covariances, factor variance set to 1. 

TBV.model <- 'TBV =~ NA*CorticalWhiteMatterVol_log

+ Brain_Stem_log

+ Left_VentralDC_log + Right_VentralDC_log

+ lhCortexVol_log + rhCortexVol_log

+ Left_Accumbens_area_log + Right_Accumbens_area_log

+ Left_Cerebellum_Cortex_log + Right_Cerebellum_Cortex_log

+ Left_Caudate_log + Right_Caudate_log

+ Left_Hippocampus_log + Right_Hippocampus_log

+ Left_Pallidum_log + Right_Pallidum_log

+ Left_Putamen_log + Right_Putamen_log

+ Right_Thalamus_Proper_log + Left_Thalamus_Proper_log

+ Left_Amygdala_log + Right_Amygdala_log

TBV ~~ 1*TBV

TBV ~ c(0, 0.053)* 1

lhCortexVol_log ~~      rhCortexVol_log

Left_Caudate_log ~~      Right_Caudate_log

Left_Putamen_log       ~~      Right_Putamen_log

Left_Cerebellum_Cortex_log     ~~      Right_Cerebellum_Cortex_log

Right_Thalamus_Proper_log    ~~      Left_Thalamus_Proper_log

Left_Accumbens_area_log ~~      Right_Accumbens_area_log

Left_VentralDC_log      ~~      Right_VentralDC_log

Left_Amygdala_log    ~~      Right_Amygdala_log

Left_Hippocampus_log  ~~      Right_Hippocampus_log

Left_Pallidum_log ~~ Right_Pallidum_log

Right_Putamen_log       ~~      Right_Amygdala_log

Right_Caudate_log     ~~      Right_Putamen_log

Left_Caudate_log~~Left_Thalamus_Proper_log'

1) Since I want to know the power of my model to detect group differences in loadings and intercepts and also the minimum effect size of the group difference in loading or intercept that can be detected by my model, I should not free my parameters? 

2) How do I obtain the VTE (or TE) and RPS (or PS) for my current model in each group?

3) So far I have loadings, intercepts and mean factor which I know where to put in the model function, but how do I include the residual correlations and observed variable variance which differs between groups? (Code for all underneath) 

a. my loadings for each group  

loadingVal_Control <-as.matrix(parameterEstimates(Control_fit_TBV1, standardized=TRUE) %>%

                                 filter(op == "=~") %>%

                                 select(Beta=std.all))

loadingVal_ASD <-as.matrix(parameterEstimates(ASD_fit_TBV1, standardized=TRUE) %>%

                                 filter(op == "=~") %>%

                                 select(Beta=std.all))

LY_not_free = list(loadingVal_Control, loadingVal_ASD)

b. my intercepts for each group

TYControl_all <-as.matrix(parameterEstimates(Configural_Model, standardized=TRUE) %>% 

                      select(Beta=std.all))

TY1 <- TYControl_all[196:217,]

TYASD_all <-as.matrix(parameterEstimates(Configural_Model, standardized=TRUE) %>% 

                      select(Beta=std.all))

TY2 <- TYASD_all[196:217,]

c. my factor mean for each group

AL1 <-c(0)

AL2 <- c(0.053)

AL_not_free= list(AL1, AL2)

d. correlation matrix for each group 

Corr_ASD <- as.matrix(parameterEstimates(ASD_fit_TBV1, standardized=TRUE) %>% 

                      filter(op == "~~") %>% 

                      select(Beta=std.all))

Corr_Control <- as.matrix(parameterEstimates(Control_fit_TBV1, standardized=TRUE) %>% 

                     filter(op == "~~") %>% 

                     select(Beta=std.all))

e. my variance for each group

Var_Control_all <- as.matrix(parameterEstimates(Control_fit_TBV1, standardized=TRUE) %>% 

                        select(Beta=std.all))

Var_Control <- Var_Control_all[173:194,]

Var_ASD_all <- as.matrix(parameterEstimates(ASD_fit_TBV1, standardized=TRUE) %>% 

                        select(Beta=std.all))

Var_ASD <- Var_ASD_all[173:194,]

Thank you for your time and your advice! I dont know where to turn! 

Camille 

[Camille Williams]

Hi Terrence,

How can i use the lavaanList() function to examine the effect size of the group difference in factor loadings and intercepts that I can identify given my sample size and model. 

I can't seem to understand how to do this with the simsem package or the lavaanList() option. 

Thank you,

Camille

On Friday, March 29, 2019 at 11:02:56 AM UTC+1, Terrence Jorgensen wrote:

[Terrence Jorgensen]

 can't seem to understand how to do this with the simsem package or the lavaanList() option. 

Then don't use them. Here is a toy example you can adapt for your own study.  If you don't understand it, then you need to hire a consultant.  I suggest asking on SEMNET:  http://www2.gsu.edu/~mkteer/semnet.html

## write functions for simulation
makeData <- function(group2loading2 = 0.7, group2intercept1 = 3) {
    pop.model <- paste0(' factor =~ c(0.7, 0.7)*x1 + c(0.7, ',
                        group2loading2, ')*x2 + c(0.8, 0.8)*x3
                        x1 ~ c(3, ', group2intercept1, ')*1
                        x2 ~ c(3, 3)*1
                        x3 ~ c(3, 3)*1 ')
    simulateData(pop.model, standardized = TRUE, sample.nobs = c(100, 100))
}

testMetric <- function(data) {
    mod.config <- ' f =~ x1 + x2 + x3 '
    mod.metric <- ' f =~ c(L1, L1)*x1 + c(L2, L2)*x2 + c(L3, L3)*x3
    f ~~ c(1, NA)*f
    '
    fit.config <- cfa(mod.config, data = data, std.lv = TRUE, group = "group")
    fit.metric <- cfa(mod.metric, data = data, std.lv = TRUE, group = "group")
    LRT <- lavTestLRT(fit.metric, fit.config)
    LRT[2, "Pr(>Chisq)"] < .05 # reject H0 at alpha = .05?
}

testScalar <- function(data) {
    mod.metric <- ' f =~ c(L1, L1)*x1 + c(L2, L2)*x2 + c(L3, L3)*x3
    f ~~ c(1, NA)*f
    '
    mod.scalar <- ' f =~ c(L1, L1)*x1 + c(L2, L2)*x2 + c(L3, L3)*x3
    f ~~ c(1, NA)*f
    x1 ~ c(T1, T1)*1
    x2 ~ c(T2, T2)*1
    x3 ~ c(T3, T3)*1
    f ~ c(0, NA)*1
    '
    fit.metric <- cfa(mod.metric, data = data, std.lv = TRUE, group = "group")
    fit.scalar <- cfa(mod.scalar, data = data, std.lv = TRUE, group = "group")
    LRT <- lavTestLRT(fit.scalar, fit.metric)
    LRT[2, "Pr(>Chisq)"] < .05 # reject H0 at alpha = .05?
}


nReps <- 100

## power analysis for difference in loadings of 0.7 - 0.4 = 0.3
reject.load.4 <- rep(NA, nReps)
set.seed(12345)
for (i in 1:nReps) {
    dat <- makeData(group2loading2 = 0.4)
    reject.load.4[i] <- testMetric(dat)
}
mean(reject.load.4) # empirical power


## power analysis for difference in intercepts of 3 - 2 = 1
reject.int.2 <- rep(NA, nReps)
set.seed(12345)
for (i in 1:nReps) {
    dat <- makeData(group2intercept1 = 2)
    reject.int.2[i] <- testScalar(dat)
}
mean(reject.int.2) # empirical power

[Camille Williams]

Thank you very much