Error message for mxCompare() with boot = T

Posted on
No user picture. Veronica_echo Joined: 02/23/2018
Hi,

I want to compare a 2-class growth mixture model to a latent growth curve model in individually-varying time points framework (i.e., definition variables), and here is the code and error message and OpenMx version:

> mxCompare(Weight2, Weight1)
base comparison ep
1 Growth Mixture Model, fixed knots 31
2 Growth Mixture Model, fixed knots Estimate a fixed knot 15
minus2LL df AIC diffLL diffdf p
1 4775.720 1464 1847.720 NA NA NA
2 4859.904 1480 1899.904 84.18401 16 2.906173e-11
> mxCompare(Weight2, Weight1, boot = T)

> mxCompare(Weight2, Weight1, boot = T, replications = 10)
Error in collectStatistics1(otherStats, ref, other, bootPair) :
Fewer than 3 replications are available.

> mxVersion()
OpenMx version: 2.9.9.248 [GIT v2.9.9-248-gd602399]
R version: R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0
MacOS: 10.13.6
Default optimiser: CSOLNP
NPSOL-enabled?: Yes
OpenMP-enabled?: Yes

May I know how to address this issue? Thank you in advance!

Replied on Fri, 07/13/2018 - 08:33
Picture of user. jpritikin Joined: 05/23/2012

You run mxCompare twice:


mxCompare(Weight2, Weight1, boot = T)
mxCompare(Weight2, Weight1, boot = T, replications = 10)

Why? And why only 10 replications the second time?

Replied on Mon, 07/16/2018 - 11:59
No user picture. Veronica_echo Joined: 02/23/2018

In reply to by jpritikin

Thanks for your reply. The
Thanks for your reply. The default setting of replication is 400 and I let it run overnight but I got nothing. So for the second number, I just tried a smaller replication number to see if I can get something.
Replied on Mon, 07/16/2018 - 12:24
Picture of user. jpritikin Joined: 05/23/2012

In reply to by Veronica_echo

What do you mean that you "got nothing?" By using the previousRun argument, you can run replication incrementally. You can run 10 then run another 15 without losing the first 10.
Replied on Mon, 06/21/2021 - 05:06
No user picture. cjvanlissa Joined: 04/10/2019

In reply to by jpritikin

I'm curious if this was ever resolved. I also "get nothing" except an error when running the following, where m1 and m2 are two mixture models:


tmp <- mxCompare(m1, m2, boot=TRUE)
Error in collectStatistics1(otherStats, ref, other, bootPair) :
Fewer than 3 replications are available.

Replied on Mon, 06/21/2021 - 16:48
Picture of user. jpritikin Joined: 05/23/2012

Looks like your models are not converging. Take a look at `attr(tmp@results, 'bootData')`. There will be a list of data.frames, one for each model. Take a look at the `statusCode` column.
Replied on Tue, 06/22/2021 - 10:32
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by cjvanlissa

I suspect you'll need to change mxOption 'Status OK' from its on-load default. This is what that option looks like on load:

> options()$mxOptions$`Status OK`
[1] OK OK/green
attr(,"mxFactor")
[1] TRUE
10 Levels: OK < OK/green < infeasible linear constraint < infeasible non-linear constraint < ... < infeasible start
> levels(options()$mxOptions$`Status OK`)
[1] "OK" "OK/green" "infeasible linear constraint"
[4] "infeasible non-linear constraint" "iteration limit/blue" "not convex/red"
[7] "nonzero gradient/red" "bad deriv" "internal error"
[10] "infeasible start"

For example, if your models involve thresholds, you might want to set the option as follows:

mxOption(NULL,"Status OK",c("OK","OK/green","not convex/red","nonzero gradient/red"))

Note that this is only a fix if your models are failing to converge in most of the bootstrap replications (which is possible, even if they do converge when fitted to your real dataset).
Replied on Tue, 06/22/2021 - 11:23
Picture of user. jpritikin Joined: 05/23/2012

> Since the call to mxCompare() stops with an error, the object tmp is never created, so I cannot examine its attributes.

Gah, that's super annoying. I removed that code. Can you install v2.19.5-26-g8a688cd29 ?

Replied on Wed, 06/23/2021 - 03:14
No user picture. cjvanlissa Joined: 04/10/2019

In reply to by jpritikin

Thank you sincerely. I am experiencing serious difficulties installing the package from GitHub; will come back to this problem when I have a bit more time to set everything up correctly and debug.
Replied on Wed, 06/23/2021 - 11:07
No user picture. cjvanlissa Joined: 04/10/2019

I think I know why this fails. mxAutoStart requires mxData() inside each component of a mixture model. mxCompare(boot = TRUE), by contrast, fails when there are mxData() inside of each component of the mixture model. Thus, one can currently have either autoStarts, or mxCompare, but not both. Please see the reproducible example in the attached file. Thank you for your support.
File attachments
Replied on Wed, 06/23/2021 - 12:36
Picture of user. jpritikin Joined: 05/23/2012

In reply to by cjvanlissa

I bet the problem is here,
```
classes = lapply(classes, mxAutoStart, type = "ULS")
coef(classes[[1]]) - coef(classes[[2]])
```
When you start both mixtures at the same point in the parameter space, the optimizer has trouble getting anywhere because it isn't clear that there is more than 1 mixture component.
Replied on Wed, 06/23/2021 - 13:41
No user picture. cjvanlissa Joined: 04/10/2019

In reply to by jpritikin

Dear Joshua, I respectfully believe that this interpretation can be excluded, for multiple reasons:
1) The model does arrive at the correct solution for the 2-class mixture, despite having the same starting values for both classes. (Correct solution as validated by Mplus and Mclust).
2) When I manually provide the same starting values for both mixtures, but also put a mxData() object inside both classes, mxCompare(boot = TRUE) also fails. Moreover, when I reverse the order of the models, mxCompare succeeds, but I notice that the number of degrees of freedom is incorrect. Cases from all the datasets are added up.
3) When I manually provide the same starting values for both mixtures, and do not put a mxData() object inside both classes, mxCompare(boot = TRUE) works. Thus proving that putting mxData() in the submodels breaks mxCompare(boot = TRUE)

See attached script for a demonstration.

File attachments
Replied on Thu, 06/24/2021 - 11:06
Picture of user. jpritikin Joined: 05/23/2012

In reply to by cjvanlissa

Of course you are correct. If there is more than one mxData in the model then each mxData is bootstrap resampled independently. So instead of a mixture model, you are fitting a nonsensical multiple group model.

I don't think there is a bug here. That's how it works by design.

Replied on Thu, 06/24/2021 - 12:26
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by jpritikin

But Joshua, isn't `mxCompare()` doing a bootstrap likelihood-ratio test? The bootstrap LRT doesn't re-sample the actual data. It generates data parametrically under the null model, and fits both the null and alternative models to each generated dataset. The _p_-value is then the proportion of replications in which the LRT statistic exceeds the real-data LRT statistic.

Remember, mixture models were the whole motivation behind implementing the bootstrap LRT in OpenMx

Replied on Thu, 06/24/2021 - 13:30
Picture of user. jpritikin Joined: 05/23/2012

In reply to by AdminRobK

> It generates data parametrically under the null model

Yeah, but that's not an important distinction. I agree that my diagnosis was not quite correct.

If you run `mxGenerateData` you will see that it only finds at the top-level mxData and ignores the mxData in the mixture components. If the mixture components contain an mxData then they will use it instead of traversing up the model tree to the top-level mxData. So the result is similar to what I suggested before, the mixture components never see the bootstrapped datasets and just fit against the same datasets over and over again.

Replied on Fri, 06/25/2021 - 03:12
No user picture. cjvanlissa Joined: 04/10/2019

In reply to by jpritikin

Thank you for the detailed explanation; it makes perfect sense that the individual classes should not contain mxData(). I guess it does not make that much sense that this should break mxCompare(boot = TRUE) though, as the presence of mxData() in classes should be irrelevant to the BLRT.

Secondly, I wonder how to set starting values for a mixture model using mxAutoStart() without adding mxData() to the individual classes. Do you think mxAutoStart() should be able to handle such cases where mxData() is only present in the top-level mxModel()?

Replied on Thu, 06/24/2021 - 13:46
Picture of user. AdminRobK Joined: 01/24/2014

In other words, `mxCompare()` with `boot=TRUE` doesn't actually use the MxData objects. It only needs the real-data -2logL, and an expectation object.

Edit: To be clear, both the null and alternative models need to have valid expectation objects.

Replied on Thu, 06/24/2021 - 15:20
Picture of user. jpritikin Joined: 05/23/2012

In reply to by AdminRobK

If you do `debug(OpenMx:::mxGenerateData)` and then call `mxCompare(..., boot=TRUE)` then you will see that `mxGenerateData` is used (for the null model). Maybe I misunderstood something when we discussed the implementation. Rob, are you sure it can be made to work without `mxData`?
Replied on Thu, 06/24/2021 - 16:08
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by jpritikin

If you do debug(OpenMx:::mxGenerateData) and then call mxCompare(..., boot=TRUE) then you will see that mxGenerateData is used (for the null model). Maybe I misunderstood something when we discussed the implementation. Rob, are you sure it can be made to work without mxData?

I'm looking over the code for `mxGenerateData()`, and yeah, I'm pretty sure that it can work without the models containing any MxData objects. It doesn't need raw data, just the expected sufficient statistics under the null model. Of course, what exactly `mxGenerateData()` does depends upon the definition of method `genericGenerateData()` for the relevant expectation class...which I notice only has definitions for MxExpectationNormal, MxExpectationLISREL, MxExpectationRAM, and MxExpectationStateSpace. There is no definition for MxExpectationMixture.