對於在 lmer+lmerTest 中具有交互作用的混合模型，summary() 和 anova() 的結果衝突
我最近遇到了一個問題，我認為包中的
anova()函數如何lmerTest計算其 F 統計量和來自混合效應模型的固定效應的相應 P 值。首先讓我說，我知道圍繞從混合效應模型計算 P 值的爭議（原因在這裡討論）。儘管如此，許多人仍然想要 P 值，因此已經開發了許多方法來適應這種情況（參見此處）。在這裡，我想展示一種常用方法的結果——即包中的anova函數lmerTest——並希望有人知道為什麼結果不太有意義。首先是我的數據。由於它的大小，我不得不鏈接到它。請注意，生物量列已標準化（平均值 = 0，sd = 1），因此為負值。這不會改變輸出。下載並指定工作目錄後，可以按如下方式讀取文件：
dat <- read.csv("StackOverflow_Data.csv", header = T)下面是我使用
lmer函數的模型lme4。在這個模型中，我將植物生物量作為響應變量和三個因素——A、B 和 C——每個都有兩個水平作為預測變量。植物基因型和空間塊作為隨機效應包括在內。model <- lmer(Biomass ~ A + B + C + A:B + A:C + B:C + A:B:C + (1 | Genotype) + (1 | Block) , data = dat, REML = T)總結上述模型，
summary(model)我們得到：Linear mixed model fit by maximum likelihood t-tests use Satterthwaite approximations to degrees of freedom [lmerMod] Formula: Biomass ~ A + B + C + A:B + A:C + B:C + A:B:C + (1 | Genotype) + (1 | Block) Data: dat AIC BIC logLik deviance df.resid 1059.7 1111.0 -518.8 1037.7 776 Scaled residuals: Min 1Q Median 3Q Max -3.04330 -0.63914 0.00315 0.69108 2.82368 Random effects: Groups Name Variance Std.Dev. Genotype (Intercept) 0.07509 0.2740 Block (Intercept) 0.01037 0.1018 Residual 0.19038 0.4363 Number of obs: 787, groups: Genotype, 50; Block, 6 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 2.27699 0.08162 47.50000 27.897 < 2e-16 *** AYes -0.02308 0.09958 99.30000 -0.232 0.81719 BReduced -0.11036 0.06232 733.00000 -1.771 0.07700 . CSupp -0.02152 0.06243 733.70000 -0.345 0.73039 AYes:BReduced 0.25113 0.08838 733.70000 2.841 0.00462 ** AYes:CSupp 0.02179 0.08854 734.50000 0.246 0.80567 BReduced:CSupp 0.19436 0.08838 733.10000 2.199 0.02817 * AYes:BReduced:CSupp -0.21746 0.12507 734.20000 -1.739 0.08251 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) AYes BRedcd CSupp AYs:BR AYs:CS BRd:CS AYes -0.607 BReduced -0.379 0.311 CSupp -0.379 0.311 0.498 AYes:BRedcd 0.269 -0.444 -0.706 -0.354 AYes:CSupp 0.268 -0.444 -0.352 -0.708 0.503 BRedcd:CSpp 0.267 -0.219 -0.706 -0.705 0.498 0.500 AYs:BRdc:CS -0.190 0.315 0.500 0.502 -0.709 -0.709 -0.708上面的摘要使用該
lmerTest包使用 Satterthwaites 對分母自由度的近似值從 t 統計量計算 P 值。由此我們看到，在 p = 0.05 的水平上，A:B和交互作用都顯著。B:C理論上，這些結果至少在質量上應該與包anova()中的函數產生的結果一致，該函數以lmerTest相同的方式計算 P 值。然而事實並非如此；這是來自 的輸出anova(model, type = 3)。注意 SS 的類型 3 測試Analysis of Variance Table of type III with Satterthwaite approximation for degrees of freedom Sum Sq Mean Sq NumDF DenDF F.value Pr(>F) A 0.09492 0.09492 1 49.87 0.4986 0.48342 B 0.66040 0.66040 1 732.66 3.4688 0.06294 . C 0.20207 0.20207 1 733.90 1.0614 0.30324 A:B 0.99470 0.99470 1 732.56 5.2247 0.02255 * A:C 0.36903 0.36903 1 733.66 1.9383 0.16427 B:C 0.35867 0.35867 1 733.20 1.8839 0.17031 A:B:C 0.57552 0.57552 1 734.23 3.0230 0.08251 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1這些結果明顯不同。
B:C交互作用不再顯著，交互作用的 P 值要高A:B得多。兩種模型都應該以相似的方式計算 P 值，因此很難想像它們會有如此不同。為什麼它們不同？
這是原始問題的一部分，但它可能會產生誤導，請參閱下面的答案。
實際上，看起來該
anova(model, type = 3)功能實際上是在使用 2 型 SS，我們可以通過運行來驗證anova(model, type = 2)。Analysis of Variance Table of type II with Satterthwaite approximation for degrees of freedom Sum Sq Mean Sq NumDF DenDF F.value Pr(>F) A 0.09526 0.09526 1 49.87 0.5004 0.48263 B 0.65996 0.65996 1 732.66 3.4665 0.06302 . C 0.19639 0.19639 1 733.91 1.0315 0.31013 A:B 0.99282 0.99282 1 732.56 5.2148 0.02268 * A:C 0.37018 0.37018 1 733.65 1.9444 0.16362 B:C 0.35523 0.35523 1 733.20 1.8659 0.17237 A:B:C 0.57552 0.57552 1 734.23 3.0230 0.08251 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1結果非常相似，考慮到模型中存在交互，情況不應如此。為了證明它
lmerTest::anova()實際上使用的是 2 型 SS，而不是它在輸出中顯示的 3 型 SS，我們可以使用包中的Anova()函數car。Anova(model, type = 2, test.statistic = 'F')產生：Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df) Response: Biomass F Df Df.res Pr(>F) A 0.4857 1 48.28 0.48917 B 3.4537 1 726.63 0.06351 . C 1.0337 1 727.77 0.30962 A:B 5.1456 1 726.54 0.02360 * A:C 1.9302 1 727.55 0.16517 B:C 1.8776 1 727.12 0.17103 A:B:C 2.9915 1 728.06 0.08413 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1請注意，對於我的數據，使用 Kenward-Roger ddf 並沒有太大改變結果。很清楚的是，2 型 SS 來自於
Car封裝的生產結果，類似於 3 型 SS 來自於lmerTest包裝的結果。這表明該lmerTest包實際上是計算類型 2 SS。我很難弄清楚為什麼會出現這種情況，除非從lmerTest包中計算 P 值存在問題。我錯過了什麼嗎？歡迎任何建議或想法。非常感謝！
編輯：2016 年 12 月 6 日，上午 11:40
一些人表示這個問題是從這裡重複的。但是我不明白這是怎麼回事。那篇文章旨在了解為什麼
aov()並lme()產生不同的 F 統計量，結果證明這與如何從不同的函數計算方差分量有關。在這裡，我只運行一個模型，使用並lmer試圖理解為什麼會產生不同的 P 值，儘管它們應該以類似的方式計算。似乎使用的是 2 型 SS 而不是報告的 3 型 SS，這僅在存在交互的情況下才重要，另一篇文章未包含在任何列出的模型中。lmerTest::anova(model)``summary(model)``lmerTest::anova()
毫無疑問，這是一個新手錯誤，我已經解決了我之前在上面詳述的問題。獲得
lmerTest::anova()和summary()同意的關鍵是確保通過輸入來改變 R 中的全局對比options(contrasts = c("contr.sum","contr.poly"))如果未設置這些對比，
summary()則將生成默認類型 I 平方和（請注意，car::Anova()也將使用此默認值）。另一方面lmerTest::anova()，無論默認對比度如何，都會產生 III 型 SS，這表明它在內部執行此操作。如果如上所述設置默認對比度，則所有三個函數都將使用所需的 III 型 SS 產生輸出。