R (2) 다중비교 - Tukey의 HSD (honestly significant difference) 검정 : TukeyHSD()

지난번 포스팅에서 일원분산분석(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 참고

 

 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡'를 꾸욱 눌러주세요. ^^