R, Python 분석과 프로그래밍의 친구 (by R Friend)

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 3개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 3개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다. 

지난번 포스팅에서 이원분산분석(two-way ANOVA) 중에서 (1) 관측값이 하나일 경우 (one observation in each cell)의 이원분산분석 에 대해서 알아보았다면, 이번 포스팅에서는 (2) 관측값이 두개 이상일 경우(more than 2 observations in each cell)의 이원분산분석에 대해서 알아보겠습니다.  단, 이때 각 집단 내 관측값의 개수는 동일합니다.

관측값이 두개 이상일 경우의 이원분산분석에서는 두 개의 요인 수준별 주 효과(main effect)와 더불어 두 요인이 서로 상호간에 영향을 주고 받으면서 나타나는 반응효과인 교호작용 효과(interaction effect)를 추가로 분석하는 것이 관측값이 하나인 경우와의 차이점입니다.

[ 관측값의 개수에 따른 이원분산분석의 분석 대상 ]

요인 A의 i번째 수준과 요인 B의 j번째 수준에서 측정된 k번째 관측값을 Yijk라고 할 때 데이터셋의 형태는 아래와 같습니다. 

[ 이원분산분석 - 데이터 형태 ]

(dataset for two-way ANOVA, factor A & B, factor levels a & b, observations k)]

요인 A의 i번째 수준과 요인 B의 j번째 수준에서 측정된 k번째 관측값을 Yijk는 요인 A와 B의 주 효과(main effect)와 함께 두 요인의 교호작용 효과(interaction effect)이 존재하게 되며, 편차 는 아래와 같이 4개의 성분합으로 나타낼 수 있습니다.

[ 편차의 4개 성분합으로의 분해 ]

위의 식의 양변을 제곱하여 모든 관측값에 대하여 더한 후에 정리를 하면

의 관계식을 얻게 됩니다.

이를 요인과 교호작용별로 제곱합, 자유도, 평균제곱, F통계량을 알기 쉽도록 정리한 이원분산분석표는 아래와 같습니다.

[ 이원분산분석표 (two-way ANOVA table) ]

 * SST : Total Sum of Squares,  SSE : Error Sum of Squares

검정통계량으로는 F 통계량(요인A 효과 검정 = MSA/MSE, 요인B 효과 검정 = MSB/MSE, 교호작용 효과 검정 = MSAB/MSE)을 사용하며, 요인 A의 효과와 요인 B, 교호작용 효과에 대한 검정 절차 및 방법은 아래와 같습니다.

[ 이원분산분석 검정 방법 ]

그럼, 아래 예제에 대해서 R의 aov() 함수를 활용해서 문제를 풀어보도록 하겠습니다.

아래에 요약통계량과 교호작용 여부를 가늠해볼 수 있는 그래프 생성 R 함수도 일부 추가하였습니다. 만약 교호작용효과 그래프(Interaction Effect Plot) 가 서로 교차를 하면 교호작용이 있다고 보면 되며, 서로 평행선을 이룬다면 교호작용이 없다고 해석하면 되겠습니다.  (교호작용도는 엑셀로 치면 "꺽은선형 차트"가 되겠습니다)

[ R aov() 함수 사용 방법 ]

> # (2) two-way ANOVA when there are more than one observation per cell (different treatment groups)
> #     (the number of observations in each cell must be equal)
> 
> gender.fac <- as.factor(c(rep("M", 9), rep("F", 9)))
> gender.fac
 [1] M M M M M M M M M F F F F F F F F F
Levels: F M
> 
> class <- c("class_1", "class_1", "class_1", "class_2", "class_2", "class_2", "class_3", "class_3", "class_3")
> class.fac <- as.factor(c(rep(class, 2)))
> class.fac
 [1] class_1 class_1 class_1 class_2 class_2 class_2 class_3 class_3 class_3 class_1 class_1
[12] class_1 class_2 class_2 class_2 class_3 class_3 class_3
Levels: class_1 class_2 class_3
> 
> score_stats <- c(71, 77, 78, 76, 77, 78, 71, 70, 69, 80, 76, 80, 79, 78, 77, 73, 71, 70)
> 
> 
> # summary statistics
> score.df <- data.frame(gender.fac, class.fac, score_stats)
> score.df
   gender.fac class.fac score_stats
1           M   class_1          71
2           M   class_1          77
3           M   class_1          78
4           M   class_2          76
5           M   class_2          77
6           M   class_2          78
7           M   class_3          71
8           M   class_3          70
9           M   class_3          69
10          F   class_1          80
11          F   class_1          76
12          F   class_1          80
13          F   class_2          79
14          F   class_2          78
15          F   class_2          77
16          F   class_3          73
17          F   class_3          71
18          F   class_3          70
> 
> install.packages("doBy")

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/doBy_4.5-13.zip'
Content type 'application/zip' length 3431466 bytes (3.3 MB)
downloaded 3.3 MB

package ‘doBy’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpWOKfji\downloaded_packages
 
 
> library(doBy)
> summaryBy(score_stats ~ gender.fac, data=score.df, FUN = c(mean, sd, min, max))
  gender.fac score_stats.mean score_stats.sd score_stats.min score_stats.max
1          F         76.00000       3.807887              70              80
2          M         74.11111       3.756476              69              78
> summaryBy(score_stats ~ class.fac, data=score.df, FUN = c(mean, sd, min, max))
  class.fac score_stats.mean score_stats.sd score_stats.min score_stats.max
1   class_1         77.00000       3.346640              71              80
2   class_2         77.50000       1.048809              76              79
3   class_3         70.66667       1.366260              69              73
> summary(score_stats, data=score.df)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  69.00   71.00   76.50   75.06   78.00   80.00 
> 

> # box plot and interaction plot
> par(mfrow = c(2, 2))
> plot(score_stats ~ gender.fac, main="box plot by gender")
> plot(score_stats ~ class.fac, main="box plot by class")
> interaction.plot(gender.fac, class.fac, score_stats, bty='l', main="interaction effect plot")
> interaction.plot(class.fac, gender.fac, score_stats, bty='l', main="interaction effect plot")
> 

> 
> # two-way ANOVA : aov()
> # replicates, interaction effect
> aov_model = aov(score_stats ~ gender.fac + class.fac + gender.fac:class.fac)
> summary(aov_model)
                     Df Sum Sq Mean Sq F value   Pr(>F)    
gender.fac            1  16.06   16.06   3.853 0.073249 .  
class.fac             2 174.11   87.06  20.893 0.000123 ***
gender.fac:class.fac  2   4.78    2.39   0.573 0.578354    
Residuals            12  50.00    4.17                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R 분석 결과 중에서 제일 아래에 있는 이원분산분석표를 가지고 해석을 해보겠습니다.  먼저 (1) 교호작용이 있는지를 살펴보면, 성별과 학급별에 따른 교호작용효과에 대한 P-value가 0.578354로서 유의수준 α 0.05보다 크므로 귀무가설 H0를 채택하여 "성별과 학급의 교호작용에 의한 효과는 없다"고 판단할 수 있습니다.

두번째로 (2) 성별 요인(factor A levels of "M", "F")에 대한 P-value는 0.073249로서 유의수준 α 0.05 보다 크므로 귀무가설 H0를 채택하게 되어 "성별에 따른 통계학 성적 차이는 없다"고 판단할 수 있습니다.

세번째로 (3) 학급 요인(factor B levels of "class_1", "class_2", "class_3")에 따른 P-value는 0.000123으로서 유의수준 α 0.05보다 작으므로 귀무가설 H0를 기각하고 대립가설 H1을 채택하여 "학급에 따른 통계학 성적의 차이가 있다"고 판단할 수 있습니다.

학급 요인에 대해서는 유의수준 α 0.05에서 차이가 있다고 나왔으므로, 어떤 학급 간에 차이가 있는지를 알아보고 싶다면, 사후 쌍을 이룬 다중 비교 분석(post-hoc pair-wise multiple comparison)을 하면 됩니다.  이에 대한 자세한 설명은 아래의 링크를 참조하시기 바랍니다.

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

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

2개의 모집단에 대한 평균을 비교, 분석하는 통계적 기법으로 t-Test를 활용하였다면, 비교하고자 하는 집단이 3개 이상일 경우에는 분산분석 (ANOVA : Analysis Of Variance)를 이용합니다. 

설명변수는 범주형 자료(categorical data)이어야 하며, 종속변수는 연속형 자료(continuous data) 일 때 3개 이상 집단 간 평균 비교분석에 분산분석(ANOVA) 을 사용하게 됩니다.

분산분석(ANOVA)은 기본적으로 분산의 개념을 이용하여 분석하는 방법으로서, 분산을 계산할 때처럼 편차의 각각의 제곱합을 해당 자유도로 나누어서 얻게 되는 값을 이용하여 수준평균들간의 차이가 존재하는 지를 판단하게 됩니다.  

[ 일원분산분석 (one-way ANOVA)의 개념 ]

분산분석 (ANOVA)은 영국의 통계학자 피셔(R.A. Fisher)가 고안한 분석 기법으로서, 최초에는 농사실험과 관련된 연구에서 사용되었으나 요즘에는 사회과학, 공학, 의학 등 폭넓게 적용되고 있습니다.

[ 분산분석 (ANOVA)의 분류 ]

3개 이상의 집단에 대한 비교 시 모집단이 정규분포를 따르지 않을 경우에는 비모수 추정과 검정 기법을 이용해야 하며, 이에 대해서는 추후에 별도로 다루도록 하겠습니다.

분산분석은 측정하고자 하는 값에 영향을 미치는 요인(Factor)의 수에 따라서 구분하는데요, 가령, 작업자(Worker)에 따라서 생산량(Volume of Production)의 차이가 있는가를 비교 분석한다고 했을 때 '작업자(Worker)'가 바로 요인(Factor)이 되겠습니다. 그리고 작업자 3명 '홍길동', '김철희', '최영희'를 요인 수준 (Factor Level) 혹은 처리 (treatment) 라고 합니다.

측정값에 영향을 미치는 요인(Factor)이 1 개인 실험은 '일원분산분석' (one-way ANOVA) 라고 하며, 측정값에 영향을 미치는 요인(Factor)이 2 개인 실험은 '이원분산분석' (two-way ANOVA) 라고 부릅니다.

[ ANOVA Model and Description in R ]

이번 포스팅에서는 요인(Factor)이 1개인 일원분산분석(one-way ANOVA)의 이론과 R의 aov() 함수의 사용법에 대해서 설명해보겠습니다.

요인(Factor)이 1개이면서 r개의 요인 수준(Factor Level)을 가지고 있고, 각 수준에서 박복 측정수가 n개인 일원분산분석(one-way ANOVA)는 다음과 같은 형태의 자료 형태와 모형을 가지게 됩니다.

[ 일원분산분석 데이터 형태 (dataset for one-way ANOVA) ]

[ 일원분산분석 모형 (one-way ANOVA model) ] 

분산분석(ANOVA)은 측정된 자료값들의 전체 변동을 비교하고자 하는 요인 수준(Factor Level) 간의 차이에 의해서 발생하는 변동과 그 밖의 요인에 의해서 발생하는 변동으로 나누어서 자료를 분석하는 것이 기본 원리가 되겠습니다.

[ 분산분석 기본 원리 1 ]

위에서 제시한 식의 양변을 제곱하여 전체 측정값들에 대해서 모두 더하여 식을 정리하면 아래와 같이 됩니다. (자세한 과정은 생략함)

[ 분산분석 기본 원리 2 ]

위의 총제곱합(SST)와 처리제곱합(SSTR), 오차제곱합(SSE) 간의 기본원리를 이용하여 분산분석에서 사용하는 통계량들을 표로 일목요연하게 정리한 것이 아래의 일원분산분석표(one-way ANOVA table) 가 되겠습니다.

[ 일원분산분석표 (one-way ANOVA table) ]

분산분석표의 제일 오른쪽에 있는 F0 통계량을 이용해서 전체 수준들간의 평균이 같은지 아니면 다른지를 검정합니다.  기본 개념을 설명하자면, F0 통계량은 처리평균제곱 (MSTR)의 크기에 영향을 받으므로 처리평균제곱 (MSTR)이 커지면 오차평균제곱(MSE)은 작아지게 되며, 따라서 F0 통계량은 분자가 커지고 분모가 작아지므로 당연히 커지게 됩니다. 즉, F0 통계량 값이 크다는 것은 수준 평균들간에 평균의 차이가 크다는 의미가 되겠습니다.

그러면, 요인효과에 대한 검정을 위해서 분산분석에서는 아래와 같은 귀무가설과 대립가설을 사용하며, 검정통계량으로는 F 를 사용하여 기각역 또는 P-value 에 의해서 검정을 하게 됩니다.  

이론 설명이 무척 길어졌는데요, 이제 드디어, 아래 문제를 R의 aov() 함수를 사용해서 풀어보도록 하겠습니다.

문제 ) 정유사에서 온도(Factor, one-way)에 따라서 휘발유 생산량에 변화가 있는지 (즉, 영향이 있는지) 알아보기 위하여 온도를 3가지 조건(3 Factor Levels)으로 실험설계를 하여 10번에 걸쳐서 휘발유 생산량을 측정하였습니다. 관찰값이 아래와 같을 때 조사되었을 때 온도의 조건에 따라서 휘발유 생산량에 차이가 있는지 유의수준 α = 10% 로 검정하시오.

R에 (1) 위의 관측값들을 벡터로 입력하고,

      (2) boxplot() 함수와 summary()를 이용하여 탐색적 분석을 해본 후에,

      (3) aov() 함수를 이용하여 one-way ANOVA 분석을 하였으며 
         : aov(y ~ group, data = dataset)

      (4) 기본 가정 중에 오차의 등분산성을 검정하기 위해 Bartlett 검정

         : bartlett.test(y ~ group, data = dataset)

을 실시하였습니다.

이때 조심해야 할 것이 있는데요, dataset을 데이터프레임으로 해서 그룹(group, factor level)에 해당하는 변수는 반드시 요인(Factor)형이어야 하므로, 만약 요인(Factor)형이 아니라면 먼저 transfrom(factor()) 함수를 이용해서 변환을 시켜주고 ANOVA 분석을 수행해야 합니다.

> ##---------------------------------------------------------- 
> ## 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  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

일원분산분석(one-way ANOVA) 결과 검정통계량 F-value가 108.8으로 나왔으며, P-value 값은 '1.2e-13'으로서 매우 작게 나와 유의수준 10% 에서 귀무가설을 기각하고 대립가설을 채택하게 되어 "온도 조건 1/2/3에 따라서 휘발유 생산량에 차이가 있다"라고 판단할 수 있겠습니다.

 오차의 등분산성 가정에 대해 Bartlett 검정 결과 P-value가 0.4368로서 유의수준 10%에서 귀무가설을 채택하여 "오차의 등분산성 가정을 만족한다"고 할 수 있겠습니다.

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

질문은 댓글로 남겨주세요.

• 쌍을 이룬 집단 간 평균 다중비교 (multiple comparisons)

☞ Tukey's HSD(honestly significant difference) test 참조

☞ Duncan's LSR(least significant range) test 참고

• 대비 (contrast)

☞ 샤페 검정법 (scheffe test) 참고

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