Justifying the Need for Mixed Effects Models (aka. LME, MLM, etc.)

Question

Firstly, I am not an expert in using multi-level modelling (MLM), and I have read this and this questions, however, my question is slightly different in the sense that method 2 below is not mentioned.

Since multi-level modelling is rather complex, I want to justify the need for it in the first place. In that regard, I know of two methods:

1) Assessing whether there is enough and significant variation across items (aka. contexts):

This method is mentioned in the book Discovering Statistics Using R, section 19.6.6. It implies comparing a baseline intercept-only generalized least squares fit by maximum likelihood to another linear mixed-effects model fit by maximum likelihood where intercepts are allowed to vary across items. If the fit significantly improves, this warrants using MLM.

My example of the two models where R is the response/outcome variable:

M1 = nlme::gls(R ~ 1, data = univariate_data, method = "ML")
M2 = nlme::lme(R ~ 1, data = univariate_data, method = "ML", random = ~1|item_id)

The ANOVA comparison:

##                    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## M1                     1  2 9181.778 9191.491 -4588.889                        
## M2                     2  3 9170.908 9185.477 -4582.454 1 vs 2 12.87025   3e-04

From the tests, we see that after addressing the variability in our items/contexts, there is a significant improvement in the log-likelihood by 12.87 at the expense of 1 degree of freedom, so: $\chi^2(1) = 12.87, p = .0003$. This necessitates using MLM.

2) Comparing Unconditional LME Models:

I read this online but I don't recall where. The two unconditional LME models are compared to each other, and if allowing the intercepts to vary across items (contexts) does significantly improve the fit, then using MLM is asserted.

My example of the two models where R is the response/outcome variable:

MN1 = lmer(R ~ 1 + (1 | subject_id), data = univariate_data, REML = FALSE,
           control = lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))

MN2 = lmer(R ~ 1 + (1 | subject_id) + (1 | item_id), data = univariate_data, REML = FALSE,
           control = lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))

The ANOVA comparison:

##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## MN1            3 7096.9 7110.7 -3545.5   7090.9                     
## MN2            4 7096.7 7115.1 -3544.4   7088.7 2.1966  1     0.1383

The fit, as you see, is not significantly different between the two LME models.

My conundrum arises from having inconsistent results: method 1 justified MLM but method 2 doesn't. How can we interpret this discrepancy? and what method is more robust in order to study the feasibility of MLM?

Note: in a previous question I came to know that visual inspection alone is a weak approach for studying the feasibility of MLM.

Answer 1

First of all, I agree with @EdM, you should first think if multilevel models in substantive terms instead of only looking for a statistical significance in model comparisons. But you may ask: how to do so? Well, asking yourself a few questions about your research hypothesis might help. Is it important to your research:

• 1. Assessing the effect of more than one second-level predictor?
• 2. Assessing the interaction between first and second level predictors (or third and second, etc.)?
• 3. The decomposition of between-effects and within-effects of first level predictors?

If you answered yes to any of these questions, then a multilevel model may be your best, if not only choice. I guess this reasoning comes from Gelman & Hill classic book on multilevel models. You can probably find it in other links as well. This classic paper Enders and Tofighi (2007) may help you with insights on these matters, and also is useful to learn how to avoid criticism for correlation problems within multilevel models (MLM). You also find similar rationale in this and this papers. I explain it a bit better below.

Second, you might be also interested in the intra-class correlation (ICC). That is, you want to know how much of your dependent variable variance is due to differences between or within second and third level groups. A high ICC is sometimes presented as an evidence that you need a MLM. It is in fact an evidence that you have a highly clustered data, which violates standard OLS regression assumptions (this is almost always the case with repeated measures/within-subject research designs). But you may correct it with fixed effects as well, as you may see in the aforementioned papers.

The problem with fixed effects is that only one variable will suck up all the variance from the higher levels groups. That is why if you are interested in more than one second level predictor, a MLM may be your best choice. In addition, the ICC may be an interesting measure itself. At the same time it provides some support for a multilevel model, and it is also an easy to understand and useful substantive interpretation of the multilevel structure of your data.

Third, although research hypothesis and substantive reasoning should come first in model choice, model fit statistics also matter. If your new model has a worse fit to your data than your previous model, you probably need to go back, or change something else in the new one.

Fourth, your question:

"method 1 justified MLM but method 2 doesn't. How can we interpret this discrepancy? and what method is more robust in order to study the feasibility of MLM?"

Well, there is not any discrepancy between the results. The problem is that in "method 1" you are comparing a one-level model with a two-level model. In "method 2" you are comparing a two-level model with a three-level model. That is why you got different results. You are comparing different stuff. The two methods are pretty much the same: a chi-square test. But what you put inside them that was different.

As a matter of fact, you are comparing four different models with only two pairwise comparisons. In M1 you have only one level. In M2 you have a two-level model, and your second level is item. In MN1 you have a two level-model, but your second level is subject. In MN2 you have a three-level model and your second level is subject, and your third level is item. You cannot say which of these four models is the best, at least not if you look only to the p-value, because you did not compare them all.

However, if you forgot a bit about the p-value, and take a look to other statistics, specially the Akaike Information Criteria (AIC) and the Bayesian Information Criteria (BIC), which are currently the standard for assessing model fit of multilevel models, you may see that your best model is between MN1 and MN2.

If you take a close look in the second test: You will see that your p-value is not super high. You also see that the AIC and the log-likelihood ratio test shows that your three-level model has a slight better fit than the two-level model. The BIC and Chi-square show the opposite. How can your p-value be higher than 0.05? Well, your p-value is based solely on the Chi-square test, which is probably the less used test (of the four shown) to compare MLM. But it is really a close call. So, you have a statistical tie, and you can decide based on your substantive questions.

Fifth, you may also consider that multilevel models have some assumptions and mathematical requirements. @EdM may be right when he says that MLM work with partial pooling and that you don’t need to worry so much about few first level observations per group. I also learned this way. However, this assumption is questioned by this paper. I never follow this last paper sample size requirements, and in fact, I don't know any multilevel model in social sciences that do. But it is more or less consolidated in the field that you need at least 50 higher level units/groups, as you can see here. If you have fewer than that, fixed effects might get you better estimates.

Answer 2

What you are testing for is the statistical "significance" of the random effect terms. As is so often the case in statistics, that can be a good deal different from their "importance." Note this part of the answer to one of the questions that you have read:

Random effects are typically included to account for the correlation of the measurements within a group/cluster.

If your experimental design involves correlations among observations, potentially invalidating the assumption of independent observations, those correlations need to be taken into account. If there are more than a few individuals/groups/clusters, modeling them with random effects provides a well documented way to do do. You might find that the random effects don't add anything "significant" to your model when you perform your tests, but your audience will expect that it is "important" to take such correlations into account.

In multi-level modeling you might even choose to incorporate predictors like gender, typically modeled as fixed effects, as random effects. Discussion on this page and this page cover circumstances in which this can make sense. If you are using frequentist tools like lmer() for multi-level modeling, this can provide an advantage. The random effects are modeled by partial pooling among all groups rather than estimating for each group separately. As this answer notes:

Partial pooling means that, if you have few data points in a group, the group's effect estimate will be based partially on the more abundant data from other groups. This can be a nice compromise between estimating an effect by completely pooling all groups, which masks group-level variation, and estimating an effect for all groups completely separately, which could give poor estimates for low-sample groups.

If you want the advantages of partial pooling, it would be "important" to model with random effects even if they don't turn out to be "significant."

So worry less about statistically justifying a choice of random-effect modeling; think more about whether such modeling will help strengthen your analysis.