机器算法验证 - 带有 Error() 的 aov() 是否与 R 中 lme4 包的 lmer() 相同？ - 吾爱随笔录

带有 Error() 的 aov() 是否与 R 中 lme4 包的 lmer() 相同？

机器算法验证 r 方差分析混合模式重复测量 lme4-nlme

2022-03-19 07:26:12

以下2个模型真的一样吗？

> library(lme4)
> library(lmerTest)
> lmod = lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy)    
> summary(lmod)
Linear mixed model fit by REML 
t-tests use  Satterthwaite approximations to degrees of freedom ['merModLmerTest']
Formula: Reaction ~ Days + (Days | Subject)
   Data: sleepstudy

REML criterion at convergence: 1743.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9536 -0.4634  0.0231  0.4634  5.1793 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.09   24.740       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error      df t value             Pr(>|t|)    
(Intercept)  251.405      6.825  17.000  36.838 < 0.0000000000000002 ***
Days          10.467      1.546  17.000   6.771           0.00000326 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
Days -0.138

和：

> laov = aov(Reaction ~ Days + Error(Subject/Days), data=sleepstudy)   
> summary(laov)

Error: Subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 17 250618   14742               

Error: Subject:Days
          Df Sum Sq Mean Sq F value     Pr(>F)      
Days       1 162703  162703   45.85 0.00000326 ***  
Residuals 17  60322    3548                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
           Df Sum Sq Mean Sq F value Pr(>F)
Residuals 144  94312   654.9

两者都显示了 Days 变量的相似 P 值。两种方法有什么区别？

1个回答

这两个模型给出相同的 $p$ 值，并且 $\sqrt{F}=\sqrt{45.85}=6.77=t$ ，是模型确实相同的强有力的先验证据。两者最终都适合模型

\begin{aligned} μ_{i t} & = β_{0} + b_{0, i} + (β_{1} + b_{1, i}) t \\ Y_{i} & \sim Normal (μ_{i t}, σ^{2}) \\ b & \sim MVN (0, Σ) \end{aligned}

$\begin{split} \mu_{it} & = \beta_0 + b_{0,i} + (\beta_1 + b_{1,i}) t \\ Y_i & \sim \textrm{Normal}(\mu_{it}, \sigma^2) \\ \boldsymbol{b} & \sim \textrm{MVN}(\boldsymbol{0},\Sigma) \end{split}$

一般来说，aov()它更简单、更快，但只适用于更简单的模型子集。

请注意，虽然nlme::lme也可以拟合此模型，并获得相同的统计量，但在这种情况下，其猜测正确分母自由度的启发式方法是不正确的（它得到 161 而不是 17）： $F$

library(nlme)
anova(lme(Reaction ~ Days,
                random = ~Days|Subject, sleepstudy))
##             numDF denDF   F-value p-value
## (Intercept)     1   161 1454.0766  <.0001
## Days            1   161   45.8534  <.0001

@amoeba 正确地指出

一个非常重要的警告是，这两个模型只是因为Days被视为数值变量而重合。如果它是一个分类因素，那么aovwithError(Subject/Days)和lmerwith(Days|Subject)将是两个不同的模型。要获得 [在这种情况下] 的 [the]lmer等价物，aov可能需要使用(1|Subject)+(1|Subject:Days)[or equivalent (1|Subject/Days)]

其它你可能感兴趣的问题

上一篇概率及其补码的乘积的术语是什么？下一篇基测度在指数族中的作用