比较线性混合效应模型中的随机效应结构

机器算法验证 混合模式 随机效应模型 aic 似然比
2022-03-27 06:58:05

最近提出的关于线性混合效应模型的问题中,有人告诉我,不应使用似然比检验在具有不同随机效应结构的模型之间进行比较。到目前为止,我一直在装有 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 也较低。任何关于这是否是一种稳健方法的建议将不胜感激。

2个回答

我是向你建议这个的人;正如我在那里的评论中提到的那样:“抱歉误导了我的大部分评论认为选择(在)而不是XZ ”。我的意思是我主要指的是固定效应而不是随机效应结构。

是的,如果您在使用 REML 拟合的模型时 ,则可以使用 LRT 。在这些情况下,您应该能够谨慎使用 AIC。这是因为如何定义与特定随机效应相关的自由度并不明显。您不应该直接使用 AIC 的“香草”版本。参阅 Greven 和 Kneib,2010 年的相关内容;他们提出了一个修正的cAIC。他们还提供了一个 R 包,用于实现他们概述的更正 cAIC。X

AIC 和 LRT 是渐近测试,但是当您需要估计可能接近样本空间边界的参数时(即,当您测试接近的方差时,事情往往会变得棘手。在这种情况下,您实际上想要一个 -分布的混合。相关的参考文献是Lindquist et al., 2012。在这方面,Morell, 1999也可以提供有关使用 ReML 的理论依据。0χ2

您询问了一种“稳健的方法”来选择您的随机效应结构;在第一个实例中,引导您的样本。使用参数引导来评估模型的渐近行为。请参阅glmm.wikidot中提到的关于随机效应是否显着的评论。正如我之前的评论中提到的那样,我会非常谨慎地开始上的模型选择;根据我的研究问题,我更喜欢“把它当作给定的”。否则,我只是简单地挑选我的错误结构,试图“从剩余的术语中挤出更多的意义”[ glmm.wikidot ]。Z

回顾一下:使用轻轨并非“不合理”;它虽然容易受到 LRT 关于其渐近行为的限制。有许多关于如何提供补救措施的参考资料。此时对您来说最简单的事情就是在第一次使用RLRsim它基于 Greven, Scheipl et al., 2008的另一件作品。

您当然不能使用 AIC、BIC 或包含基于模型中参数数量计算的显式惩罚项的类似标准。正如我在本主题中指出的那样,与随机效应相关的参数的有效数量是未知的。当我发布这个问题时,我不确定我是否正确,但没有人质疑我。

同样,要根据 LR 统计量计算 p 值,必须知道模型之间参数数量的差异。我有一种不祥的感觉,就像 AIC 一样,这种差异必须是有效的,而不是名义上的。