在最近提出的关于线性混合效应模型的问题中,有人告诉我,不应使用似然比检验在具有不同随机效应结构的模型之间进行比较。到目前为止,我一直在装有 REML 的嵌套模型上使用这种方法,在该模型中,固定效应保持不变,以此作为找到最佳随机效应结构的一种方式。我的方法基于 Alain Zuur(2009 年)第 5 章为生态学家广泛使用的关于统计建模的书“混合效应模型和生态学扩展”。Pinheiro & Bates 的另一本关于 LME 的书也支持这种方法(2000) 即第 83 页。
我想就这是否确实是一种不健全的方法寻求进一步的建议,如果是,请在 R 中找到一个更强大的可行替代方案。
我在下面给出了两个嵌套模型的示例(使用 R 中的 lme() 函数创建),以及我目前如何将它们与 LRT 或 AIC 进行比较:
# 模型 1:随机截距模型 # > M1 = lme(dtim ~ dd, random = ~1 | fInd,数据=df,方法=“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 表明在模型中添加随机斜率项是一项重大改进。AIC 也较低。任何关于这是否是一种稳健方法的建议将不胜感激。