simsem: Problems to conduct simulation study with missing structure and using ML family estimator

[Jaakko Tammilehto]

Hi,

I have recently worked with simsem that is a great R package. I utilized lavaan syntax to form population and analysis models and then conducted MC simulation with 1000 replications using MLR as an estimator (which is my plan also in the actual study). The basic simulations without specifying the missing structure worked just fine.  However, more recently, I have tried to take into account the known missing structure in these simulations. I have used both miss(logit=script, m=0) as well as miss(pmMCAR=0.3, m=0) functions to define the missing structure (see R script below). For some reason, conducting MC simulations is not possible if the ML family estimator is used. The error message is as follow: 

Progress: 1 / 1000 

Error in lav_options_set(opt) : 

 lavaan ERROR: estimator ML for ordered data is not supported yet. Use WLSMV instead.

Using WLSMV as an estimator, I manage to run the MC simulations; however, as I did not specify the thresholds for the indicators (they should be handled as continuous), I got following warning message:

1: In lav_partable_check(lavpartable, categorical = lavoptions$categorical,  : lavaan WARNING: parameter table does not contain thresholds 

So, the questions that I have wondered without solution: 

1. Is it really so that it is not possible to conduct MC simulation with missing structure when using ML family estimator? When I tested this with this example https://github.com/simsem/simsem/blob/master/SupportingDocs/Examples/Version05/ex8/ex8.R I manage the run simulation with ML estimator. 

2. Or is the problem that data is handled as ordered? If so, how can I tell lavaan that all (endogenous) variables are continuous (I assumed that this should be automatic and lavaan only needs to know what variables are categorical (ordered=))? Or is there some other problem in my code?

Thank you for your help!

Best regards,

Jaakko Tammilehto

Here is the code:

#Used packages

install.packages("lavaan")

install.packages("simsem")

install.packages("simstandard")

library("lavaan")

library("simsem")

library("simstandard")

#POPULATION MODEL

AUTmodelPOP.CHILDHOOD<-'

#MOTHER-CHILD RELATIONSHIP AUTONOMY in infancy, middle childhood, and late adolescence

MAUTONOMY1=~0.8*maut1_1+0.7*maut1_2+0.8*maut1_3+0.7*maut1_4 

MAUTONOMY2=~0.8*maut2_1+0.7*maut2_2+0.8*maut2_3+0.7*maut2_4

MAUTONOMY3=~0.8*maut3_1+0.7*maut3_2+0.8*maut3_3+0.7*maut3_4

#FATHER-CHILD RELATIONSHIP AUTONOMY in infancy, middle childhood, and late adolescence

PAUTONOMY1=~0.7*paut1_1+0.8*paut1_2+0.7*paut1_3+0.8*paut1_4

PAUTONOMY2=~0.7*paut2_1+0.8*paut2_2+0.7*paut2_3+0.8*paut2_4

PAUTONOMY3=~0.7*paut3_1+0.8*paut3_2+0.7*paut3_3+0.8*paut3_4

#cross-sectional correlations between MOTHER-CHILD and FATHER-CHILD RELATIONSHIP AUTONOMY

MAUTONOMY1~~0.2*PAUTONOMY1

MAUTONOMY2~~0.15*PAUTONOMY2

MAUTONOMY3~~0.1*PAUTONOMY3

REAPPRAISAL=~1*REAPPRAISAL1

SUPPRESSION=~1*SUPPRESSION1

RUMINATION=~1*RUMINATION1

RUMINATION1~~0.20*RUMINATION1

REAPPRAISAL1~~0.25*REAPPRAISAL1

SUPPRESSION1~~0.25*SUPPRESSION1

#error correlations

maut1_1~~.07*maut2_1

maut1_1~~.03*maut3_1

maut2_1~~.07*maut3_1

maut1_2~~.07*maut2_2

maut1_2~~.03*maut3_2

maut2_2~~.07*maut3_2

maut1_3~~.07*maut2_3

maut1_3~~.03*maut3_3

maut2_3~~.07*maut3_3

maut1_4~~.07*maut2_4

maut1_4~~.03*maut3_4

maut2_4~~.07*maut3_4

paut1_1~~.07*paut2_1

paut1_1~~.03*paut3_1

paut2_1~~.07*paut3_1

paut1_2~~.07*paut2_2

paut1_2~~.03*paut3_2

paut2_2~~.07*paut3_2

paut1_3~~.07*paut2_3

paut1_3~~.03*paut3_3

paut2_3~~.07*paut3_3

paut1_4~~.07*paut2_4

paut1_4~~.03*paut3_4

paut2_4~~.07*paut3_4

#path coefficients

PAUTONOMY2~0.4*PAUTONOMY1+0.15*MAUTONOMY1+0*SEX+0.1*HOITO+0.2*SES

MAUTONOMY2~0.15*PAUTONOMY1+0.4*MAUTONOMY1+0*SEX+0.1*HOITO+0.2*SES

PAUTONOMY3~0.35*PAUTONOMY2+0.15*MAUTONOMY2+0*SEX+0*HOITO+0.1*SES

MAUTONOMY3~0.15*PAUTONOMY2+0.35*MAUTONOMY2+0*SEX+0*HOITO+0.1*SES

REAPPRAISAL~0.13*PAUTONOMY1+0.13*MAUTONOMY1+0.13*PAUTONOMY2+0.13*MAUTONOMY2+0.13*PAUTONOMY3+0.13*MAUTONOMY3+0*SEX+0*HOITO+0.1*SES

SUPPRESSION~(-0.13)*PAUTONOMY1+(-0.13)*MAUTONOMY1+(-0.13)*PAUTONOMY2+(-0.13)*MAUTONOMY2+(-0.13)*PAUTONOMY3+(-0.13)*MAUTONOMY3+(-0.2)*SEX+0*HOITO+(-0.1)*SES

RUMINATION~(-0.13)*PAUTONOMY1+(-0.13)*MAUTONOMY1+(-0.13)*PAUTONOMY2+(-0.13)*MAUTONOMY2+(-0.13)*PAUTONOMY3+(-0.13)*MAUTONOMY3+0.2*SEX+0*HOITO+(-0.1)*SES

#correlations between ER strategies

REAPPRAISAL~~(0.0)*SUPPRESSION

REAPPRAISAL~~(-0.1)*RUMINATION

SUPPRESSION~~0.1*RUMINATION

#correlations between covariates

SES~~0*HOITO

SEX~~0*SES

SEX~~0*HOITO

# Variances

maut1_1 ~~ 0.36 * maut1_1

maut2_1 ~~ 0.36 * maut2_1

maut1_2 ~~ 0.51 * maut1_2

maut2_2 ~~ 0.51 * maut2_2

maut1_3 ~~ 0.36 * maut1_3

maut2_3 ~~ 0.36 * maut2_3

maut1_4 ~~ 0.51 * maut1_4

maut2_4 ~~ 0.51 * maut2_4

paut1_1 ~~ 0.51 * paut1_1

paut2_1 ~~ 0.51 * paut2_1

paut1_2 ~~ 0.36 * paut1_2

paut2_2 ~~ 0.36 * paut2_2

paut1_3 ~~ 0.51 * paut1_3

paut2_3 ~~ 0.51 * paut2_3

paut1_4 ~~ 0.36 * paut1_4

paut2_4 ~~ 0.36 * paut2_4

REAPPRAISAL ~~ 0.7273162365475 * REAPPRAISAL

SUPPRESSION ~~ 0.6873162365475 * SUPPRESSION

RUMINATION ~~ 0.6873162365475 * RUMINATION

HOITO ~~ 1 * HOITO

SEX ~~ 1 * SEX

maut3_1 ~~ 0.36 * maut3_1

maut3_2 ~~ 0.51 * maut3_2

maut3_3 ~~ 0.36 * maut3_3

maut3_4 ~~ 0.51 * maut3_4

paut3_1 ~~ 0.51 * paut3_1

paut3_2 ~~ 0.36 * paut3_2

paut3_3 ~~ 0.51 * paut3_3

paut3_4 ~~ 0.36 * paut3_4

SES ~~ 1 * SES

MAUTONOMY1 ~~ 1 * MAUTONOMY1

MAUTONOMY2 ~~ 0.7435 * MAUTONOMY2

MAUTONOMY3 ~~ 0.791607375 * MAUTONOMY3

PAUTONOMY1 ~~ 1 * PAUTONOMY1

PAUTONOMY2 ~~ 0.7435 * PAUTONOMY2

PAUTONOMY3 ~~ 0.791607375 * PAUTONOMY3

'

#calculating variances for standardized parameters (variances have been added to model syntax)

m_complete <- model_complete(AUTmodelPOP.CHILDHOOD)

cat(m_complete)

#ANALYSIS MODEL

AUTmodel.CHILDHOOD<-'

#MOTHER-CHILD RELATIONSHIP AUTONOMY in infancy, middle childhood, and late adolescence

MAUTONOMY1=~maut1_1+maut1_2+maut1_3+maut1_4

MAUTONOMY2=~maut2_1+maut2_2+maut2_3+maut2_4

MAUTONOMY3=~maut3_1+maut3_2+maut3_3+maut3_4

#FATHER-CHILD RELATIONSHIP AUTONOMY in infancy, middle childhood, and late adolescence

PAUTONOMY1=~paut1_1+paut1_2+paut1_3+paut1_4

PAUTONOMY2=~paut2_1+paut2_2+paut2_3+paut2_4

PAUTONOMY3=~paut3_1+paut3_2+paut3_3+paut3_4

REAPPRAISAL=~1*REAPPRAISAL1

SUPPRESSION=~1*SUPPRESSION1

RUMINATION=~1*RUMINATION1

RUMINATION1~~0.20*RUMINATION1

REAPPRAISAL1~~0.25*REAPPRAISAL1

SUPPRESSION1~~0.25*SUPPRESSION1

#error correlations

maut1_1~~maut2_1

maut1_1~~maut3_1

maut2_1~~maut3_1

maut1_2~~maut2_2

maut1_2~~maut3_2

maut2_2~~maut3_2

maut1_3~~maut2_3

maut1_3~~maut3_3

maut2_3~~maut3_3

maut1_4~~maut2_4

maut1_4~~maut3_4

maut2_4~~maut3_4

paut1_1~~paut2_1

paut1_1~~paut3_1

paut2_1~~paut3_1

paut1_2~~paut2_2

paut1_2~~paut3_2

paut2_2~~paut3_2

paut1_3~~paut2_3

paut1_3~~paut3_3

paut2_3~~paut3_3

paut1_4~~paut2_4

paut1_4~~paut3_4

paut2_4~~paut3_4

#cross-sectional correlations between MOTHER-CHILD AND FATHER-CHILD RELATIONSHIP AUTONOMY

MAUTONOMY1~~PAUTONOMY1

MAUTONOMY2~~PAUTONOMY2

MAUTONOMY3~~PAUTONOMY3

#path coefficients

PAUTONOMY2~PA1*PAUTONOMY1+MA1*MAUTONOMY1+SEX+HOITO+SES

MAUTONOMY2~PA2*PAUTONOMY1+MA2*MAUTONOMY1+SEX+HOITO+SES

PAUTONOMY3~PA3*PAUTONOMY2+MA3*MAUTONOMY2+SEX+HOITO+SES

MAUTONOMY3~PA4*PAUTONOMY2+MA4*MAUTONOMY2+SEX+HOITO+SES

REAPPRAISAL~PARE1*PAUTONOMY1+MARE1*MAUTONOMY1+PARE2*PAUTONOMY2+MARE2*MAUTONOMY2+PARE3*PAUTONOMY3+MARE3*MAUTONOMY3+SEX+HOITO+SES

SUPPRESSION~PASU1*PAUTONOMY1+MASU1*MAUTONOMY1+PASU2*PAUTONOMY2+MASU2*MAUTONOMY2+PASU3*PAUTONOMY3+MASU3*MAUTONOMY3+SEX+HOITO+SES

RUMINATION~ PARU1*PAUTONOMY1+MARU1*MAUTONOMY1+PARU2*PAUTONOMY2+MARU2*MAUTONOMY2+PARU3*PAUTONOMY3+MARU3*MAUTONOMY3+SEX+HOITO+SES

#correlation between ER strategies

REAPPRAISAL ~~ SUPPRESSION

REAPPRAISAL ~~ RUMINATION

SUPPRESSION ~~ RUMINATION

'

#missing data structure

script<-'maut1_1~p(0)

maut1_2~p(0.0)

maut1_3~p(0.0)

maut1_4~p(0.0)

maut2_1~p(0.112)

maut2_2~p(0.112)

maut2_3~p(0.112)

maut2_4~p(0.112)

maut3_1~p(0.174)

maut3_2~p(0.174)

maut3_3~p(0.174)

maut3_4~p(0.174)

paut1_1~p(0.073)

paut1_2~p(0.073)

paut1_3~p(0.073)

paut1_4~p(0.073)

paut2_1~p(0.452)

paut2_2~p(0.452)

paut2_3~p(0.452)

paut2_4~p(0.452)

paut3_1~p(0.346)

paut3_2~p(0.346)

paut3_3~p(0.346)

paut3_4~p(0.346)

REAPPRAISAL1~p(0.196)

SUPPRESSION1~p(0.196)

RUMINATION1~p(0.196)

SEX~p(0.0)

HOITO~p(0.0)

SES~p(0.0)'

missingness<-miss(logit=script, m=0)

AUTanalysis.POPCHILDHOOD.CHILDHOOD <-sim(1000, n =546, model=AUTmodel.CHILDHOOD, seed=565, generate=AUTmodelPOP.CHILDHOOD, lavaanfun = "sem", missingness, std.lv=T, cilevel=0.90, estimator="MLR")

Progress: 1 / 1000 

Error in lav_options_set(opt) : 

  lavaan ERROR: estimator ML for ordered data is not supported yet. Use WLSMV instead.

Progress: 2 / 1000 

Error in lav_options_set(opt) : 

  lavaan ERROR: estimator ML for ordered data is not supported yet. Use WLSMV instead.

#the same thing with using MCAR

miss.model <- miss(pmMCAR=0.3, m=0)

AUTanalysis.POPCHILDHOOD.CHILDHOOD <-sim(1000, n =546, model=AUTmodel.CHILDHOOD, seed=565, generate=AUTmodelPOP.CHILDHOOD, lavaanfun = "sem", miss.model, std.lv=T, cilevel=0.90, estimator="MLR")

Progress: 1 / 1000 

Error in lav_options_set(opt) : 

  lavaan ERROR: estimator ML for ordered data is not supported yet. Use WLSMV instead.

Progress: 2 / 1000 

Error in lav_options_set(opt) : 

  lavaan ERROR: estimator ML for ordered data is not supported yet. Use WLSMV instead.

[Alex Schoemann]

Hi Jaakko,

This is a bit odd. Looking over your code, the first thing I'd check is how you're calling the missing data object in the sim function. You don't include the name of the argument for either of the missing data templates (missingness or miss.model), thus I'm not sure exactly which argument the missing data information is being passed to.  I'd add miss = missingness (or miss.model) to your sim command and see if that fixes things.

Alex

[Jaakko Tammilehto]

Thank you for your help Alex,

It seems to help! However, at the same time, the simulation is really slow (i.e., 4 replication in 5 minutes). But perhaps this is because of the model complexity, and there is nothing to do about it.

Best,
Jaakko

--
You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.
To post to this group, send email to lav...@googlegroups.com.
Visit this group at https://groups.google.com/group/lavaan.
To view this discussion on the web visit https://groups.google.com/d/msgid/lavaan/ce9564b9-3f54-486b-9b92-7b068b201d42%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Alex Schoemann]

Hmmm, that does sound slow. It's possible that everything is working fine and that's what's happening, but usually if things are going that slowly it implies that there are convergence difficulties with each model. I'd recommend running a small number of replications (e.g 5-10) and looking at the results before diving into a full 1000. That way if there's something wonky with the models that is causing convergence issues you can diagnose it more quickly.

Alex

To unsubscribe from this group and stop receiving emails from it, send an email to lav...@googlegroups.com.