Search code examples
rpurrrr-lavaan

R - use simsem within map

I have a list of more than 100 SEM models computed in lavaan. For this example, I'm going to use just two models:

fit.model1.cfa <- '
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    dem65 =~ y5 + y6 + y7 + y8
    ind60 ~~ ind60
    dem60 ~~ dem60
    dem65 ~~ dem65' %>%
       sem(PoliticalDemocracy)

fit.model2.cfa <- '
    ind60 =~ x1 + x2
    dem60 =~ y1 + y2 + y3
    dem65 =~ y5 + y6 + y7
    ind60 ~~ ind60
    dem60 ~~ dem60
    dem65 ~~ dem65' %>%
       sem(PoliticalDemocracy)

It's not that much important what the models actually are. Now we will join them in a list:

models <- list(fit.model1.cfa, fit.model2.cfa)

I want to do post-hoc power analysis on my models and for that purpose I am using the following for loop:

pwr <- list()

for(i in 1:2){
   sim(nRep = 1000, model = models[[i]], n = models[[i]]@SampleStats@ntotal, generate = models[[i]]) %>%
       getPower() -> pwr[[i]]
}

However, I would like to avoid the for loop and instead use the map function. I was trying the following:

models %>%
   map(., sim, nRep = 1000, model = ., n = pluck(pluck(., "SampleStats"), "ntotal"), generate = .) %>%
   map(getPower)

However, that is not working. So what am I doing wrong?

Solution

  • This is the first time I heard about lavaan, simsem and its models, but I'm trying to help you with the map. I made a simple solution with nRep = 10 just to test and you can replicate it with nRep = 1000. Also, lots of warnings and progress were returned, but I deleted them in the example below for the sake of cleanliness.

    library(tidyverse, verbose = F)
library(lavaan)
library(simsem)

fit.model1.cfa <- '
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    dem65 =~ y5 + y6 + y7 + y8
    ind60 ~~ ind60
    dem60 ~~ dem60
    dem65 ~~ dem65' %>%
  sem(PoliticalDemocracy)

fit.model2.cfa <- '
    ind60 =~ x1 + x2
    dem60 =~ y1 + y2 + y3
    dem65 =~ y5 + y6 + y7
    ind60 ~~ ind60
    dem60 ~~ dem60
    dem65 ~~ dem65' %>%
  sem(PoliticalDemocracy)

models <- list(model1 = fit.model1.cfa, model2 = fit.model2.cfa)


power <- map(models, function(x){

  n <- x@SampleStats@ntotal

  z <- sim(model = x, nRep = 10, n = n, generate = x)

  getPower(z)

})

power

#> $model1
#>    ind60=~x2    ind60=~x3    dem60=~y2    dem60=~y3    dem60=~y4    dem65=~y6 
#>          1.0          1.0          1.0          1.0          1.0          1.0 
#>    dem65=~y7    dem65=~y8 ind60~~ind60 dem60~~dem60 dem65~~dem65       x1~~x1 
#>          1.0          1.0          1.0          1.0          1.0          1.0 
#>       x2~~x2       x3~~x3       y1~~y1       y2~~y2       y3~~y3       y4~~y4 
#>          0.5          1.0          1.0          1.0          1.0          1.0 
#>       y5~~y5       y6~~y6       y7~~y7       y8~~y8 ind60~~dem60 ind60~~dem65 
#>          1.0          1.0          1.0          1.0          1.0          1.0 
#> dem60~~dem65 
#>          1.0 
#> 
#> $model2
#>    ind60=~x2    dem60=~y2    dem60=~y3    dem65=~y6    dem65=~y7 ind60~~ind60 
#>          1.0          1.0          1.0          1.0          1.0          1.0 
#> dem60~~dem60 dem65~~dem65       x1~~x1       x2~~x2       y1~~y1       y2~~y2 
#>          1.0          1.0          0.4          0.4          0.8          1.0 
#>       y3~~y3       y5~~y5       y6~~y6       y7~~y7 ind60~~dem60 ind60~~dem65 
#>          1.0          1.0          1.0          1.0          0.9          1.0 
#> dem60~~dem65 
#>          1.0