지난번 포스팅에서 일원분산분석(one-way ANOVA) 이론과 R 의 aov() 함수에 대해서 알아보았습니다. 그러면 이번 포스팅에서는 일원분산분석 후에 이어서 하게 되는 다중비교(Multiple Comparison)에 대해서 알아보겠습니다.
F-검정 통계량을 사용해서 한개의 요인(Factor) 내 세 개 이상의 수준(Levels)의 집단 간에 평균의 차이가 있는지 없는지를 검정하였는데, 만약 분석 결과가 "인자 수준들의 평균 μi 는 모두 같다"는 귀무가설 H0를 기각하고 "인자 수준들의 평균 μi 는 모두 같지 않다" 대립가설 H1을 채택할 경우가 생겼다고 합시다. 이럴 경우 수준(Levels)의 집단 간에 어디에서 차이가 생겼는지 알고 싶을 때 사용하는 분석법이 다중비교(multiple comparison)이 되겠습니다.
수준의 수가 2일 때 일원분산분석에서 귀무가설 H0가 기각될 경우에는 두 수준 평균간의 차이에 대해서는 t-Test를 적용하면 됩니다. 하지만 수준의 수(number of Levles, groups)가 3개 이상인 경우에는 t-Test를 사용하면 안됩니다. 왜냐하면 비교해야할 수준수가 커질수록 두 수준간의 짝들 간의 비교(pair-wise comaprison)를 하게 되면 제 1종 오류(참을 거짓으로 판정할 오류)가 애초 정한 유의수준 α(significance level α)가 아니라 이보다 훨씬 커지게 되기 때문입니다.
t-Test를 사용할 경우 무슨 일이 벌어지는지 예를 들어서 설명해보겠습니다. 만약 수준의 수(number of Levels, groups)가 5개일 경우 각 각 쌍을 지어서 비교를 한다고 했을 때 수준평균짝들의 총 갯수는 10개(5 Combination 2 = 4 + 3 + 2 + 1 = 10) (R의 combination 함수를 사용하면, choose(5,2) 하면 10이 나옴) 이므로 t-Test를 10번 적용해야 합니다. 만약 유의수준 α (significance level α) = 0.05라고 한다면, 10번의 t-Test를 각각 독립적으로 하고 나면 제1종 오류가 생길 확률은 1 - (0.95)^10 = 0.4012631 가 됩니다. 이는 처음의 유의수준 α 0.05보다 약 8배 크며, 1종 오류가 생길 확률도 약 8배 정도 크다는 소리입니다.
이러한 오류를 피하기 위해서는 유의수준을 설정된 크기로 유지하면서 모든 두 수준짝들의 평균을 동시에 비교할 수 있도록 고안된 검정방법인 다중비교(multiple comparison)으로 Tukey's HSD(honestly significant difference) test와 Duncan's LSR(least significant range) test 가 있는데요, 이번 포스팅에서는 Tukey's HSD test 에 대해서 소개하겠습니다.
Tukey's HSD(honestly significant difference) test는 studentized range distribution을 이용하여 모든 가능한 두 수준들의 평균간의 차이가 있는지를 검정(pairwise post-hoc testing using Tukey HSD test)하는 방법입니다. Tukey's HSD test 의 통계량 및 검정 방법은 아래와 같습니다.
[ Tukey's HSD(honestly significant difference) test ]
그러면, (1) R의 aov() 함수를 사용해서 유의수준 0.05 하에 일원분산분석(one-way ANOVA)를 먼저 해보고 (<= 이것은 지난번 포스팅과 동일합니다), (2) Tukey's HSD test를 통해서 다중비교(multiple comparison)을 하여 쌍을 이룬 집단간의 평균 비교 중에서 어디서 차이가 발생한건지 알아보도록 하겠습니다.
> ##----------------------------------------------------------
> ## One-way ANOVA : aov(), oneway.test
> ##----------------------------------------------------------
>
> ##--- Are there any daily outcome differences among temperature conditions?
> # group 1 : temperature condition 1
> # group 2 : temperature condition 2
> # group 3 : temperature condition 3
>
> # daily outcome by tmep condition (group 1/2/3)
> y1 <- c(50.5, 52.1, 51.9, 52.4, 50.6, 51.4, 51.2, 52.2, 51.5, 50.8)
> y2 <- c(47.5, 47.7, 46.6, 47.1, 47.2, 47.8, 45.2, 47.4, 45.0, 47.9)
> y3 <- c(46.0, 47.1, 45.6, 47.1, 47.2, 46.4, 45.9, 47.1, 44.9, 46.2)
>
> y <- c(y1, y2, y3)
> y
[1] 50.5 52.1 51.9 52.4 50.6 51.4 51.2 52.2 51.5 50.8 47.5 47.7 46.6 47.1 47.2 47.8 45.2 47.4 45.0
[20] 47.9 46.0 47.1 45.6 47.1 47.2 46.4 45.9 47.1 44.9 46.2
>
> n <- rep(10, 3)
> n
[1] 10 10 10
>
> group <- rep(1:3, n)
> group
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3
>
> # combining into data.frame
> group_df <- data.frame(y, group)
> group_df
y group
1 50.5 1
2 52.1 1
3 51.9 1
4 52.4 1
5 50.6 1
6 51.4 1
7 51.2 1
8 52.2 1
9 51.5 1
10 50.8 1
11 47.5 2
12 47.7 2
13 46.6 2
14 47.1 2
15 47.2 2
16 47.8 2
17 45.2 2
18 47.4 2
19 45.0 2
20 47.9 2
21 46.0 3
22 47.1 3
23 45.6 3
24 47.1 3
25 47.2 3
26 46.4 3
27 45.9 3
28 47.1 3
29 44.9 3
30 46.2 3
>
> sapply(group_df, class)
y group
"numeric" "integer"
>
> # transform from 'integer' to 'factor'
> group_df <- transform(group_df, group = factor(group))
> sapply(group_df, class)
y group
"numeric" "factor"
>
>
> # boxplot
> attach(group_df)
The following objects are masked _by_ .GlobalEnv:
group, y
> boxplot(y ~ group,
+ main = "Boxplot of Daily Outcome by Temperature condition 1/2/3",
+ xlab = "Factor Levels : Temperature condition 1/2/3",
+ ylab = "Daily Outcome")
>
>
> # descriptive statistics by group
> tapply(y, group, summary)
$`1`
Min. 1st Qu. Median Mean 3rd Qu. Max.
50.50 50.90 51.45 51.46 52.05 52.40
$`2`
Min. 1st Qu. Median Mean 3rd Qu. Max.
45.00 46.72 47.30 46.94 47.65 47.90
$`3`
Min. 1st Qu. Median Mean 3rd Qu. Max.
44.90 45.92 46.30 46.35 47.10 47.20
> detach(group_df)
>
> # one-wayANOVA
> aov(y ~ group, data = group_df)
Call:
aov(formula = y ~ group, data = group_df)
Terms:
group Residuals
Sum of Squares 156.302 19.393
Deg. of Freedom 2 27
Residual standard error: 0.8475018
Estimated effects may be unbalanced
> summary(aov(y ~ group, data = group_df))
Df Sum Sq Mean Sq F value Pr(>F)
group 2 156.30 78.15 108.8 1.2e-13 ***
Residuals 27 19.39 0.72
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1>
> # Bartlett test to test the null hypothesis of equal group variances
> bartlett.test(y ~ group, data = group_df)
Bartlett test of homogeneity of variances
data: y by group
Bartlett's K-squared = 1.6565, df = 2, p-value = 0.4368
> > ##------------------------------------------------------------------
> ## multiple comparison
> ## - Tukey's HSD(honestly significant difference) test : TukeyHSD()
> ##------------------------------------------------------------------ > group_aov <- aov(y ~ group, data = group_df)
> TukeyHSD(group_aov)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = y ~ group, data = group_df)
$group
diff lwr upr p adj
2-1 -4.52 -5.459735 -3.5802652 0.0000000
3-1 -5.11 -6.049735 -4.1702652 0.0000000
3-2 -0.59 -1.529735 0.3497348 0.2813795
|
위의 일원분산분석(one-way ANOVA) 결과를 보면 P-value가 1.2e-13 으로서 대립가설 H1을 채택하게 되어 평균이 서로 다르다고 판단할 수 있으므로, 자연스럽게 그 다음 질문은 "다르다고 하는데 그러면 쌍을 이룬 평균 비교 중에서 어디에서 차이가 난다는 말인가?"가 될 것입니다.
이 질문에 답하기 위해 위의 Tukey's HSD test 결과를 보면, "multiple comparisons of means, 95% family-wise confidence level'이라는 설명이 나오는데요, 2개씩 쌍을 이룬 수준간 평균의 다중 비교를 95% 신뢰수준으로 상한(upr)과 하한(lwr)의 신뢰계수 구간을 구했고, P-value도 구해주었습니다.
위의 결과를 보면 'group 2'와 'group 1'은 평균 차이가 -4.52 이고 adj. P-value가 0.0000000이고, 'group 3'과 'group 1'의 평균 차이는 -5.11이고 adj. P-value가 0.0000000 으로서, 유의수준 0.05보다 훨씬 작으므로 이들 group (수준, factor levels) 간에는 평균의 차이가 있다고 유의수준 0.05 기준 하에 판단할 수 있겠습니다. 반면에, 'group 3'과 'group 2'는 평균 차이가 -0.59이고 adjusted P-value 가 0.2813795 로서 유의수준 0.05보다 크므로 귀무가설 H0 '두 집단 간 평균차이는 없다'를 채택하게 됩니다. 이는 저 위에 있는 Boxplot을 보면 좀더 쉽게 이해가 될 것 같습니다.
위 분석 결과를 보기에 좋게 표로 나타내면 아래와 같습니다.
|
pair-wise comparison |
group 1 |
group 2 |
group 3 |
|
group 1 |
- |
mean diff. -4.51 ***
(95% confidence interval
: -5.45 ~ -3.58) |
mean diff. -5.11 ***
(95% confidence intervl
: -6.04 ~ -4.17) |
|
group 2 |
|
- |
mean diff. -0.59
(95% confidence interval
: -1.52 ~ 0.34) |
|
group 3 |
|
|
- |
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey's HSD test의 이론에 대해서는 알쏭달쏭한 말들을 잔뜩 써놨는데요, R 함수는 TukeyHSD() 딱 한줄이어서 미친 듯이 간단합니다. ^^' 다음번 포스팅에서는 Duncan's LSR(least significant range) test에 대해서 알아보도록 하겠습니다.
많은 도움 되었기를 바랍니다.
☞ one-way ANOVA 참고
☞ 다중비교(multiple comparison) : Duncan's LSR(least significant range) test 참고
☞ 대비(contrast) : Scheffe test 참고
이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^
