sim() function does not work (anymore) with my lavaan models

[doro]

Dear all,

I'm trying to run several scripts for power analysis which were working perfectly last year (June) - but seem completely broken now. I have the newest versions of R, lavaan, and simsem installed, but I didn't change anything in the syntax. I've tried to find any hints in the version history but with no success. For example, I've tried to run

act.power <- sim(nRep=10000, generate=mod1.pop2, model=mod1, n =104, lavaanfun = "cfa", multicore=TRUE, seed=565)

where mod.pop2 is a CFA with values from a pilot study. The very same syntax worked fine in June. Now I get the error:

Error in x$timing : $ operator is invalid for atomic vectors

For another syntax (power analysis for a latent change model), I get the error:

Error in names(coef) <- lab : 

  'names' attribute [19] must be the same length as the vector [15]

 

Help and any ideas would be very much appreciated!

Best,

Dorota

[Terrence Jorgensen]

I have the newest versions of R, lavaan, and simsem installed, but I didn't change anything in the syntax. 

I'm not sure if you mean the latest versions available on CRAN or the development versions.  Could you please provide your sessionInfo() output, and the full R script so I can reproduce your error and see your version numbers?

Terrence D. Jorgensen

Postdoctoral Researcher, Methods and Statistics

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

UvA web page: http://www.uva.nl/profile/t.d.jorgensen

[doro]

Hi Terrence,

here you are:

# Data generation model

mod1.pop2 <- ' bem =~ 0.5*m1 + 0.2*m2+ 0.6*m3+ 0.7*m4

nabl =~ 0.8*m5 + 0.1*m6 + 0.9*m7

sksm =~ 0.6*m8 + 0.7*m9 + 0.4*m10

ar =~ 0.8*m11 + 0.6*m12 + 0.7*m13 + 0.8*m14 + 0.8*m15 + 0.7*m16 + 1.0*m17

eg =~ 0.5*m18 + 0.4*m19 + 0.9*m20 + 0.9*m21 + 0.9*m22

sr =~ 0.8*m23 + 0.9*m24 + 0.8*m25 + 0.7*m26

adlb =~ 0.9*m27 + 0.8*m28 + 0.8*m29

vert =~ 1.2*m30 + 1.1*m31 + 0.7*m32

bem ~~ 1*bem

nabl ~~ 1*nabl

sksm ~~ 1*sksm

ar ~~ 1*ar

eg ~~ 1*eg

sr ~~ 1*sr

adlb ~~ 1*adlb

vert ~~ 1*vert

m1 ~~ 1.4*m1

m2 ~~ 1.0*m2

m3 ~~ 1.0*m3

m4 ~~ 1.1*m4

m5 ~~ 1.0*m5

m6 ~~ 0.9*m6

m7 ~~ 0.5*m7

m8 ~~ 1.0*m8

m9 ~~ 0.8*m9

m10 ~~ 0.8*m10

m11 ~~ 0.9*m11

m12 ~~ 0.9*m12

m13 ~~ 1.0*m13

m14 ~~ 0.6*m14

m15 ~~ 0.7*m15

m16 ~~ 0.9*m16

m17 ~~ 0.5*m17

m18 ~~ 0.9*m18

m19 ~~ 0.8*m19

m20 ~~ 0.4*m20

m21 ~~ 1.0*m21

m22 ~~ 0.4*m22

m23 ~~ 0.4*m23

m24 ~~ 0.6*m24

m25 ~~ 1.1*m25

m26 ~~ 0.8*m26

m27 ~~ 0.4*m27

m28 ~~ 0.8*m28

m29 ~~ 0.7*m29

m30 ~~ 0.5*m30

m31 ~~ 0.5*m31

m32 ~~ 0.6*m32

bem ~~ -0.3*nabl

bem ~~ 0.05*sksm

bem ~~ 0.9*ar

bem ~~ 1.3*eg

bem ~~ 1.1*sr

bem ~~ 0.9*adlb

bem ~~ 0.6*vert

nabl ~~ 0.3*sksm

nabl ~~ -0.1*ar

nabl ~~ -0.4*eg

nabl ~~ 0.01*sr

nabl ~~ 0.2*adlb

nabl ~~ 0.3*vert

sksm ~~ 0.4*ar

sksm ~~ 0.03*eg

sksm ~~ 0.1*sr

sksm ~~ -0.2*adlb

sksm ~~ 0.1*vert

ar ~~ 0.8*eg

ar ~~ 0.8*sr

ar ~~ 0.7*adlb

ar ~~ 0.7*vert

eg ~~ 0.8*sr

eg ~~ 0.7*adlb

eg ~~ 0.4*vert

sr ~~ 0.9*adlb

sr ~~ 0.7*vert

adlb ~~ 0.7*vert'

# analysis model

mod1 <- 'bem =~ m1 + m2+ m3+ m4

        nabl =~ m5 + m6 + m7

sksm =~ m8 + m9 + m10

ar =~ m11 + m12 + m13 + m14 + m15 + m16 + m17

eg =~ m18 + m19 + m20 + m21 + m22

sr =~ m23 + m24 + m25 + m26

adlb =~ m27 + m28 + m29

vert =~ m30 + m31 + m32

'   

detach("package:psych", unload=TRUE)

require(simsem)

act.power <- sim(nRep=10000, generate=mod1.pop2, model=mod1, n =104, lavaanfun = "cfa", multicore=TRUE, seed=565)

And my session info:

R version 3.3.3 (2017-03-06)

Platform: x86_64-apple-darwin13.4.0 (64-bit)

Running under: macOS Sierra 10.12.3

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] semTools_0.4-14    Amelia_1.7.4       Rcpp_0.12.9        quantreg_5.29      simsem_0.5-13      lavaan_0.5-23.1097

loaded via a namespace (and not attached):

 [1] quadprog_1.5-5     lattice_0.20-34    psych_1.6.12       MASS_7.3-45        grid_3.3.3         MatrixModels_0.4-1

 [7] stats4_3.3.3       SparseM_1.74       Matrix_1.2-8       pbivnorm_0.6.0     tools_3.3.3        foreign_0.8-67    

[13] parallel_3.3.3     mnormt_1.5-5       MBESS_4.2.0 

Thanks for checking!!

Best,

Dorota

[Terrence Jorgensen]

OK, this script "works" for me, in the sense that I don't get the errors you saw when I run it using the latest versions available.  In a fresh R session, install all 3 packages:

install.packages("lavaan", repos = "http://www.da.ugent.be", type = "source")
devtools::install_github("simsem/simsem/simsem")
devtools::install_github("simsem/semTools/semTools")

Now, I still got an error, but not from simsem.  Your population model is not possible, in the sense that the model-implied covariance matrix is not positive definite.  A good trick to check whether you have specified a plausible set of parameters is to check the model-implied population covariance matrix:

(impliedSigma <- fitted(lavaan(mod1.pop2)))

Using your script above, you will notice this error:

Error in lav_start_check_cov(lavpartable = lavpartable, start = START) : 

  lavaan ERROR: please provide better fixed values for (co)variances;

                variables involved are: bem eg

In addition: Warning message:

In lav_start_check_cov(lavpartable = lavpartable, start = START) :

  lavaan WARNING: starting values imply a correlation larger than 1;

                  variables involved are: bem eg

The first problem it notices are with bem and eg, but I don't think those are the only problematic parameters.  After changing these 2 parameters:

bem ~~ .3*eg

bem ~~ .1*sr

I still got a warning that Sigma is not positive definite, so it is still not possible to draw samples from that population. You can check yourself that the determinant is negative:

det(fitted(lavaan(mod1.pop2))$cov)

You can keep adjusting the population values until you get a positive determinant.  Since you are specifying the population using lavaan syntax, you can also check whether the population is acceptable by using lavaan's simulateData() function to see whether you can draw one sample from it. 

simulateData(mod1.pop2, sample.nobs = 10)

If you can, then you should be able to provide that same script to simsem.

Good luck,

[doro]

Hi Terrence,

thanks so much for all the help!

Yes, indeed, I realize there are issues with the CFA. But the other script (power for the latent change model) was even more important at the moment. Now it works like a charm with the devtools versions! Next time I will install those first - before trying for hours to figure out why something doesn't work anymore...

Thanks again!

Best, Dorota

[Elisabeth Zureck]

Dear Terrence,
We are trying to do a power analysis with the simsem package. As you suggested we first installed the three packages. As I got an error message with my complete model I tried the following reduced model:

For example we had the following analysis model:

The analysis model was:
analysis-model <- 'CM1 =~ CM1.Ind1 + a*CM1.Ind2
            CM2 =~ CM2.Ind1 + a*CM2.Ind2
            CM1.Ind1 + CM2.Ind1 ~ 0*1
            CM1.Ind2 + CM2.Ind2 ~ c*1
            CM1 + CM2 ~ 1'

For the population model we inserted the parameter estimates resulting from the analysis model:
Pop.Mod <- 'CM1 =~0.95*CM1.Ind1 + 1.06*CM1.Ind2
            CM2 =~ 0.95*CM2.Ind1 + 1.06*CM2.Ind2

            # Intercepts
            CM1.Ind1 ~ 0*1
            CM2.Ind1 ~ 0*1
            CM1.Ind2 ~ -0.31*1
            CM2.Ind2 ~ -0.31*1
 
            CM1 ~ 2.93*1
            CM2 ~ 2.94*1

            CM1 ~~ 1*CM1
            CM2 ~~ 1*CM2

            CM1.Ind1 ~~ 0.37*CM1.Ind1 
            CM1.Ind2 ~~ 0.41*CM1.Ind2 
            CM2.Ind1 ~~ 0.45*CM2.Ind1 
            CM2.Ind2 ~~ 0.31*CM2.Ind2 

            CM1 ~~ 0.53*CM2'
and afterwards applied the sim() function on these models
lc.n <- sim(model=LC.CmGs, n = rep(seq(50,250,50),5), generate = Pop.Mod, lavaanfun = "cfa", std.lv = TRUE)

At the end of the summary-output we got the following 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.

I just can't see were the difference in between the two models might be. The summaryPopulation() Funktion also only lists the parameters I specified in my analysis- and my population model.

Help and any ideas would be very much appreciated!

Best,

Elisabeth

[Alex Schoemann]

Hi Elisabeth,

It looks like the sim command isn't using the analysis model you posted. In the sim command you're using the model LC.CmGs as the analysis model, not the model you posted in the syntax. Changing the code to use the model you specified above worked for me without that error.

Hope this helps,

Alex

[doro]

Hi Alex,

Thanks for the prompt answer! The different analysis model is a typo (we have estimated so many models...). I've tried just now with the correct analysis model and developer versions of all packages installed (as suggested by Terrence) but it still says: 

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.

But the summaryPopulation function gives correct parameters.

Any ideas? We are quite desperate :(

Best,

Dorota

[Alex Schoemann]

I think this might be because in your  analysis model has equality constraints in it and your generating model does not have the constraints in the syntax. If you add the constraints to the syntax for the generating model, I think, the warning should go away. Try the generating model below (which has the equality constraints in place)

Pop.Mod <- 'CM1 =~0.95*CM1.Ind1 + 1.06*CM1.Ind2 + a*CM1.Ind2
            CM2 =~ 0.95*CM2.Ind1 + 1.06*CM2.Ind2 + a*CM2.Ind2

            # Intercepts
            CM1.Ind1 ~ 0*1
            CM2.Ind1 ~ 0*1
            CM1.Ind2 ~ -0.31*1
            CM2.Ind2 ~ -0.31*1
 
            CM1 ~ 2.93*1
            CM2 ~ 2.94*1

            CM1 ~~ 1*CM1
            CM2 ~~ 1*CM2

            CM1.Ind1 ~~ 0.37*CM1.Ind1  
            CM1.Ind2 ~~ 0.41*CM1.Ind2  
            CM2.Ind1 ~~ 0.45*CM2.Ind1  
            CM2.Ind2 ~~ 0.31*CM2.Ind2  

            CM1 ~~ 0.53*CM2'

[Elisabeth Zureck]

Hi Alex,
Thank you for your suggestion! Indeed, the critical issue was the constraints we set by using letters in the analysis model that were not included in the population model.
Including the letters in the population model just like we did it in the analysis model and defining them in an extra line worked for the "small" and the more complex model.
e.g. CM1 =~ 0.95*CM1.Ind1 + b*CM1.Ind2
CM2 =~ 0.95*CM2.Ind1 + b*CM2.INd2
b == 1.06

Thank you very much for your help!

Best,
Elisabeth

Am Freitag, 10. März 2017 16:52:19 UTC+1 schrieb doro: