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)) )  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.