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)이 존재하게 되며, 편차 Yijk - Y^... 는 아래와 같이 4개의 성분합으로 나타낼 수 있습니다.

 

 

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

 

 

 

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

 

SST = SSA + SSB + SSAB + SSE 

 

의 관계식을 얻게 됩니다.

 

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

 

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

 

 요인

제곱합

(squared sum)

자유도

(degrees of freedom) 

평균제곱

(mean squared) 

F statistics 

 요인 A

(factor A)

 SSA

a-1

MSA

MSA/MSE

 요인 B

(factor B)

 SSB

b-1

MSB

MSB/MSE

 교호작용

(interaction effect of A, B)

 SSAB

 (a-1)(b-1)

MSAB

MSAB/MSE 

 오차

(error)

 SSE

ab(n-1)

MSE

 

 계

(total)

 SST

 nab-1

 

 

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

 

 

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

 

 

[ 관측값이 두개 이상인 경우의 이원분산분석 절차 (procedure of two-way ANOVA when there are more than 2 observations in each cell) ]

 

 

 

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

 

 

(1) 교호작용 효과 (interaction effect)

   - 귀무가설 H0 : (αβ)ij = 0

   - 대립가설 H1 : 모든 (αβ)ij는 0이 아니다

   - 검정통계량 : F0 = MSAB/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

(2) 요인 A의 주 효과 (main effect of factor A)

   - 귀무가설 H0 : α1 = ... = αa = 0

   - 대립가설 H1 : 모든 αi 는 0이 아니다

   - 검정통계량 : F0 = MSA/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

(3) 요인 B의 주 효과 (main effect of factor B)

   - 귀무가설 H0 : α1 = ... = αb = 0

   - 대립가설 H1 : 모든 αj 는 0이 아니다

   - 검정통계량 : F0 = MSB/MSE

   - 판정 : P-value가 유의수준(α)보다 작으면 귀무가설 H0를 기각하고, 대립가설 H1을 채택

 

 

 

 

 

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

 

 

(문제)  학급과 성별에 따른 통계학 성적이 아래와 같다고 할 때(각 cell 별 3명의 학생 성적 측정), 유의수준 α = 0.05 에서 학급과 성별 요인의 주효과와 교호작용효과에 대해 검정하시오.

 

[학급과 성별에 따른 통계학 점수표]

 

         Factor B (class)

 

Factor A (gender)

Class 1 

Class 2 

Class 3 

Mean 

남성 (M) 

 71, 77, 78

76, 77, 78 

71, 70, 69 

 75.33

 여성 (F)

 76, 76, 80

79, 78, 77

71, 71, 70 

 74.11

 Mean

76.33 

77.50 

70.33 

74.72 

 

 

 

아래에 요약통계량과 교호작용 여부를 가늠해볼 수 있는 그래프 생성 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)을 하면 됩니다.  이에 대한 자세한 설명은 아래의 링크를 참조하시기 바랍니다.

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

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

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

 

 

  • 대비 (contrast)

샤페 검정법 (scheffe test) 참고

 

 

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

 

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

 

Posted by R Friend R_Friend

댓글을 달아 주세요

  1. Kwon 2016.04.06 00:47  댓글주소  수정/삭제  댓글쓰기

    R 과 관련된 포스팅 지난주부터 쭉 정독해서 보고 있습니다. 제게 큰 도움이 된 것 같습니다. 정말 감사합니다 ^^

    혹시, 비모수 일원분산분석 / 비모수 이원분산분석 / 회귀분석 등과 관련한 포스팅 도 해주실 수 있으신지요 ?

    • R Friend R_Friend 2016.04.06 06:38 신고  댓글주소  수정/삭제

      좋게 봐주셔서 감사합니다. ^^ 비모수분석쪽은 포스팅 몇개 써놓은게 있으니 참고하시구요, 회귀분석은 앞으로 쓰긴 할텐데요...일단 요즘 연재하고 있는 선형대수 먼저 끝내고 쓰려고 합니다. 직장인이다보니 주말에만 글쓰곤 해서 시간이 오래 걸리네요

  2. 2016.05.22 18:17  댓글주소  수정/삭제  댓글쓰기

    rm anova라고 하려면, 한 개체에 대해 반복적으로 측정하는 변수가 있어야 합니다. 예를 들자면 환자 한명에 대해 매달 1회씩 5달에 거쳐 측정을 한다면 시간이란 변수가 반복측정변수가 됩니다.. 여하간에 잘못된 포스팅으로 보입니다

  3. 통계병아리 2016.10.25 14:52  댓글주소  수정/삭제  댓글쓰기

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    분산분석 시행시 결과 출력물 맨 아래에 이 문장이 뜨는데, 이것을 어떻게 분석하면 되나요? 아니면 분석할 필요가 없는 단순한 index 인가요?

    그리고, class <- c("class_1", "class_1", "class_1", "class_2", "class_2", "class_2", "class_3", "class_3", "class_3")를
    class <- c(rep("class_1",3),rep("class_2",3),rep("class_3",3)) 으로 나타내서 풀어보았는데
    factor화 하니 숫자로 나오더군요, 왜 이러한 결과가 나오는 걸까요?

    • R Friend R_Friend 2016.10.25 17:24 신고  댓글주소  수정/삭제

      1) Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      는 분석결과를 해석할 때 보기에 편하라고 매 분석때마다 동일하게 나타나는 안내문입니다. 회귀분석을 해보면 각 x변수별로 p-value 가 0.001 보다 작으면 별이 3개(***), p-value가 0.001~0.01 이면 별이 2개(**), p-value가 0.01~0.05 이면 별이 1개(*), p-value가 0.05보다 크면 별이 없게 보여줍니다. p-value 숫자를 일일이 안봐도 별 개수만 보면 빠르고 편하게 해석할 수 있습니다.

      2) 두번째 질문은 저도 왜 그런지 모르겠네요. rep() 함수를 사용해서 해도 결과가 동일할텐데요... character 를 factor로 변환하는거라 에러가 끼어들 여지가 없어보이는데요... 왜 그런지 저도 모르겠습니다. ^^;

    • 통계병아리 2016.10.25 21:45  댓글주소  수정/삭제

      아하, 그런 것이었군요. 이제 이해가 되네요.
      rep 함수 문제는 한 번 다시 해봐야겠네요, 감사합니다.

  4. 통계왕 2017.01.26 18:05  댓글주소  수정/삭제  댓글쓰기

    R을 이용한 글 정말정말 잘보고있습니다.!!!
    제가 궁금한건 aov를 적용할때 + 와 * 의 차이를 여쭤보고 싶습니다.
    +만을 쓸 경우 교호작용이 출력되지 않고(trt:trt를 쓰지 않는다는 가정에서)
    *를 쓰게되면 교호작용이 함께 출력되는데 이 차이 말고는 똑같은건가요 ?
    그럼 굳이 trt:trt 이런방법을 쓰지 않아도 되지 않는가 궁금합니다.!!
    그리고 anova의 자유도에 대해서도 궁금합니다 ㅜ ㅜ (사실 이게 더 궁금합니다....)
    교호작용여부를 출력하고 싶은데 "Estimated effects may be unbalanced" 이러한
    문구와 함께 p value 가 출력되지 않는 경우가 있습니다. 구글링을 해보면 자유도의 문제라고
    하는것 같은데 정확한 해답을 얻지 못했습니다..ㅠ

    새해 복 많이 받으시고 앞으로도 좋은 R 글 기대하겠습니다. !!

    • R Friend R_Friend 2017.01.26 18:29 신고  댓글주소  수정/삭제

      통계왕님, 반갑습니다.

      (1) R에서 'x*y' notation은 'x', 'y', 'x와 y의 교호작용변수' 의 3개를 모두 포함한다는 뜻입니다.

      반면 'x:y' notation 은 'x와 y의 교호작용변수' 만을 뜻합니다.

      본문에서 x + y + x:y 라고 쓴 것은 x변수 주효과, y변수 주효과, x와 y변수의 교호작용효과 각 각을 직접 하나씩 명기해줬다고 보시면 됩니다.

      (2) 'Estimated effects may be unbalanced' 메시지는 ANOVA 분석 대상 cell 의 관측치 개수가 동일하지 않을 때 뜹니다.

      이원분산분석은 "각 cell의 관측치 개수가 동일"할 때 사용합니다. 본문의 예에서는 각 cell 별로 3개씩의 동일한 관측치가 들어가 있습니다.

      cell 관측치 개수 확인해보시기 바랍니다.

    • 통계왕 2017.01.26 19:34  댓글주소  수정/삭제

      아 그럼.. 단순히 데이터프레임의 오류인건가요??? 자유도와는 상관없는 내용인건가요,,,?!

  5. SMA 2017.08.10 14:33  댓글주소  수정/삭제  댓글쓰기

    안녕하세요
    R 공부 하며서 정말 큰 도움을 받고 있어요. 감사합니다. ^-----^

    이 포스팅을 읽고

    RM ANOVA 로 분석을 하고 자 할때는
    aov를 어떻게 적용해야할 지에 대해 여쭤보고 싶어서 댓글드려요

    감사합니다. ~~

    • R Friend R_Friend 2017.08.10 15:34 신고  댓글주소  수정/삭제

      안녕하세요 .
      Repeated Measure 데이터도 포스팅 R 예제와 동일한 장식으로 anova 분석하시면 됩니다.

      데이터 형태가 똑같습니다.

    • SMA 2017.08.13 16:41  댓글주소  수정/삭제

      답변 감사합니다. ^----^

      전 '객' 님의 댓글에 대한 답변글을 보고
      '바로 잡았다'고 하는 부분이 헷갈렸어요

      궁금한 부분을 좀더 자세히 말한다면
      다음과 같아요

      저는 특정 치료를 받은 군과 치료를 받지 않은 군으로 나누었고, 특정 test를 '실험전, 실험중간, 실험후'에 반복측정하여 이 데이터를 RM ANOVA로 돌려
      치료여부효과, 시간효과 , 치료와 시간의 교호작용을 보려했어요


      위의 포스터에서 예제로 올려주신 데이터는 한개체에서 시간에 따른 반복측정을 한 것이 아니어서
      동일하게 적용해도 오류가 없을까 의문이 들었어요

      ( 예를 들면 위의 예제
      성별*class의 데이터 값 각각이 아무런 연관성이 없지만

      치료*시간의 데이터 값 즉
      한사람에서 시간에 따라 3번 측정한 데이터 값은 처음 데이터 값이 두번째, 마지막에 측정한 데이터 값에 영향을 미쳐서요.)


      답변처럼 동일하게 aov를 적용해도 될것 같기도 하고, 오류가 있을 것 같기도 해서 한번 더 질문드려요 ^^

      감사합니다.

  6. 2018.05.03 09:21  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

  7. 임카터 2018.05.03 09:23  댓글주소  수정/삭제  댓글쓰기

    의학 관련 댓글 남겼는데 혹시 댓글 보는 방법을 몰라서 혹시라도 보셨더라면, carter.lim09@gmail.com 이쪽으로 한번만 확인 해주시면 진짜 감사하겠습니다. ㅠㅠ !!

    좋은 하루 보내세요! 오늘은 날씨가 제 마음과 같네요. ㅎㅎ

  8. 닥터케빈 2018.10.23 11:21  댓글주소  수정/삭제  댓글쓰기

    안녕하세요? 항상 유익한 포스팅으로 공부 잘하고 있습니다.
    한가지 궁금한 것이 있는데요. 편차의 분해에서 교호작용과 오차 부분이 조금 이상한 것 같습니다.
    제 생각에는 교호작용 부분은 y_ij.의 평균 - y_i..의 평균 - y_.j.의 평균 + y_...의 평균이 될 것 같고,
    오차는 y_ijk - y_ij.의 평균될 것 같습니다.

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 ]

 

n-way ANOVA

Model

Description

one-way ANOVA

y ~ x1

  y is explained by x1 only

two-way ANOVA

y ~ x1 + x2

  y is explained by x1 and x2

two-way ANOVA

y ~ x1 * x2

  y is explained by x1, x2 and the interaction   
  between them

three-way ANOVA

y ~ x1 + x2 + x3

  y is explained by x1, x2 and x3

 

 

 

이번 포스팅에서는 요인(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) ]

 

 구분

 제곱합

(Sum of Squares)

자유도

(degrees of freedom) 

평균제곱

(Mean Square Error) 

검정통계량 F0

(F statistics)

 처리

(Treatment)

 SSTR

r-1 

MSTR 

MSTR/MSE 

 오차

(Error)

 SSE

 nT - r

 MSE

 

 전체

(Total)

 SST

 nT - 1

 

 

 

 

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

 

 

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

 

 

  [ 가설 ]

   귀무가설 H0 : μ1 = μ2 = ... = μr

   대립가설 H1 : 모든 μi 는 같지 않다.  i = 1, 2, ..., r

 

  [ 기각역 혹은 P-value에 의한 검정 ]

   표본으로부터 계산된 검정통계량의 값 f0가 유의수준(significance level) α 에서

 

   f0 >= Fα(r - 1, nT - r) 또는 P-value = P(F >= f0) < α 이면, H0 기각 (H0 reject)

 

   f0 < Fα(r - 1, nT - r) 또는 P-value = P(F >= f0) >= α 이면, H0 채택 (H0 accept)

 

 

 

 

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

 

 

 


 

 

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

 

 

 

 10번 실험에 의한 측정값

 요인 수준

(Factor Level)

 1

2

 10

 온도 조건 1

(Temp condition 1)

 50.5

52.1

51.9

52.4 

50.6

51.4 

51.2 

52.2 

51.5 

50.8 

 온도 조건 2

(Temp condition 2)

 47.5

47.7 

46.6 

47.1 

47.2 

47.8 

45.2 

47.4 

45.0 

47.9 

 온도 조건 3

(Temp condition 3)

 46.0

47.1

45.6

47.1

47.2

46.4

45.9 

47.1

44.9

46.2 

 

 

 


 

 

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

 

 

 

일원분산분석(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) 참고

 

 

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

 

 

 

Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 2016.03.14 14:49  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend R_Friend 2016.03.14 21:15 신고  댓글주소  수정/삭제

      하루에 400여명 방문자가 가운데 대부분의 분들은 그냥 포스팅 글만 보고 가는데요, 이렇게 도움이 되었다는 댓글 남겨주시니 감사합니다. ^^

  2. 쿡북 2016.08.05 09:42 신고  댓글주소  수정/삭제  댓글쓰기

    아 진짜 잘 봤습니다. 이렇게 R 코드 레벨까지 깔끔하게 정리해주시다니...
    다른 글들도 지금 챙겨보는 중입니다. 성의 있는 포스팅 감사합니다.

  3. 쿡북 2016.08.09 21:10 신고  댓글주소  수정/삭제  댓글쓰기

    네 글좀 페북에 공유차 올렸습니다 아마 많이들 들어오실듯 ^^ 좋은글 감사합니다

    • R Friend R_Friend 2016.08.09 21:26 신고  댓글주소  수정/삭제

      JY_P님, 페북 링크해주셔서 감사합니다. 어제 580명 방문이었는데 오늘 현재 2,300여명 방문이네요. 역대 최고 기록입니다. 페북 영향력이 굉장하시네요. ^^b

  4. 쿡북 2016.08.09 21:28 신고  댓글주소  수정/삭제  댓글쓰기

    ㅎㅎ 그렇죠? 좋은글은 공유가 잘 되어야죠 ㅎㅎ 혹시 페북하시면 여러 활동하셔도 좋을듯 하네요 앞으로 많이 드나들겠습니다 감사합니다

  5. 통계야친구하자 2016.08.09 22:05  댓글주소  수정/삭제  댓글쓰기

    통계의 필요성을 막 깨닫고 여기저기 기웃거리고 있는 일인입니다. 페이스북에 개념글로 소개된것 보고 여기까지 왔네요ㅎ정성 깊은 포스팅 ,개념정립 하는데 많은 도움 받겠습니다. 감사합니다

  6. !! 2016.08.09 22:14  댓글주소  수정/삭제  댓글쓰기

    페북에서 보고 왔는데. 와우! R뿐 아니라 통계설명도 너무 훌륭하네요! 널리 알려야겠습니다.

  7. freepsw 2016.08.10 00:20  댓글주소  수정/삭제  댓글쓰기

    JY_P 와 R_Friend 두분은 서로 아는사이인데 ㅎㅎ
    모른척 하시는건가요 ^^
    이차장님 일도 바쁘신데 블로그까지 대단하세요!!

  8. sfsfsf 2016.08.10 23:02  댓글주소  수정/삭제  댓글쓰기

    감사합니다!
    마지막 예제에서 anova table에서 궁금한 것이 있습니다.
    SSTR의 d.f가 어째서 1이 되는거죠? factor level이 1,2,3 이므로 3-1=2 여야 하지 않나요??
    SSE의 d.f는 30-3=27 이고
    합해서 Total의 d.f가 29가 되어야 하지 않나 하는데 왜 이렇게 출력되는 건지 궁금합니다.

    • R Friend R_Friend 2016.08.10 23:21 신고  댓글주소  수정/삭제

      sfsfsf님,
      aov(y ~ group, data = group_df) 에서는 말씀하신대로 제대로 분석결과가 나왔는데요,

      summary(aov(y ~ group, data = group_df)) 에서 제가 괄호안에 'data = group_df' 를 실수로 빼먹었었네요. 죄송합니다. ^^;

      summary(aov(y ~ group, data = group_df)) 에 data = group_df 추가해서 R script 수정해서 이제 분석결과 제대로 나왔습니다.

      잘못된 부분 지적해주셔서 감사합니다.

    • sfsfsf 2016.08.10 23:22  댓글주소  수정/삭제

      summary에서 aov table 부르실때 data= 옵션이 누락되어서 그렇게 나온 것 같네요!
      summary(aov(y ~ group, data=group_df))로 바꾸니 df가 제대로 나오네요... 왜 이렇게 된 건지. data= 생략했을 때 R이 어떻게 처리를 하는지는 저도 잘 모르겠습니다 ㅠ

    • sfsfsf 2016.08.10 23:24  댓글주소  수정/삭제

      앗ㅎㅎ 답글 달던 중이었는데 빠른 답장 감사합니다!! 오늘 처음 왔는데 자주 찾아뵙겠습니다 ㅎㅎ

    • R Friend R_Friend 2016.08.10 23:30 신고  댓글주소  수정/삭제

      ㅋㅋ, 채팅하는거 같네요.

      aov{stats} 에서 Arguments 찾아보니

      data : A data frame in which the variables specified in the formula will be found. If missing, the variables are searched for in the standard way.

      data 자리를 입력하지 않으면 standard way(?)로 변수를 찾는다고 나와있네요. y랑 group이 객체 values에 저장이 되어있어서 실행이 되긴 했나본데요, 결과가 왜 다르게 나왔는지는 저도 잘 모르겠습니다.

      자주 들러주시고, 댓글도 남겨주시면 저도 덕분에 배우고 좋겠습니다. ^^

  9. LDH 2016.09.01 11:34  댓글주소  수정/삭제  댓글쓰기

    항상 R 분석과 프로그래밍에 대해 많은 지식을 얻어 갑니다. 정말 감사합니다. 궁금한 점이 있는데 P 값이 유의수준보다 작으면 귀무가설을 채택하는 의미가 아닌가요??

  10. LDY 2016.09.06 22:56  댓글주소  수정/삭제  댓글쓰기

    꽤 오래 공부했음에도 정확하게 이해하지 못하고 느낌만 이해하다가, 이번에 완전히 이해했습니다. 양질의 자료 정말 감사드립니다!!

  11. cori 2017.04.23 13:43  댓글주소  수정/삭제  댓글쓰기

    안녕하세요,
    좋은 글 과제하는데 많은 참고가 되었습니다~^^
    근데, 한가지 질문이 있는데 유의수준 10%라는 부분은 어떻게 알 수 있는지 궁금합니다.
    기본적으로 R은 유의수준 5% 검정으로 알고 있는데,
    별도로 옵션을 설정한 곳이 없어서요.
    t.test에서는 conf.level로 설정이 가능한데, aov함수에서는 아무리 찾아도 안 보이네요.
    최종 결론 도출을 할 때 필요해서 여쭤봅니다.

  12. cori 2017.04.26 19:26  댓글주소  수정/삭제  댓글쓰기

    안녕하세요,
    빨리 댓글을 달아주셔서 감사합니다.
    근데, 제가 이제 막 R과 통계에 입문한 상태여서
    조금만 더 자세하게 설명해 주시면 감사하겠습니다.
    말씀하신 summary 부분이 아래 부분 같은데,
    이 곳에서 어떻게 10% 유의수준을 확인할 수 있는지와,
    별도로 유의수준을 5%로 변경 가능한지 시간 되실 때
    확인 부탁 드리겠습니다~ 죄송합니다~^^;;

    > 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


    출처: http://rfriend.tistory.com/131 [R, Python 분석과 프로그래밍 (by R Friend)]

    • R Friend R_Friend 2017.04.26 23:27 신고  댓글주소  수정/삭제

      안녕하세요 cori님,

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

      에서 순서대로 풀어보면요,
      (1) Df : 2
      (2) Sum Sq : 156.30
      (3) Mean Sq : 78.15 ( =(2)/(1) )
      (4) F value : 108.8
      (5) Pr(>F) : 1.2e-13***

      (4) 번의 F value가 F-test 값입니다.
      그리고 (5)번이 P-value 입니다. 1.2e-13 이면 0.00000000000012 이므로 유의 수준 5% 미만을 충족하고도 남을 만큼 매우 매우 작은 값임을 알 수 있습니다.

      뒤에 *** asterisk 는 아래 설명처럼 P-value 값에 따라서 '***', '**', '*', ' ' 로 가독성을 높여주기 위해서 뒤에 붙여주는 것입니다. *** 세개면 P-value가 0~0.001 로서 매우 유의하다고 판단할 수 있습니다.

      Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


      ========================
      Bartlett test는 아래 결과 처럼 p-value = 0.4368 이라고 명확하게 표기되어 있으므로 혼란이 없을 것으로 생각합니다.
      -------------
      Bartlett test of homogeneity of variances

      data: y by group
      Bartlett's K-squared = 1.6565, df = 2, p-value = 0.4368 <= 바로 이 값이 p-value

    • aov 2019.04.06 16:25  댓글주소  수정/삭제

      group_df <- transform(group_df, group = factor(group))

      factor 처리 안하면 그렇게 나오네요

  13. zobiet 2017.06.22 17:48  댓글주소  수정/삭제  댓글쓰기

    정말 좋은 자료 올려주셔서 감사합니다^^

  14. cymen 2017.07.23 18:11  댓글주소  수정/삭제  댓글쓰기

    글 너무 잘 읽었습니다 예제가 잘되있어서 따라하기 너무 좋아요 !!
    예제를 따라하다가 궁금증이 생겼는데요
    예제를 따라하고 나온 결과입니다.
    Df Sum Sq Mean Sq F value Pr(>F)
    group 1 130.56 130.56 81 9.34e-10 ***
    Residuals 28 45.13 1.61


    선생님의 예제는
    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


    "검정통계량 F-value가 81으로 나왔으며, P-value 값은 '1.2e-13'으로서 매우 작게 나와..."

    여기서 궁금증이 생겼는데 선생님의 예제에는 F-value가 108.8 이 나온것으로 보이는데 F-value 가 제가 실행하고 나온 예제 결과는 81이 나왔는데 어떤게 맞는건가요?? 아니면 무슨 차이가 숨어
    있는건지요?? 이부분만 다른결과가 나와서 궁금하네요 ㅠㅠ

    • R Friend R_Friend 2017.07.23 22:20 신고  댓글주소  수정/삭제

      안녕하세요 cymen 님,
      반갑습니다.

      본문 R script 다시 보면서 실행해보니 F통계량은 108.8이 맞네요. 본문 내용 수정하였습니다.

      -----------
      아래 처럼 R script 쓰는게 맞구요, 이러면 F-value가 108.8 이 나옵니다.
      summary(aov(y ~ group, data = group_df)) <== 올바른 R script
      -----------

      -----------
      반면에, 만약 아래처럼 괄호안에 data = group_df 를 안쓰면 F-value로 81이 나오게 됩니다. 이게 예전에 제가 실수로 잘못 썼던 R script 인데요, 이렇게 하면 안된답니다. ㅜ_ㅜ

      summary(aov(y ~ group)
      <== 잘못된 R script (괄호 안에 data = group_df 명기 필요)
      -----------

      도움 되었기를 바랍니다.

  15. Lawine 2019.03.06 18:20  댓글주소  수정/삭제  댓글쓰기

    안녕하세요 도움 많이 받고 있습니다.

    수식을 보다 이상한 점이 보여서 덧글 답니다.

    글로 표현하려니 이상하긴 한데

    [분산분석 기본 원리 2]의 ② 처리제곱합(SSTR)에서 Y_i가 아니라 Y ̅_i가 아닌가 싶네요

    확인 부탁드립니다.

    • R Friend R_Friend 2019.04.10 21:42 신고  댓글주소  수정/삭제

      안녕하세요 Lawine님,

      댓글에 남겨주신 말씀처럼 제가 실수를 했네요. 댓글 남겨주신 덕분에 포스팅 본문의 수식 수정하였습니다.

      수식 꼼꼼히 봐주시고 댓글 남겨주셔서 고맙습니다.

  16. Sanochi1031 2019.11.29 15:38  댓글주소  수정/삭제  댓글쓰기

    경제학 공부하면서 R을 처음 써보는데, 선생님 홈페이지 보면서, 상당히 편리하게 이해했습니다. 좋은 글과 정보에 정말 감사드립니다.