지난번 포스팅에서는 ANOVA 검정 시 집단 간 평균에 유의미한 차이가 있다고 판단될 때 사후 분석으로 다중비교((post-hoc pair-wise multiple comparison)을 하는 방법 중에 Tukey's HSD(honestly significant difference) test에 대해서 알아보았습니다.

 

이번 포스팅에서는 다중비교(multiple comparison)을 하는 두번째 방법으로 Duncan's LSR(Least Significant Range) test에 대해서 알아보도록 하겠습니다.

 

 

 

[ 사후검정 중비교(post-hoc pair-wise multiple comparison) ]

 

 

 

 

 

 

기본 개념은 짝을 이룬(pair-wise) 수준평균간의 차이가 아래에 제시된 Duncan의 최소유의범위(Least Significant Range) 보다 크다면 "유의수준 α에서 두 수준평균간에는 차이가 있다"고 판단하겠다는 것입니다.

 

 

 

[ Duncan's LSR(least significant range) test procedure ]

 

 

Duncan's LSR test는 앞서 포스팅했던 Tukey's HSD test보다도 약간 더 보수적입니다.  그래서 Tukey's HSD test에서는 차이가 없다고 나오는 짝을이룬 수준평균간에도 Duncan's LSR test에서는 차이가 있다고 판정하는 경우도 있습니다.

 

R을 가지고 Duncan's LSR test를 하는데 있어 

 - duncan.test() function of agricolase package

 - LDuncan() function of laercio package

두가지 함수를 차례로 사용해보겠습니다.

 

온도의 조건 3가지에 따른 생산량을  각10번씩 측정한 아래의 데이터를 가지고 (1) ANOVA test를 먼저 해보고, (2) 사후 분석으로 Duncan's LSR test를 위의 2가지 패키지의 함수를 가지고 각각 실시해보겠습니다.

 

 

> ##---------------------------------------------------------- > ## multiple comparison > ## - Duncan's LSR(least significant range) 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" > > > # (1) ANOVA test > aov_model <- aov(y ~ group, data = group_df) > summary(aov_model) 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 > > > # (2-1) Duncan's LSR(least significant range) test

> # with duncan.test() function of agricolase package > install.packages("agricolae")

Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/agricolae_1.2-3.zip'
Content type 'application/zip' length 921635 bytes (900 KB)
downloaded 900 KB

package ‘agricolae’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpYhI0Fo\downloaded_packages

 

> library(agricolae) > duncan.test(aov_model, "group", alpha = 0.05, console = TRUE) Study: aov_model ~ "group" Duncan's new multiple range test for y Mean Square Error: 0.7182593 group, means y std r Min Max 1 51.46 0.6834553 10 50.5 52.4 2 46.94 1.0415800 10 45.0 47.9 3 46.35 0.7763876 10 44.9 47.2 alpha: 0.05 ; Df Error: 27 Critical Range 2 3 0.7776731 0.8170522 Means with the same letter are not significantly different. Groups, Treatments and means a 1 51.46 b 2 46.94 b 3 46.35 > >

>

> # (2-2) Duncan's LSR(least significant range) test

> # with LDuncan() function of laercio package > install.packages("laercio")

Installing package into ‘C:/Users/user/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
trying URL 'http://cran.rstudio.com/bin/windows/contrib/3.2/laercio_1.0-1.zip'
Content type 'application/zip' length 20726 bytes (20 KB)
downloaded 20 KB

package ‘laercio’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpYhI0Fo\downloaded_packages

 

> library(laercio) > aov_model <- aov(y ~ group, data = group_df) > LDuncan(aov_model, "group") DUNCAN TEST TO COMPARE MEANS Confidence Level: 0.95 Dependent Variable: y Variation Coefficient: 1.75648 % Independent Variable: group Factors Means 1 51.46 a 2 46.94 b 3 46.35 b Restarting R session...

 

결과를 해석해보면, (1) ANOVA 검정 결과 P-value가 1.2e-13 로서 유의수준 0.05보다 매우 작으므로 "집단간 평균에 차이가 있다"는 대립가설을 채택할 수 있습니다.

 

(2) 이에 사후분석으로 "어느 집단 간에 차이가 있는가?"를 알아보기 위해 Duncan's LSR test 를 해보니,

온도의 조건에 따른 3가지 factor level 별 생산량의 평균은 아래와 같았고, factor level '2'와 '3'은 서로 차이가 없고, factor level '1'만 짝을 이룬 수준평균간에 차이가 있다고 나왔습니다.

 

 factor level

mean 

 group

significantly different?

 1

 51.46

 a ***

 different

 2

 46.94

 b

 not different

 3

 46.35

 b

 

 

차이가 있다고 판단한 근거를 보면, MSE(mean squared error) 가 0.7182593 이고, 이 값을 위에서 제시했던 Duncan's LSR(Least Significant Range)를 계산하는 공식에 넣어서 계산을 해보면, 위에 R분석결과에 Critical Range 라고 해서 LSR2 = 0.7776731, LSR3 = 0.8170522 가 나왔습니다.  이 값들보다 수준집단간 평균의 차이가 더 크면 차이가 있다고 판단하게 되는 것입니다. 

 

3개 수준평균을을 크기 순서대로 배열한 후에 각 차이를 구한 것을 표로 제시해보면 아래와 같으며, 이 값들이 LSR2 = 0.7776731, LSR3 = 0.8170522 보다 크면 *** 표시를 했습니다.

(LSR2 는 인접한 수준들간의 평균 차이를 의미, LSR3는 2단계 떨어진 수준들간의 평균 차이를 의미)

 

 

mean of factor level 3: 

46.35

mean of factor level 2:

46.94

mean of factor level 1:

51.46 

mean of factor level 3:

 46.35

 0.59

5.11 ***

mean of factor level 2: 

46.94

 

4.52 ***

mean of factor level 1:

 51.46

 

 

 

많은 도움 되었기를 바랍니다.

 

☞ one-way ANOVA 참고

☞ 다중비교(multiple comparison) Tukey's HSD(honestly significant difference) test 참고

☞ 대비(contrast) Scheffe test 참고

 

 

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

 

 

Posted by R Friend R_Friend

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

 

 

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

 

Posted by R Friend R_Friend