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

71 views
Skip to first unread message

Jaakko Tammilehto

unread,
Jul 14, 2019, 1:21:09 PM7/14/19
to lavaan
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

unread,
Jul 14, 2019, 8:50:33 PM7/14/19
to lavaan
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

unread,
Jul 15, 2019, 1:07:35 AM7/15/19
to lav...@googlegroups.com
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

unread,
Jul 16, 2019, 9:08:07 AM7/16/19
to lavaan
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.
Reply all
Reply to author
Forward
0 new messages