混合方差分析正态性：应该检查哪些变量？（在 stats::aov 的普遍和实际应用中）

机器算法验证 r 方差分析 重复测量 残差 正态假设
2022-03-30 21:28:44

我花了很多时间阅读书籍章节、文章、在线教程等，但没有明确的答案（主要是因为它们只描述了单向方差分析或其他非常具体的应用程序）。该站点上也有许多类似的问题，但对于我的目的而言，再次没有令人满意的答案。

本质上，我想知道在给定任意数量的主题内或主体间因素（具有任意数量的水平）。

（注意：这里唯一的问题是应该检查哪些变量，而不是如何检查它们。通过“检验/检验正态性”，我不一定指统计假设检验，它也可以基于密度或 QQ 图，等等，没关系。唯一的问题是是否可能需要多元正态性检验，在这种情况下，问题又是应该包含哪些变量。）

至少本教程和这个答案建议检查每个单元格的正态性，即每个因素的每个级别的每个可能组合——但没有给出参考或详细的推理，这对于复杂的设计似乎相当极端。但大多数其他人（例如this或this或this answer或this book chapter或this video tutorial）建议只检查残差（无论因子内/因子之间如何）。即使我假设这是后者，问题仍然存在：应该检查哪些残差？

在下文中，我使用R函数stats:aov输出来举例说明一些可能的答案。

我准备了一个发明的数据集用于说明。每个单独的主题都用“ subject_id”表示。有两个主体间因素：“ btwn_X”和“ btwn_Y”。还有两个主体内因素：“ wthn_X”和“ wthn_Y”。

# preparing some invented data    
dat_example = data.frame(
    subject = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
    btwn_X = c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2),
    btwn_Y = c(1, 2, 1, 2, 2, 1, 1, 1, 2, 1),
    measure_x1_yA = c(36.2, 45.2, 41, 24.6, 30.5, 28.2, 40.9, 45.1, 31, 16.9),
    measure_x2_yA = c(-14.1, 58.5, -25.5, 42.2, -13, 4.4, 55.5, -28.5, 25.6, -37.1),
    measure_x1_yB = c(83, 71, 111, 70, 92, 75, 110, 111, 110, 85),
    measure_x2_yB = c(8.024, -14.162, 3.1, -2.1, -1.5, 0.91, 11.53, 18.37, 0.3, -0.59),
    measure_x1_yC = c(27.4,-17.6,-32.7, 0.4, 37.2, 1.7, 18.2, 8.9, 1.9, 0.4),
    measure_x2_yC = c(7.7, -0.8, 2.2, 14.1, 22.1, -47.7, -4.8, 8.6, 6.2, 18.2)
)
dat_example$subject = as.factor(as.character(dat_example$subject))
dat_example$btwn_X = as.factor(as.character(dat_example$btwn_X))
dat_example$btwn_Y = as.factor(as.character(dat_example$btwn_Y))    
vars = c(
    'measure_x1_yA',
    'measure_x2_yA',
    'measure_x1_yB',
    'measure_x2_yB',
    'measure_x1_yC',
    'measure_x2_yC'
)
dat_l = stats::reshape(
    dat_example,
    direction = 'long',
    varying = vars,
    idvar = 'subject',
    timevar = "within_factor",
    v.names = "values",
    times = vars
)    
dat_l$wthn_X = sapply(strsplit(dat_l$within_factor, split = '_', fixed =
                                   TRUE), `[`, 2)
dat_l$wthn_Y = sapply(strsplit(dat_l$within_factor, split = '_', fixed =
                                   TRUE), `[`, 3)
dat_l$wthn_X = as.factor(as.character(dat_l$wthn_X))
dat_l$wthn_Y = as.factor(as.character(dat_l$wthn_Y))

# performing the ANOVA    
aov_BBWW = aov(values ~ btwn_X * btwn_Y * wthn_X * wthn_Y +
                   Error(subject / (wthn_X * wthn_Y)), data = dat_l)

（另请参阅此处lme4::lmer的扩展版本，其中包含各种因子变化和模型之间/内部。）

aov 对象aov_BBWW返回以下内容：

Grand Mean: 23.6847

Stratum 1: subject

Terms:
                  btwn_X   btwn_Y btwn_X:btwn_Y Residuals
Sum of Squares    61.549  351.672        18.969  3221.628
Deg. of Freedom        1        1             1         6

Residual standard error: 23.17192
15 out of 18 effects not estimable
Estimated effects may be unbalanced

Stratum 2: subject:wthn_X

Terms:
                   wthn_X btwn_X:wthn_X btwn_Y:wthn_X btwn_X:btwn_Y:wthn_X Residuals
Sum of Squares  23432.120       612.948       712.387              773.779   513.165
Deg. of Freedom         1             1             1                    1         6

Residual standard error: 9.248106
8 out of 12 effects not estimable
Estimated effects may be unbalanced

Stratum 3: subject:wthn_Y

Terms:
                   wthn_Y btwn_X:wthn_Y btwn_Y:wthn_Y btwn_X:btwn_Y:wthn_Y Residuals
Sum of Squares  19262.400       982.159      1561.578             1836.188  5860.787
Deg. of Freedom         2             2             2                    2        12

Residual standard error: 22.09975
8 out of 16 effects not estimable
Estimated effects may be unbalanced

Stratum 4: subject:wthn_X:wthn_Y

Terms:
                wthn_X:wthn_Y btwn_X:wthn_X:wthn_Y btwn_Y:wthn_X:wthn_Y
Sum of Squares      20248.558              159.421              986.331
Deg. of Freedom             2                    2                    2
                btwn_X:btwn_Y:wthn_X:wthn_Y Residuals
Sum of Squares                      604.163  4789.399
Deg. of Freedom                           2        12

Residual standard error: 19.9779
Estimated effects may be unbalanced

我可以访问以下残差（有关更多详细信息，请参见此处）：

aov_BBWW$subject$residuals
aov_BBWW$`subject:wthn_X`$residuals
aov_BBWW$`subject:wthn_Y`$residuals
aov_BBWW$`subject:wthn_X:wthn_Y`$residuals
aov_BBWW$`(Intercept)`$residuals

根据上面引用的一些资料，这些残差应用于正态性检验，但尚不清楚是全部还是仅一个（以及在这种情况下是哪一个）。

编辑：

经过大量挖掘（并在 EdM 的回答和评论的帮助下），最权威的解决方案似乎是，在只有主体间因素的 ANOVA 的情况下，正确的变量只是residuals来自 aov 对象的向量（例如aov_BB$residuals） ，而如果有任何主题内变量，我应该这样做：

aov_proj = proj(aov_BBWW)
aov_proj[[length(aov_proj)]][,"Residuals"]

后者是要检查正态性和其他相关假设的变量。为什么会这样，我无法理解，但几个看似自信的消息来源给出了这个解决方案：this and this R mailing list replys, this and this and this CV answers（讽刺的是后两个不是被接受的），本教程和MASS文档. 大多数或所有这些来源都来自Venables 和 Ripley (2002)，但我认为他们不会盲目地复制不正确的东西。

尽管如此，问题仍然悬而未决：我很乐意收到关于此事的进一步验证（或反驳）和解释。

（顺便说一句，如果要信任上述来源，则拟合值显然可以通过以下方式访问fitted(aov_BBWW[[length(aov_BBWW)]])：）

1个回答

TL;DR：ANOVA 在所有观察中汇集信息，以获得对固定效应、随机效应和误差方差的最佳估计。如果您想检查 ANOVA 残差的正态性，那么在考虑所有固定和随机效应之后这样做是最有意义的。可靠的 ANOVA 估计不需要残差的正态性；问题是测试统计的分布。在重复测量方差分析中，相关结构的不平衡或错误指定等问题可能是可靠统计测试的更大障碍。

ANOVA 只是一种特定类型的线性模型，例如在从问题链接的其中一个站点的此页面上进行了描述，并在此处进行了广泛讨论。与所有线性模型一样，ANOVA 将来自预测值组合的信息结合起来，将结果值建模为预测值加上误差项的函数。假设误差项在所有情况下都有一定的分布，标准方差分析为零均值的高斯分布。有关误差项分布的信息是通过汇集所有观察结果来获得的，消除在 ANOVA 设计的单个单元格中偶然发生的变幻莫测。因此，标准正态 qq 诊断图检查所有残差值，而不是单个单元格内的残差值。

尽管 ANOVA 模型中通常假设高斯误差，但显着性检验并不一定要求满足该假设。ANOVA 中的显着性检验是对回归系数的检验。因此，当执行标准参数测试时，这些回归系数的抽样分布必须充分满足假设。

正如@whuber 在一个至关重要的评论中所说：

您真正想知道的是 ANOVA 检验统计量的假设分布是否足够准确以计算您感兴趣的 p 值。

如果满足模型假设并且共享误差项具有高斯分布，那么您就知道回归系数的检验将是有效的。* 但误差项的严格正态性对于回归系数的检验是有效的不是必需的。对于线性模型回归系数（包括方差分析）进行充分可靠的显着性检验，将正态分布的误差项视为足够但并不总是必要的。

这并不是说检查包含所有案例信息的模型预测周围的残差分布是没有用的。例如，Rlme4包提供了一个正常的 qq 图作为其诊断图之一；请参阅小插图的第 33 页。但是，您经常会发现，在这样的残差图中，与正态性的显着偏差意味着模型本身的指定很差。这可能是这样一个情节中最有用的信息。

对于仅具有固定分类预测变量并包括所有交互作用的混合 ANOVA 模型，您不必担心固定效应预测变量本身的线性。但是可能存在对结果变量的错误处理（例如，如果它基本上是对数正态而不是正态），遗漏与结果和包含的预测变量相关的关键协变量，或者随机效应结构的错误指定。修复诊断图暴露的那些问题，而不是沉迷于常态本身。

为了评估模型，应检查所有诊断图：不仅是残差正态性的 qq 图，还有拟合残差图、比例位置图和各种剖面图（参见小插图的第 36 页）混合模型及其随机效应。检查特定观察的过度影响，例如使用 R 中的影响.ME 包。这个过程，而不是简单的正态性检查，对于评估和提高模型规范的质量至关重要。

如果模型被正确指定，那么回归系数抽样分布的正态性假设可以是相当可靠的。有了足够的数据，尽管有非正态残差，中心极限定理可以帮助解决这个问题，尽管有多少数据“足够”取决于具体情况。例如，请参阅此答案。如果您不想依赖该假设，则自举提供了一种获取非参数置信区间的方法。但是只有在模型本身被充分指定时才应该这样做。

作为对问题注释的编辑，可以从由 分析的重复测量数据生成一些诊断图aov，根据其手册页，这些数据适合“通过调用lm每个层来分析方差模型”。每个层都是逐步复杂模型的观测均值的一部分，从总体均值开始。正如Venables 和 Ripley在第 283 页关于更简单的裂区设计所说：

多层模型可以使用 来拟合aov，并由以下形式的模型公式指定

响应〜mean.formula +错误（strata.formula）

在我们的示例中，strata.formula是 B/V，指定地层 2 和 3；第四层被自动包含为“内部”层，即层公式中的剩余层。

对于更复杂的模型，最后一层因此是自动包含的“内部”层。继续第 284 页：“不可能将 [最后一层的拟合值和残差] 与原始实验的图唯一地关联起来。” 您需要“将原始数据向量投影到方差分析表中每行定义的子空间”中的残差。可以检查每个层的残差，但只有最终层会考虑模型的所有方面。这个答案显示了 Venables 和 Ripley 示例的代码，其中第四层是“内部”层。

但是，在继续之前aov，请注意其帮助页面中的以下引用：

笔记

aov是为平衡设计而设计的，如果没有平衡，结果可能很难解释：请注意，响应中的缺失值可能会失去平衡。如果有两个或多个错误层，则使用的方法在没有平衡的情况下在统计上是低效的，使用lmepackage可能会更好nlme。

*这对于混合模型更为复杂，对于测试中使用的自由度数量存在争议。但是这个争议不会通过检查残差的分布来解决。混合模型的测试还可能涉及有关相关观察的协方差结构的假设。

其它你可能感兴趣的问题
上一篇对于变分自动编码器，重建损失应该计算为输入的总和还是平均值？ 下一篇在 R 中解释负二项式混合模型输出