Hello! I’m trying to follow along with the “Multilevel Models in Julia” tutorial by Yichi Zhang (I can’t include a link, but it’s publicly available) using my own data. In this exercise, you implement a random intercepts model (mm1), one with only the level-2 predictor (mm2), one with only the level-1 predictor (mm3), and then one with a cross-level interaction (mm4). I was able to do all of that successfully, but when I try to call MixedModels.likelihoodratiotest(mm1,mm2,mm3,mm4) I get the following error:
Through trial and error, I’ve found that the likelihood ratio test will run successfully if I exclude mm3, which (as in the tutorial) has just the level-1 variables. I’m a little confused as to why it’s causing an error when mm4 has those variables in it as well. Does anyone have an inkling as to what I might be doing wrong here?
Hi and welcome!
Here are the 4 models from the blogpost:
- mm1:
happy ~ (1|country)
- mm2:
happy ~ gm_GNP + (1|country)
- mm3:
happy ~ income_cmc + income_mean + (income_cmc|country)
- mm4:
happy ~ income_cmc * gm_GNP + income_mean + (income_cmc|country)
A model is nested in a larger model if the smaller model can be obtained by imposing constraints on the larger model’s parameters. For nested models, the likelihood-ratio test statistic is asymptotically chi-squared, with degrees of freedom equal to the difference in the number of estimated parameters, under the usual regularity conditions (which is what MixedModels.likelihoodratiotest uses).
Let’s check the models:
mm1 with mm2? Yes, if you constrain gm_GNP coefficient to be 0 they are the same model.
mm2 with mm4? Recall A * B in formula syntax means main effects with interaction, and is equivalent to A + B + A:B. Thus, if coefficient for income_cmc, income_cmc:gm_GNP, and income_mean is set to 0 we get the mean structure of mm2. The covariance structure can be made the same by constraining Var(income_cmc slope) = 0 and Cov(intercept, slope) = 0.
Technical details
Note that constraining it to the boundary of the parameter space like we did with the variance doesn’t satisfy the usual regularity conditions, so the asymptotic chi-square with 5 dof is wrong here (which MixedModels uses). Instead, the correct asymptotic distribution of the LRT is a mixture of chi-squared distributions (Self & Liang 1987; Stram & Lee 1994) - simulation methods, such as parametic bootstrap, may perform better in practice with finite samples (Crainiceanu & Ruppert 2004). The naive way that MixedModels uses will be conservative.
What about mm2 with mm3? Since there is no gm_GNP variable in mm3 they are not nested.
MixedModels.likelihoodratiotest checks if models are nested before performing an otherwise incorrect likelihood ratio test, which is why you see the error. Since that blog post is rather old (2021), it didn’t check for that back then (probably?).
