ANOVA
黃哲璿
2022/07/17
什麼是ANOVA (Analysis of variance)
參考: Biostatistics for Dummies, P164~168)
The R book, P498~536)
上一堂課我們介紹了T test, 用於比較兩組之間的平均。
但是在實驗中我們常常會擁有更多組別(e.g.,安慰劑,A藥品,B藥品)。
分析這樣的資料時第一個直覺會想多次使用T test來比較。
例:安慰劑vs.A; A vs. B; 安慰劑vs.B 想想看這樣的比較會有什麼問題?
A: Inflated alpha (increase type 1 error)
因為上述多重比較的問題讓我們不能直接使用T test。
相對的,我們的第一個要問的問題是:所有組別的有一樣的平均嗎? (i.e., H0: 所有的組別平均相同; HA:其中有一或多個組別平均不同)
ANOVA並不會成對的比較兩組之間的差異,而是回答上述的問題並產生相對應的P值。
要回答這個問題,ANOVA計算了組間variance和組內variance的比率:F ratio。 (所以被稱作Analysis of variance) 直觀的來想,組間的variance >> 組內的variace時,可以想見並非所有的組別的平均都是相同的。
相反的,當虛無假說為真時, F ratio會接近於1,而它的分布會接近Fisher F distribution。 因此,我們來利用F分布跟F ratio來計算P值。
在介紹計算過程之前，我們先來看看在R之中如何使用ANOVA。
set.seed(123)
# 創造四組數字，其中只有b的平均不同
a <- rnorm(10,10)
b <- rnorm(10,20)
c <- rnorm(10,10)
d <- rnorm(10,10)
# data 1 中包含了平均不同的B
data1 <- tibble(a = a, b = b, c = c) %>% gather(key = "group")
# aov是執行anova最簡單的function
aov1 <- aov(value ~ group, data1)
summary(aov1)## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 720.0 360 378.5 <2e-16 ***
## Residuals 27 25.7 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# anova是線性回歸的一種特例,因此我們可以用lm() + anova()得到相同結果
lm1 <- lm(value ~ group, data1)
anova(lm1)## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 720.04 360.02 378.48 < 2.2e-16 ***
## Residuals 27 25.68 0.95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# data 2 中全部平均相同
data2 <- tibble(a = a, d = d, c = c) %>% gather(key = "group")
aov2 <- aov(value ~ group, data2)
summary(aov2)## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 2.893 1.4464 2.112 0.141
## Residuals 27 18.487 0.6847只有數字跟結果可能很難理解，讓我來視覺化看看:
# 點到黑線(組間平均)的距離明顯小於點到藍線(整體平均)的距離
ggplot(data1, aes(group, value)) +
geom_jitter(height = 0) +
stat_summary(
fun.y = mean,
geom = "errorbar",
aes(ymax = ..y.., ymin = ..y..),
position = position_dodge(width = 0.3),
width = 1) +
geom_hline(yintercept = mean(data1$value), color = "blue")
# 點到黑線(組間平均)的距離大約等於點到藍線(整體平均)的距離
ggplot(data2, aes(group, value)) +
geom_jitter(height = 0) +
stat_summary(
fun.y = mean,
geom = "errorbar",
aes(ymax = ..y.., ymin = ..y..),
position = position_dodge(width = 0.3),
width = 1) +
geom_hline(yintercept = mean(data2$value), color = "blue")
如何計算ANOVA與理解variance table
Variance table:
anova(lm1)## Analysis of Variance Table
##
## Response: value
## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 720.04 360.02 378.48 < 2.2e-16 ***
## Residuals 27 25.68 0.95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sum Sq = Sum of squares = \(\sum(y - \overline{y})^2\)
Total Sum of squares (SST) = \(\sum\) (y - overall mean)^2 = \(\sum(y - \overline{\overline{y}})^2\)
(SST <- data1 %>% mutate(overall_mean = mean(value)) %>%
summarise(SST = sum((value - overall_mean)^2)) %>% .[[1]])## [1] 745.722
Error Sum of squares (SSE): 組內的Sum of squares
(SSE <- data1 %>% group_by(group) %>% mutate(group_mean = mean(value)) %>%
summarise(SSE = sum((value - group_mean)^2)) %>% ungroup() %>%
summarise(SSE = sum(SSE)) %>% .[[1]])## [1] 25.68336
Sum of squares absorbed by group (SSA) = SST - SSE: 組間的Sum of squares
(SSA <- SST-SSE)## [1] 720.0387由於一組內的N如果增加,sum of squares也會增加,為了去除N的影響我們要計算sum of squares 的平均: mean square (Mean Sq), 事實上 mean square就是variation
計算MS: SS/df
組間的df = number of group -1 (i.e., 3 - 1 = 2)
組內的df = total number of observations minus the number of groups (i.e., 30 - 3 = 27)
因此:
df_n <- 3 - 1
# 組間的variation
(MSA <- SSA/df_n) ## [1] 360.0193# 組內的variation
df_d <- 30 - 3
(MSE <- SSE/df_d) ## [1] 0.9512354知道組內和組間的variation後,我們就可以計算F ratio
(F_ratio <- MSA/MSE) ## [1] 378.4755知道F ratio 後,我們可以用F 分布計算P值
F 分布的形狀會由兩個數值決定，剛好是我們剛剛計算的
numerator degrees of freedom: number of groups - 1
denominator degrees of freedom: number of observations - number of groups
pf(F_ratio, df_n, df_d, lower.tail = F)## [1] 1.780368e-20anova(lm1)$`Pr(>F)`## [1] 1.780368e-20 NAggplot() +
geom_function(fun = ~ df(x = .x, df_n, df_d)) +
xlim(c(0, 400)) +
geom_vline(aes(xintercept = F_ratio)) +
labs(x = "F-ratio", y = "density")
順道看看各種F分布的形狀吧~
# 自由度對F分布的影響
F_curve <- function(df1, df2){
ggplot() +
geom_function(fun = df, args = list(df1 = df1, df2 = df2))
}
F_curve <- function(df1, df2){
ggplot() +
geom_function(fun = "df", args = list(df1 = df1, df2 = df2)) +
scale_x_continuous(limits = c(0, 5), breaks = 0:5) +
geom_vline(aes(xintercept = qf(p = 0.95, df1 = df1, df2 = df2))) +
labs(x = "F-ratio", y = "density",
subtitle = str_c("df1 = ", df1, ", df2 = ", df2))
}
d <- expand.grid(1:3,c(10,100,1000))
map2(d[[1]],d[[2]], F_curve) %>%
patchwork::wrap_plots()
Post hoc test
但回到一開始的問題，我們要怎麼知道A組跟安慰劑之間有沒有差異呢?
這時候我們要進行post-hoc test(事後檢驗)。
顧名思義，這是在ANOVA檢驗後確定有一組的平均明顯不同，用來檢驗是如何的不同。
Post-hoc的出發點是為了避免type-1 error 的增加，而根據你比較的方式，使用的post hoc test也不同。
以下介紹幾種常見的post hoc test 並運用在我們的例子之中。
The Bonferroni test
最保守但也是最泛用的test，它進行了多次的t test來比較平均，再事後調整p_value (乘上總共比較了幾次)
means <- emmeans(lm1, ~ group)
contrast(means,method = "pairwise", adjust = "bonferroni")## contrast estimate SE df t.ratio p.value
## a - b -10.134 0.436 27 -23.234 <.0001
## a - c 0.499 0.436 27 1.144 0.7874
## b - c 10.633 0.436 27 24.378 <.0001
##
## P value adjustment: bonferroni method for 3 testst.test(a, c, var.equal = T)##
## Two Sample t-test
##
## data: a and c
## t = 1.1845, df = 18, p-value = 0.2516
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3862263 1.3845954
## sample estimates:
## mean of x mean of y
## 10.074626 9.575441The Tukey’s test
利用Student化全距分布(Studentized range distribution)計算正確的P值。
是進行pairwise比較(A-B, A-C, B-C)時，最正確的post-hoc test，但不應用於其他比較方式(見Dunnette’s test & Scheffe’s test)。
contrast(means, method = "pairwise", adjust = "tukey")## contrast estimate SE df t.ratio p.value
## a - b -10.134 0.436 27 -23.234 <.0001
## a - c 0.499 0.436 27 1.144 0.4958
## b - c 10.633 0.436 27 24.378 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimatesThe Fisher’s LSD (least significant difference) test:
產生的P值不正確，請勿使用。
The Scheffe’s test:
它允許你將兩個組別綁再一起跟另外一個組別比較。
例如: 藥B和藥C 是否優於安慰劑。
是相當保守(不易有type I error)的test。
A = c(1, 0, 0)
B = c(0, 1, 0)
C = c(0, 0, 1)
contrast(means, method = list("A - (B+C)/2" = A - (B + C)/2),
adjust = "scheffe")## contrast estimate SE df t.ratio p.value
## A - (B+C)/2 -4.82 0.378 27 -12.753 <.0001
##
## P value adjustment: scheffe method with rank 1Note: 參考: https://aosmith.rbind.io/2019/04/15/custom-contrasts-emmeans/#treatment-vs-control-comparisons
The Dunnett’s test:
用於一對多的比較。例如:安慰劑vs.藥B; 安慰劑vs.藥C。
由於比較的次數變少(不比較B跟C)，它的檢定力會比pairwise高。
contrast(means, method = "trt.vs.ctrl", adjust = "dunnett", ref = "a")## contrast estimate SE df t.ratio p.value
## b - a 10.134 0.436 27 23.234 <.0001
## c - a -0.499 0.436 27 -1.144 0.4286
##
## P value adjustment: dunnettx method for 2 tests更多關於多重比較在R的操作方法可以參照(需先有迴歸分析的基礎)
https://cran.r-project.org/web/packages/emmeans/emmeans.pdf
https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html
由於多重比較是每個研究中一定會用到的分析技巧，個人強力推薦好好弄懂這個觀念。