invoking lavaan() via lavaan.mi(): Error ... subscript out of bounds

[Christian Seubert]

Hi all,

this is my first post and I hope this is the right place for my problem. I want to estimate a structural equation model with multiply imputed data using lavaan.mi() from the semTools package. The data were previously imputed with the mice package, so I extracted the five imputed data sets with:

imputed<-list()

for(i in 1:5) {

  imputed[[i]] <- complete(imp, i)

}

Then I tried to fit the model:

sem.m3.fit<-lavaan.mi(sem.m3, data=imputed, chi="all", likelihood="wishart", int.ov.free=T,

int.lv.free=F, auto.fix.first=T, auto.fix.single=T, auto.var=T, auto.cov.lv.x=F, auto.cov.y=F, auto.th=T, auto.delta=T)

This leaves me with an error message:

Error in `[<-`(`*tmp*`, rhs, lhs, value = 1) : subscript out of bounds

However, when I fit the model with just one of the imputed data sets,

sem.m3.fit<-lavaan(sem.m3, data=imputed[[1]], likelihood="wishart", int.ov.free=T,

int.lv.free=F, auto.fix.first=T, auto.fix.single=T, auto.var=T, auto.cov.lv.x=F, auto.cov.y=F, auto.th=T, auto.delta=T)

the model is fitted without any error. I varied the arguments in lavaan.mi() that are passed to lavaan() and tried sem.mi() as an alternative - to no avail. Therefore, if I didn't make a mistake, I suspect the problem might be in lavaan.mi() [or the function it wraps, runMI()] rather than lavaan(). Still, as Yves Rosseel also contributes to the semTools package, I suspected this might be the right place to ask for help. So: Can anybody help me out? Thanks in advance.

Christian

[Alex Schoemann]

Hmmm, this is a bit hard to diagnose without a but more information. There are a couple of things you might try:

1. Does running lavaan work on each imputed data set (you tried with the first but the problem might crop up in another one).

2. When you run lavaan.mi, try setting chi="none" If this works, it will help us pinpoint the issue.

3. If you're willing, posting the code you're using to define the model is always helpful.

Alex

[Christian Seubert]

Thanks for the quick response.

1. lavaan does indeed work on each imputed data set, it also works on the non-imputed dataset (complete cases only).

2. I tried all variants of chi="...", all with the same result.

3. Here's my model:

sem.m3<-'

outcome1 ~ 2ndOrderF1 + 2ndOrderF2 + mediator1

outcome2 ~ 2ndOrderF3 + 2ndOrderF2 + mediator2

mediator1 ~ 2ndOrderF1 + 2ndOrderF2

mediator2 ~ 2ndOrderF3 + 2ndOrderF2

2ndOrderF1 =~ 1stOrderF1.a + 1stOrderF1.b

2ndOrderF2 =~ 1stOrderF2.a + 1stOrderF2.b

2ndOrderF3 =~ 1stOrderF3.a + 1stOrderF3.b

1stOrderF1.a =~ ind1+ind2+ind3

1stOrderF2.a =~ ind4+ind5+ind6

1stOrderF3.a =~ ind7+ind8+ind9

1stOrderF1.b =~ ind10+ind11+ind12

1stOrderF2.b =~ ind13+ind14+ind15

1stOrderF3.b =~ ind16+ind17+ind18

mediator1 =~ ind19+ind20+ind21

mediator2 =~ ind22+ind23+ind24

outcome1 =~ ind25+ind26+ind27

outcome2 =~ ind28+ind29+ind30

2ndOrderF1 ~~ 2ndOrderF2

2ndOrderF1 ~~ 2ndOrderF3

2ndOrderF3 ~ 2ndOrderF2

mediator1 ~~ mediator2

'

[Christian Seubert]

Hi Alex,

regarding your remarks in the post that has now been deleted:

When I try to run the code you suggest, I get this error

test<-semTools:::runlavaanMI( data=imputed, syntax=sem.m3, fun = "lavaan", likelihood="wishart", int.ov.free=T, int.lv.free=F, auto.fix.first=T, auto.fix.single=T, auto.var=T, auto.cov.lv.x=F, auto.cov.y=F, auto.th=T, auto.delta=T)
> Error in list(model = syntax, data = MIdata) : 'MIdata' is missing

Does the list containing the imputed data sets have to have a specific structure?

And, yes, thank you for the remark, the second to last line should read ~~ instead of ~. It happened when I renamed the variables to copy them here, in R, it is specified correctly.

[yrosseel]

On 10/18/2013 08:46 AM, Christian Seubert wrote:
> Hi Alex,
>
> regarding your remarks in the post that has now been deleted:

Administrative note: I never delete a post, but sometimes, the google
groups system 'marks' a post as 'spam' (for reasons that go beyond me).
I get a notice, and I have to 'unmark' that post, after which it gets
posted anyway (but with some delay).

Yves.

[Christian Seubert]

Thanks, Yves, for the clarification. Even though censorship was not on my mind when I wrote that, it's informative to be aware of the behavior of google groups.

[Alex Schoemann]

Ah, there was an error in the code I sent to you (maybe that's why Google thought it was spam:)). The data command below should read data=imputed[[1]] (runlavaanMI only works on a single dataset). 

I believe that the list of imputed data sets should simply be a list of data frames. In fact we use the same code you did to pull data frames out of mice. So I don't think that's an issue. 

[Christian Seubert]

Ok, I had already tried that and verified it again right now - using a single data set - whether one of the five imputed or a non-imputed, complete-cases-only data set - as input still gives me this error message: "Error in list(model = syntax, data = MIdata) : 'MIdata' is missing." I use exactly the syntax you gave me. The data entered is of class "data.frame"..

[Alex Schoemann]

Oh man, instead of  data=imputed[[1]], it should read MIdata=imputted[[1]] 

I need to drink more coffee right now. I am not great in the morning!

[Christian Seubert]

Darn! I could have figured that out, too :-) It works now, with any of the imputed data sets. What's next?

[Alex Schoemann]

Okay, I've pinpointed the problem playing with some of my own data. The issue is in the specification of the higher order factors. Currently, when higher order factors are specified using the =~ operator to define the higher order factor kicks back an error (this stems from the specification of the null and saturated models used when computing model fit indices). I'll put this on the list of issues to fix. 

Fortunately there's an easy fix. Just switch the =~ operator to ~ for the second order factors and (with a few other adjustments) you should be all set.

So this:

2ndOrderF1 =~ 1stOrderF1.a + 1stOrderF1.b

2ndOrderF2 =~ 1stOrderF2.a + 1stOrderF2.b

2ndOrderF3 =~ 1stOrderF3.a + 1stOrderF3.b

Becomes:

#Specify 2nd order latents with a loading of 0

#Trick to estimate the model

2ndOrderF1 =~ 0*ind10

2ndOrderF2 =~ 0*ind10

2ndOrderF3 =~ 0*ind10

#Use regression instead of loading operator

1stOrderF1.a ~  1*2ndOrderF1 

1stOrderF1.b ~  2ndOrderF1 

1stOrderF2.a ~ 1*2ndOrderF2 

1stOrderF2.b ~ 2ndOrderF2

1stOrderF3.a ~ 1*2ndOrderF3

1stOrderF3.b ~ 2ndOrderF3

#don't allow 1st order factors to covary (as is default)

1stOrderF1.a ~~ 0*1stOrderF1.b + 0*1stOrderF2.a + 0*1stOrderF2.b + 0*1stOrderF3.a + 0*1stOrderF3.b

1stOrderF1.b ~~ 0*1stOrderF2.a + 0*1stOrderF2.b + 0*1stOrderF3.a + 0*1stOrderF3.b

1stOrderF2.a ~~ 0*1stOrderF2.b + 0*1stOrderF3.a + 0*1stOrderF3.b

1stOrderF2.b ~~ 0*1stOrderF3.a + 0*1stOrderF3.b

1stOrderF3.a ~~ 0*1stOrderF3.b

[Christian Seubert]

It seems like by now, you've had enough coffee ;-) lavaan.mi() fits the model successfully and I've verified that the models are equivalent using a single data set.

Awesome! Thanks, Alex, for the great support. It is very much appreciated.

Best regards,
Christian