Comparing between random effects structures in a linear mixed-effects model

Question

During a recently asked question about linear mixed-effects models I was told that one should not compare between models with different random effects structures using likelihood ratio tests. Up until now, I had used this approach on nested models fitted with REML in which the fixed effects were kept constant as a way to find the optimal random effects structure. My method was based on a widely used book on statistical modelling for ecologists "Mixed effects models and extensions in ecology with R" written by Alain Zuur (2009) chapter 5. This approach was also backed up in another book on LMEs by Pinheiro & Bates (2000) i.e. pg 83.

I would like to seek further advice on whether this is indeed an unsound method, and if so, find a workable alternative within R that is more robust.

I give examples of two nested models below (created using the lme() function in R) and how I would currently compare between them with LRTs or AICs:
# Model 1: Random intercept model # > M1 = lme(dtim ~ dd, random= ~1 | fInd, data=df, method="REML")

Linear mixed-effects model fit by REML
 Data: df 
       AIC      BIC    logLik
  47344.74 47373.58 -23668.37

Random effects:
 Formula: ~1 | fInd
        (Intercept) Residual
StdDev:   0.5244626 2.574662

Fixed effects: dtim ~ dd
                 Value  Std.Error   DF   t-value p-value
(Intercept) -0.8681514 0.17048746 9988  -5.09217       0
    dd       2.2424996 0.01260611 9988 177.88982       0
 Correlation: 
       (Intr)
   dd -0.203

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-5.7610884 -0.4620287 -0.1732839  0.2395293 13.0981698 

Number of Observations: 10000
Number of Groups: 11


# Model 2: random intercept and slope model
# > M2 = lme(dtim ~ dd, data=df, random= ~1 + dd|fInd, method="REML")

Linear mixed-effects model fit by REML
 Data: df 
       AIC      BIC    logLik
  47041.82 47085.08 -23514.91

Random effects:
 Formula: ~1 + dd | fInd
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 0.4860448 (Intr)
dd          0.3231004 -0.687
Residual    2.5314343       

Fixed effects: dtim ~ dd 
                 Value  Std.Error   DF   t-value p-value
(Intercept) -0.5568345 0.15839434 9988 -3.515495   4e-04
 dd          2.0912224 0.09974746 9988 20.965168   0e+00
 Correlation: 
       (Intr)
 dd   -0.676

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-4.6988351 -0.4460439 -0.1848166  0.2296023 12.9419866 

Number of Observations: 10000
Number of Groups: 11 

# Compare the two models using LRTs
> anova(M1,M2)
   Model df      AIC      BIC    logLik   Test  L.Ratio p-value
M1     1  4 47344.74 47373.58 -23668.37                        
M2     2  6 47041.82 47085.08 -23514.91 1 vs 2 306.9191  <.0001

# L ratio test statistic: to get correct p-value from L ratio test I would then use the below formula (due to testing the boundary effect)
# 0.5 * (  (1 - pchisq(L.ratio, 1)) +  (1 - pchisq(L.ratio, 2))  )
> 0.5 * (  (1 - pchisq(306.9191, 1)) +  (1 - pchisq(306.9191, 2))  )
[1] 0

L.Ratio suggests that adding a random slope term to the model is a significant improvement. The AIC is also lower. Any advice on whether this is a robust method would be appreciated.

Answer 1

I was the one suggesting this to you; as I mentioned to my comments there though: "Apologies for being misleading most of my comment regarded selection (on) $X$ not $Z$". By that I mean that I was referring mostly to the fixed effects rather than the random effects structure.

Yes, you can use LRT if you have the same $X$ while using a model fitted by REML. You should be able to use AIC in these cases with caution. This is because it is not obvious how to define the degrees of freedom associated with a specific random effect. You should not use AIC's "vanilla" version directly. Please look at Greven and Kneib, 2010 regarding this; they present a corrected cAIC. They also provide an R package implementing the corrected cAIC they outline.

AIC and LRT are asymptotic tests but things tend to get hairy when you need to estimate parameters that might be close to the boundary of your sample space (ie. when you are testing for variances being close to $0$. In that case you actually want a mixture of $\chi^2$-distributions. A relevant reference of that is Lindquist et al., 2012. To that extent Morell, 1999 can also help if a theoretical justification regarding the use of ReML.

You inquired for a "robust method" to select your random effects structure; on first instance, bootstrap your sample. Use parametric bootstrap to evaluate the asymptotic behavior of your model. Please see the comments mentioned in glmm.wikidot regarding whether a random effect is significant. As mentioned to you in my earlier comment I would be extremely cautious to start model-selection on $Z$; I prefer to "treat it as given" based on my research question. Otherwise I simply cherry-pick my error structure trying to "squeeze more significance out of the remaining terms" [glmm.wikidot].

To recap: using LRT is not "unsound"; it though prone to the limitations of LRTs regarding their asymptotic behavior. There are a number of references on how to provide a remedy. The easiest thing for you at this point would be to simply use RLRsim at first instance. It is based on another piece of work of Greven, Scheipl et al., 2008.

Answer 2

You certainly can't use AIC, BIC or similar criteria that contain an explicit penalty term computed based on the number of parameters in the model. As I pointed out in this topic, the effective number of parameters associated with random effects is unknown. I wasn't sure I was right when I posted that question, but no one challenged me.

Likewise, to compute a p-value based on the LR statistic, one has to know the difference in the number of parameters between the models. I have an ominous feeling that, just like for AIC, that difference has to be in effective terms as opposed to nominal terms.