机器算法验证 - 如何在线性混合模型中选择随机效应和固定效应结构？ - 吾爱随笔录

如何在线性混合模型中选择随机效应和固定效应结构？

机器算法验证混合模式重复测量模型选择 lme4-nlme 似然比

2022-01-31 06:00:11

从受试者设计的双向考虑以下数据：

df <- "http://personality-project.org/r/datasets/R.appendix4.data"
df <- read.table(df,header=T)
head(df)

Observation Subject Task Valence Recall
1           1     Jim Free     Neg      8
2           2     Jim Free     Neu      9
3           3     Jim Free     Pos      5
4           4     Jim Cued     Neg      7
5           5     Jim Cued     Neu      9
6           6     Jim Cued     Pos     10

我想使用混合线性模型对此进行分析。考虑到所有可能的固定效应和随机效应，有多种可能的模型：

# different fixed effects with random-intercept
a0 <- lmer(Recall~1 + (1|Subject), REML=F,df)
a1 <- lmer(Recall~Task + (1|Subject), REML=F,df)
a2 <- lmer(Recall~Valence + (1|Subject), REML=F,df)
a3 <- lmer(Recall~Task+Valence + (1|Subject), REML=F,df)
a4 <- lmer(Recall~Task*Valence + (1|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope
b0 <- lmer(Recall~1 + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b1 <- lmer(Recall~Task + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b2 <- lmer(Recall~Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b3 <- lmer(Recall~Task+Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b4 <- lmer(Recall~Task*Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope including variance-covariance matrix
c0 <- lmer(Recall~1 + (1 + Valence + Task|Subject), REML=F,df)
c1 <- lmer(Recall~Task + (1 + Valence + Task|Subject), REML=F,df)
c2 <- lmer(Recall~Valence + (1 + Valence + Task|Subject), REML=F,df)
c3 <- lmer(Recall~Task+Valence + (1 + Valence + Task|Subject), REML=F,df)
c4 <- lmer(Recall~Task*Valence + (1 + Valence + Task|Subject), REML=F,df)

在这种情况下选择最佳拟合模型的推荐方法是什么？使用对数似然比检验时，推荐的程序是什么？向上（从零模型到最复杂的模型）或向下（从最复杂的模型到零模型）生成模型？逐步包含或排除？还是建议将所有模型进行一次对数似然比检验并选择 p 值最低的模型？如何比较未嵌套的模型？
是否建议先找到适当的固定效应结构，然后再找到适当的随机效应结构，或者反过来（我已经找到了这两个选项的参考资料......）？
报告结果的推荐方式是什么？报告来自对数似然比检验的 p 值，将完整混合模型（有问题的影响）与简化模型（没有问题的影响）进行比较。还是使用对数似然比检验找到最佳拟合模型，然后使用 lmerTest 报告最佳拟合模型中效果的 p 值是否更好？

2个回答

我不确定这是否真的有一个规范的答案，但我会试一试。

在这种情况下选择最佳拟合模型的推荐方法是什么？使用对数似然比检验时，推荐的程序是什么？向上（从零模型到最复杂的模型）或向下（从最复杂的模型到零模型）生成模型？逐步包含或排除？还是建议将所有模型进行一次对数似然比检验并选择 p 值最低的模型？如何比较未嵌套的模型？

这取决于你的目标是什么。

一般来说，您应该非常非常小心地选择模型（参见例如这个答案，或者这个帖子，或者只是谷歌“Harrell stepwise”......）。
如果您有兴趣让您的 p 值有意义（即您正在进行验证性假设检验），则不应进行模型选择。但是：如果您对模型的非焦点部分进行模型选择，例如，如果您的主要兴趣是推断固定效应，则对随机效应进行模型选择对我来说不是很清楚。
没有“将所有模型放在一个似然比测试中”这样的事情——似然比测试是一个成对的过程。如果您想对随机效应进行模型选择（例如），我可能会推荐使用本示例中的信息标准的“一次性”方法——至少避免了逐步方法的一些问题（但不是模型选择更普遍）。
巴尔等人。2013 年“保持最大”记忆和语言杂志(doi:10.1016/j.jml.2012.11.001) 建议使用最大模型（仅）。
Shravan Vasishth 不同意，他认为除非数据集非常大（并且信噪比很高），否则此类模型将动力不足并因此存在问题
另一种合理辩护的方法是拟合一个大但合理的模型，然后，如果拟合是奇异的，则删除项，直到不再存在
有一些注意事项（在GLMM FAQ中列举），您可以使用信息标准来比较具有不同随机效应的非嵌套模型（尽管 Brian Ripley 不同意：请参见此处第 6 页的底部）

是否建议先找到适当的固定效应结构，然后再找到适当的随机效应结构，或者反过来（我已经找到了这两个选项的参考资料......）？

我想没有人知道。更一般地参见关于模型选择的先前答案。如果您可以足够清楚地定义您的目标（很少有人这样做），那么这个问题可能是可以回答的。如果您对这两个选项都有参考，那么编辑您的问题以包含它们会很有用......（对于它的价值，这个例子（上面已经引用过）使用信息标准来选择随机效应部分，然后避开选择模型的固定效应部分。

报告结果的推荐方式是什么？报告来自对数似然比检验的 p 值，将完整混合模型（有问题的影响）与简化模型（没有问题的影响）进行比较。还是使用对数似然比检验找到最佳拟合模型，然后使用 lmerTest 报告最佳拟合模型中效果的 p 值是否更好？

这是（唉）另一个难题。如果您报告所报告的边际效应lmerTest，则必须担心边际性（例如，当模型中存在逐-交互作用时，对A和的主要影响的估计是否有意义）；这是一个巨大的蠕虫罐头，但如果您按照. 平衡的设计也有一点帮助。如果您真的想解决所有这些裂缝，我想我会推荐，它会为您提供类似于的输出，但会尝试处理这些问题。BABcontrasts="sum"afex::mixed()afex::mixedlmerTest

2017 年 5 月更新：事实证明，我在这里写的一些内容有点错误。在整个帖子中进行了一些更新。

我非常同意 Ben Bolker 已经说过的话（感谢您对的大喊afex::mixed()），但让我在这个问题上添加一些更一般和具体的想法。

关注固定效应与随机效应以及如何报告结果

对于 Jonathan Baron 的示例数据集中表示的实验研究类型，您使用的重要问题通常是操纵因素是否具有整体影响。例如，我们是否发现的整体主效应或交互作用Task？重要的一点是，在这些数据集中，通常所有因素都处于完全的实验控制之下并且是随机分配的。因此，人们关注的焦点通常是固定效应。
相比之下，随机效应分量可以被视为捕捉系统方差（即效应大小的个体间差异）的“讨厌”参数，这些参数对于主要问题不一定重要。从这个角度来看，建议使用 Barr 等人提倡的最大随机效应结构。顺其自然。很容易想象，实验操作不会以完全相同的方式影响所有个体，我们希望对此进行控制。另一方面，因子或水平的数量通常不会太大，因此过度拟合的危险似乎相对较小。

因此，我会遵循 Barr 等人的建议。并指定一个最大随机效应结构并报告固定效应的测试作为我的主要结果。为了测试固定效应，我还建议使用afex::mixed()它报告效应或因素的测试（而不是参数的测试）并以某种合理的方式计算这些测试（例如，对所有模型使用相同的随机效应结构去除单一效应，使用总和到零的对比度，提供不同的方法来计算p值，...）。

示例数据呢

您提供的示例数据的问题在于，对于此数据集，最大随机效应结构会导致模型过饱和，因为设计的每个单元格只有一个数据点：

> with(df, table(Valence, Subject, Task))
, , Task = Cued

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

, , Task = Free

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

因此，lmer对最大随机效应结构的扼流圈：

> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable

不幸的是，据我所知，没有商定的方法来处理这个问题。但让我草拟并讨论一些：

第一个解决方案可能是移除最高随机斜率并测试此模型的效果：
```
require(afex)
mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
        Effect    F ndf  ddf F.scaling p.value
1         Task 6.56   1 4.00      1.00     .06
2      Valence 0.80   2 3.00      0.75     .53
3 Task:Valence 0.42   2 8.00      1.00     .67
```
~~然而，这个解决方案有点临时性，并没有过度激励。~~

2017 年 5 月更新：这是我目前认可的方法。请参阅此博客文章和我正在合着的章节的草稿， “传统 ANOVA 设计的随机效应结构”部分。
另一种解决方案（并且可以被视为 Barr 等人的讨论所倡导的解决方案）可能是始终移除随机斜率以获得最小的影响。但是这有两个问题：（1）我们使用哪种随机效应结构来找出最小的效应是什么，以及（2）如果更高阶的效应（例如这种效果的相互作用是存在的（见这里）。因此，需要手动设置这种随机效应结构，并将如此构建的模型矩阵传递给 lmer 调用。
第三种解决方案可能是使用随机效应部分的替代参数化，即对应于该数据的 RM-ANOVA 模型的参数化。不幸的是（？），lmer不允许“负方差”，所以这个参数化并不完全对应于所有数据集的 RM-ANOVA ，请参阅此处和其他地方的讨论（例如此处和此处）。这些数据的“lmer-ANOVA”将是：
```
> mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
        Effect    F ndf  ddf F.scaling p.value
1         Task 7.35   1 4.00      1.00     .05
2      Valence 1.46   2 8.00      1.00     .29
3 Task:Valence 0.29   2 8.00      1.00     .76
```

考虑到所有这些问题，我根本不会用它lmer来拟合设计的每个单元格只有一个数据点的数据集，除非有更一致的解决方案来解决最大随机效应结构的问题。

~~相反，我~~仍然可以使用经典的 ANOVA。car::Anova()在结果中使用其中一个包装器afex将是：

> aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
        Effect         df  MSE      F  ges   p
1      Valence 1.44, 5.75 4.67   1.46  .02 .29
2         Task       1, 4 4.08 7.35 +  .07 .05
3 Valence:Task 1.63, 6.52 2.96   0.29 .003 .71

请注意，afex现在还允许返回适合的模型，该模型aov可以传递给lsmeans事后测试（但对于效果测试，报告的效果car::Anova仍然更合理）：

> require(lsmeans)
> m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
> lsmeans(m, ~Task+Valence)
 Task Valence lsmean       SE   df lower.CL upper.CL
 Cued Neg       11.8 1.852026 5.52  7.17157 16.42843
 Free Neg       10.2 1.852026 5.52  5.57157 14.82843
 Cued Neu       13.0 1.852026 5.52  8.37157 17.62843
 Free Neu       11.2 1.852026 5.52  6.57157 15.82843
 Cued Pos       13.6 1.852026 5.52  8.97157 18.22843
 Free Pos       11.0 1.852026 5.52  6.37157 15.62843

Confidence level used: 0.95

其它你可能感兴趣的问题

上一篇检测时间序列的变化（R 示例）下一篇精确召回曲线的良好 AUC 是多少？