机器算法验证 - 重复测量 anova 模型和混合模型之间的等价性：lmer vs lme，以及复合对称性 - 吾爱随笔录

重复测量 anova 模型和混合模型之间的等价性：lmer vs lme，以及复合对称性

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

2022-03-09 11:14:29

aov我在重复测量模型和lmer混合模型之间获得等效结果时遇到了一些麻烦。

我的数据和脚本如下所示

data=read.csv("https://www.dropbox.com/s/zgle45tpyv5t781/fitness.csv?dl=1")
data$id=factor(data$id)
data
   id  FITNESS      TEST PULSE
1   1  pilates   CYCLING    91
2   2  pilates   CYCLING    82
3   3  pilates   CYCLING    65
4   4  pilates   CYCLING    90
5   5  pilates   CYCLING    79
6   6  pilates   CYCLING    84
7   7 aerobics   CYCLING    84
8   8 aerobics   CYCLING    77
9   9 aerobics   CYCLING    71
10 10 aerobics   CYCLING    91
11 11 aerobics   CYCLING    72
12 12 aerobics   CYCLING    93
13 13    zumba   CYCLING    63
14 14    zumba   CYCLING    87
15 15    zumba   CYCLING    67
16 16    zumba   CYCLING    98
17 17    zumba   CYCLING    63
18 18    zumba   CYCLING    72
19  1  pilates   JOGGING   136
20  2  pilates   JOGGING   119
21  3  pilates   JOGGING   126
22  4  pilates   JOGGING   108
23  5  pilates   JOGGING   122
24  6  pilates   JOGGING   101
25  7 aerobics   JOGGING   116
26  8 aerobics   JOGGING   142
27  9 aerobics   JOGGING   137
28 10 aerobics   JOGGING   134
29 11 aerobics   JOGGING   131
30 12 aerobics   JOGGING   120
31 13    zumba   JOGGING    99
32 14    zumba   JOGGING    99
33 15    zumba   JOGGING    98
34 16    zumba   JOGGING    99
35 17    zumba   JOGGING    87
36 18    zumba   JOGGING    89
37  1  pilates SPRINTING   179
38  2  pilates SPRINTING   195
39  3  pilates SPRINTING   188
40  4  pilates SPRINTING   189
41  5  pilates SPRINTING   173
42  6  pilates SPRINTING   193
43  7 aerobics SPRINTING   184
44  8 aerobics SPRINTING   179
45  9 aerobics SPRINTING   179
46 10 aerobics SPRINTING   174
47 11 aerobics SPRINTING   164
48 12 aerobics SPRINTING   182
49 13    zumba SPRINTING   111
50 14    zumba SPRINTING   103
51 15    zumba SPRINTING   113
52 16    zumba SPRINTING   118
53 17    zumba SPRINTING   127
54 18    zumba SPRINTING   113

基本上，3 x 6 名受试者 ( id) 分别接受了三种不同的FITNESS锻炼计划，并PULSE在执行三种不同类型的耐力后对其进行测量TEST。

然后我安装了以下aov模型：

library(afex)
library(car)
set_sum_contrasts()
fit1 = aov(PULSE ~ FITNESS*TEST + Error(id/TEST),data=data)
summary(fit1)
Error: id
          Df Sum Sq Mean Sq F value   Pr(>F)    
FITNESS    2  14194    7097   115.1 7.92e-10 ***
Residuals 15    925      62                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: id:TEST
             Df Sum Sq Mean Sq F value   Pr(>F)    
TEST          2  57459   28729   253.7  < 2e-16 ***
FITNESS:TEST  4   8200    2050    18.1 1.16e-07 ***
Residuals    30   3397     113                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我使用得到的结果

set_sum_contrasts()
fit2=aov.car(PULSE ~ FITNESS*TEST+Error(id/TEST),data=data,type=3,return="Anova")
summary(fit2)

与此相同。

混合模型运行 usingnlme给出直接等效的结果，例如使用lme：

library(lmerTest)    
lme1=lme(PULSE ~ FITNESS*TEST, random=~1|id, correlation=corCompSymm(form=~1|id),data=data)
anova(lme1)
             numDF denDF   F-value p-value
(Intercept)      1    30 12136.126  <.0001
FITNESS          2    15   115.127  <.0001
TEST             2    30   253.694  <.0001
FITNESS:TEST     4    30    18.103  <.0001


summary(lme1)
Linear mixed-effects model fit by REML
 Data: data 
       AIC      BIC    logLik
  371.5375 393.2175 -173.7688

Random effects:
 Formula: ~1 | id
        (Intercept) Residual
StdDev:    1.699959 9.651662

Correlation Structure: Compound symmetry
 Formula: ~1 | id 
 Parameter estimate(s):
       Rho 
-0.2156615 
Fixed effects: PULSE ~ FITNESS * TEST 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   81.33333  4.000926 30 20.328628  0.0000
FITNESSpilates                 0.50000  5.658164 15  0.088368  0.9308
FITNESSzumba                  -6.33333  5.658164 15 -1.119327  0.2806
TESTJOGGING                   48.66667  6.143952 30  7.921069  0.0000
TESTSPRINTING                 95.66667  6.143952 30 15.570868  0.0000
FITNESSpilates:TESTJOGGING   -11.83333  8.688861 30 -1.361897  0.1834
FITNESSzumba:TESTJOGGING     -28.50000  8.688861 30 -3.280062  0.0026
FITNESSpilates:TESTSPRINTING   8.66667  8.688861 30  0.997446  0.3265
FITNESSzumba:TESTSPRINTING   -56.50000  8.688861 30 -6.502579  0.0000

或使用gls：

library(lmerTest)    
gls1=gls(PULSE ~ FITNESS*TEST, correlation=corCompSymm(form=~1|id),data=data)
anova(gls1)

但是，我使用lme4's获得的结果lmer是不同的：

set_sum_contrasts()
fit3=lmer(PULSE ~ FITNESS*TEST+(1|id),data=data)
summary(fit3)
Linear mixed model fit by REML ['lmerMod']
Formula: PULSE ~ FITNESS * TEST + (1 | id)
   Data: data

REML criterion at convergence: 362.4

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept)  0.00    0.0     
 Residual             96.04    9.8     
...

Anova(fit3,test.statistic="F",type=3)
Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)

Response: PULSE
                    F Df Df.res    Pr(>F)    
(Intercept)  7789.360  1     15 < 2.2e-16 ***
FITNESS        73.892  2     15 1.712e-08 ***
TEST          299.127  2     30 < 2.2e-16 ***
FITNESS:TEST   21.345  4     30 2.030e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lmer有人认为我在模型上做错了什么吗？或者差异来自哪里？它是否必须与lmer不允许负类内关联或类似的事情做任何事情？但是，鉴于nlme'gls和lme确实返回了正确的结果，我想知道gls和有什么不同lme？是不是该选项correlation=corCompSymm(form=~1|id)导致他们直接估计类内相关性，它可以是正的或负的，而lmer估计一个不能为负的方差分量（在这种情况下最终被估计为零）？

1个回答

正如 Ben Bolker 在评论中已经提到的那样，问题正如您所怀疑的那样：lmer()模型被绊倒了，因为它试图估计方差分量模型，而方差分量估计值被限制为非负数。我将尝试做的是对导致这种情况的数据集有什么直观的理解，以及为什么这会导致方差分量模型出现问题。

这是您的数据集的图。白点是实际观察值，黑点是主体均值。

在此处输入图像描述

为了使事情更简单，但不改变问题的精神，我将减去固定效应（即FITNESS和TEST效应，以及大均值）并将残差数据作为单向随机效应问题处理. 所以这是新数据集的样子：

在此处输入图像描述

仔细观察这个情节中的模式。想想从同一主题获得的观察结果与从不同主题获得的观察结果有何不同。具体来说，请注意以下模式：由于受试者的一个观察值高于（或低于）受试者平均值，因此来自该受试者的其他观察值往往位于受试者平均值的另一侧。并且观察值与主体均值的距离越远，其他观测值往往与相反侧的主体均值越远。这表明负的类内相关性。平均而言，从同一主题获得的两次观察实际上往往比从数据集中随机抽取的两次观察更不相似。

考虑这种模式的另一种方法是根据主体间和主体内方差的相对大小。与受试者之间的差异相比，受试者内的差异似乎要大得多。当然，我们预计这会在某种程度上发生。毕竟，主体内方差是基于单个数据点的变化，而主体间方差是基于单个数据点（即主体均值）的变化，我们知道随着被平均的事物数量的增加，平均值将趋于减少。但在这个数据集中，差异非常显着：有办法主体内的变化多于主体间的变化。实际上，这种差异正是出现负类内相关性的原因。

好的，问题就在这里。方差分量模型假设每个数据点是一个主题效应和一个误差的总和：，其中是第个主题的效应。因此，让我们考虑一下如果主体效应的方差真正为 0 会发生什么——换句话说，如果真正的主体间方差分量为 0。给定在此模型下生成的实际数据集，如果我们要计算样本均值对于每个受试者的观察数据，这些样本均值仍然会有一些非零方差，但它们只会反映误差方差，而不是任何“真实”的受试者方差（因为我们假设没有）。 $y_{ij}=u_j+e_{ij}$ $u_j$ $j$

那么我们期望这些主题手段有多大的可变性呢？好吧，基本上每个估计的主体效应都是一个平均值，我们知道平均值方差的公式：，其中是被平均的事物的数量。现在让我们将此公式应用于您的数据集，看看如果真实的受试者间方差分量恰好为 0，我们期望在估计的受试者效果中看到多少方差。 $\text{var}(\bar{X})=\text{var}(X_i)/n$ $n$

主体内方差为，每个主体效应计算为 3 个观测值的平均值。因此，主题均值的预期标准差——假设真实的主题间方差为 0——计算结果约为。现在将其与我们实际观察到的主题均值的标准差进行比较：！当我们假设受试者间方差为 0 时，观察到的变化远小于预期的变化。对于方差分量模型，可以预期观察到的变异与我们实际观察到的一样低的唯一方法是，如果真实的主体间方差在某种程度上为负 $348$ $10.8$ $4.3$ . 这就是问题所在。数据表明存在某种负方差分量，但软件（明智地）不允许对方差分量进行负估计，因为方差实际上永远不会是负数。您拟合的其他模型通过直接估计类内相关性而不是假设简单的方差分量模型来避免这个问题。

如果您想了解如何实际获得数据集隐含的负方差分量估计，您可以使用我在我的另一个最近的答案中R说明的过程（附带代码）。该过程并非完全微不足道，但也不是太难（对于像这样的平衡设计）。

其它你可能感兴趣的问题

上一篇如何检验具有不同样本量的两组的方差同质性？下一篇如何评估生存函数的拟合优度