r - 嵌套方差分析唯一因子水平

Question

我正在运行具有以下设置的嵌套方差分析:2 个区域，一个是引用，一个是曝光(名为 CI = Control/Impact 的列)。两个时间段(影响之前和之后，名为 BA 的列)，之前时间段为 1 年，之后时间段为 3 年。岁月是嵌套的。

我的问题是:如果我使用原始年份(在玩具数据集中的 Time2 列中)，我会得到一个结果。如果我重命名年份，以便它们在 Before 和 After 中分别为 1 和 1-3，我会得到不同的结果。

问题:

由于年份是唯一的，嵌套不应该考虑正确的数据结构吗？

对于使用类型 1 或类型 2 SS 的两个模型，我的结果是相同的。为什么它们在不同类型的 SS 之间没有区别？我本来希望 BA(但不是 CI)会改变。

使用唯一名称的类型 3 SS 的结果与类型 1 和 2 SS 的结果相同。但是使用重命名的值，我得到了不同的结果。到底是怎么回事？

玩具数据集:

toy <- structure(list(BA = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
), .Label = c("A", "B"), class = "factor"), Time = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("1", "2", "3"), class = "factor"), 
Time2 = structure(c(4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L
), .Label = c("11", "12", "13", "15", "16", "17"), class = "factor"), 
Lake = c("Area 1", "Area 1", "Area 1", "Area 1", "Area 1", 
"Area 2", "Area 2", "Area 2", "Area 2", "Area 2", "Area 1", 
"Area 1", "Area 1", "Area 1", "Area 1", "Area 2", "Area 2", 
"Area 2", "Area 2", "Area 2", "Area 1", "Area 1", "Area 1", 
"Area 1", "Area 1", "Area 2", "Area 2", "Area 2", "Area 2", 
"Area 2", "Area 1", "Area 1", "Area 1", "Area 1", "Area 1", 
"Area 2", "Area 2", "Area 2", "Area 2", "Area 2"), CI = structure(c(2L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L), .Label = c("C", "I"), class = "factor"), 
Response = c(78.3, 75.3, 69.4, 75.1, 71.1, 49.7, 61, 59.6, 
35.3, 26.5, 80.9, 81.4, 67.6, 73.6, 73, 46.4, 73.6, 67.1, 
34, 45.5, 86.6, 78, 68.2, 76.8, 69.6, 52.1, 61.9, 50.8, 39.2, 
49.6, 72, 74, 71, 68, 58, 40, 41, 34, 54, 61)), .Names = c("BA", 
"Time", "Time2", "Lake", "CI", "Response"), row.names = c(NA, 
40L), class = "data.frame")

使用类型 1 SS 进行分析:

mod <- lm(Response ~ BA + CI + BA*CI + BA/Time + BA/Time*CI, data = toy)
mod1 <- lm(Response ~ BA + CI + BA*CI + BA/Time2 + BA/Time2*CI, data = toy)
# results are the same
anova(mod)
anova(mod1)

现在尝试使用类型 2

library(car)
options(contrasts=c("contr.sum", "contr.poly"))
mod <- lm(Response ~ BA + CI + BA*CI + BA/Time + BA/Time*CI, data = toy)
mod1 <- lm(Response ~ BA + CI + BA*CI + BA/Time2 + BA/Time2*CI, data = toy)
Anova(mod, type = "II", singular.ok = TRUE)
Anova(mod1, type = "II", singular.ok = TRUE)

和类型 3

Anova(mod, type = "III", singular.ok = TRUE)
Anova(mod1, type = "III", singular.ok = TRUE)

Answer 1

一开始就警告，这其中的一些内容可能并不完全正确——我承认此时我对 ANOVA 已经生疏了。

在 III 类 SS 分析中，我们说主效应由相互作用限定。简而言之，这意味着在存在高阶相互作用的情况下，低阶相互作用和主效应的可解释性较差。传统将我们推向 III 型分析，但由于这个原因，它们确实很痛苦。反正...

让我们快速浏览一下您的重新编码。您重新编码的效果是 Time2在 Time 中从具有四个不同的值变为具有三个不同的值.您可能会保留一些可解释性，因为 BA 的组合在级别 B 的时间值是唯一的，以前是 13，但现在是 1。

让我们回到你的数据。现在，BA:Time2 携带与 BA:CI 相同的所有信息。这在您的 Time2 方差分析结果中看起来如何...

> Anova(mod1, type = "III", singular.ok = TRUE)
Anova Table (Type III tests)

Response: Response
            Sum Sq Df  F value    Pr(>F)    
(Intercept)  73945  1 712.5963 < 2.2e-16 ***
BA             246  1   2.3678    0.1337    
CI            2484  1  23.9401 2.713e-05 ***
BA:CI            0  1   0.0046    0.9462    
BA:Time2        95  2   0.4570    0.6372    
BA:CI:Time2     37  2   0.1797    0.8364    
Residuals     3321 32                 

... BA:CI 的 SS 与预期的一样(或多或少)为 0。将此与时间模型进行对比...

> Anova(mod, type = "III", singular.ok = TRUE)
Anova Table (Type III tests)

Response: Response
            Sum Sq Df   F value    Pr(>F)    
(Intercept) 107772  1 1038.5835 < 2.2e-16 ***
BA             209  1    2.0099    0.1659    
CI            4220  1   40.6655 3.661e-07 ***
BA:CI            9  1    0.0907    0.7653    
BA:Time         95  2    0.4570    0.6372    
BA:CI:Time      37  2    0.1797    0.8364    
Residuals     3321 32                        

BA:CI 必须保持一点差异……模型的其余部分似乎也做得更好。

我的感觉是，您在所有编码方案下的 ANOVA 数据集都指定得很差。在两种编码方式下，AB 的 B 级与您的时间组的一个级别混淆。还要特别注意 package:car for Anova 中的注释，尤其是关于 single.ok 的争论。简而言之，它说“对于 II 类测试默认为 TRUE，对于 III 类测试默认为 FALSE(其中具有混叠系数的模型的测试将无法直接解释)”......在这里你看到了看似混叠的系数。

……别介意福克斯对他在 II 型和 III 型测试中所做的描述是……难以解释。这提醒了我为什么我总是使用 package:ez for ezANOVA() 回到我不得不容忍非 I 型测试的那一天。