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.