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

 

 

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

 

728x90
반응형
Posted by Rfriend
,