[R] 13. 분산분석(ANOVA)

예시로 쓰일 데이터 예제

set.seed(2021)

# 임의로 데이터를 생성한다. (100명의 유저가 특정 곡을 스트리밍한 이력)
temp <- tibble(
  user_id = c(10000:10149),
  user_age = sample(x = c("20대", "30대"), size = 150, replace = TRUE),
  user_gender = sample(x = c("남성", "여성"), size = 150, prob = c(0.5, 0.5), replace = TRUE),
  song_id = sample(x = letters[1:5], size = 150, replace = TRUE),
  streaming_count = rpois(n = 150, lambda = 20),
  download_count = rpois(n = 150, lambda = 5)
)

temp

## # A tibble: 150 x 6
##    user_id user_age user_gender song_id streaming_count download_count
##      <int> <chr>    <chr>       <chr>             <int>          <int>
##  1   10000 20대     여성        c                    15              3
##  2   10001 30대     여성        a                    23              4
##  3   10002 30대     여성        a                    14              5
##  4   10003 30대     남성        d                    24              2
##  5   10004 20대     여성        b                    17              3
##  6   10005 30대     남성        c                    20              1
##  7   10006 30대     남성        b                    23              5
##  8   10007 30대     여성        d                    24              6
##  9   10008 30대     남성        c                    21              8
## 10   10009 30대     남성        d                    18              6
## # … with 140 more rows

분산분석(ANOVA, analysis of variance)

• 범주형 변수의 수준이 세 수준 이상일 때 연속형 변수인 종속변수의 평균을 비교하는 통계적 기법입니다.
  • t-Test와 차이가 있다면 범주형 변수의 수준이 이표본이 아닌 세 가지 이상일 때 쓰일 수 있습니다.
• 분산분석에 공변량(covariate) 변수를 추가로 통제할 경우 이를 공분산분석 이라고 부릅니다.
  • 공변량: 여러 변수들이 공통적으로 함께 공유하고 있는 변량
  • 독립변수 이외에 종속변수에 영향을 줄 수 있는 노이즈를 통제하고자 하는 변수라고 이해하시면 됩니다.
• 분산분석은 aov() 함수를 이용하며 기존에 보았던 것과 같이 종속변수 ~ 독립변수와 같은 formula 포맷을 파라미터로 받습니다.
• 그리고 이 결과는 summary() 함수를 통해 확인할 수 있습니다.

 

• 예시 데이터에서 곡(song_id) 별 스트리밍 횟수의 평균 차이를 분산분석을 통해 살펴보겠습니다.

#aov() 함수를 적용한 후 결과를 불러오기 위해 summary() 함수 사용
temp %>% 
  aov(streaming_count ~ song_id, data = .) %>% 
  summary()

##              Df Sum Sq Mean Sq F value Pr(>F)
## song_id       4    171   42.79    1.72  0.149
## Residuals   145   3608   24.88

• 마찬가지로 broom 라이브러리의 tidy() 함수를 이용하여 tibble 포맷의 결과를 출력해볼 수 있습니다.

#aov() 함수를 적용한 후 tidy() 함수 사용
temp %>% 
  aov(streaming_count ~ song_id, data = .) %>% 
  tidy(.)

## # A tibble: 2 x 6
##   term         df sumsq meansq statistic p.value
##   <chr>     <dbl> <dbl>  <dbl>     <dbl>   <dbl>
## 1 song_id       4  171.   42.8      1.72   0.149
## 2 Residuals   145 3608.   24.9     NA     NA

• 결과를 해석하자면 p-value 값이 유의수준 0.05보다 크므로 유의수준 0.05 하에 귀무가설을 기각할 수 없습니다.
  • 다섯 개의 곡의 평균 스트리밍 횟수 차이는 없다고 볼 수 있습니다.
• 하지만 위의 결과는 햑심적인 정보를 제공해 주지는 못합니다.
  왜냐하면 어떤 변수들의 평균값이 서로 다른 것인지 구체적으로 확인할 수 없기 때문입니다.
• 각 수준별 평균차이를 살펴보기 위해서는 흔히 사후비교를 실시합니다.
• 많이 사용되는 방법으로 Tukey’s HSD(Honest significance difference) 기법을 적용해보겠습니다.
  • 함수는 TukeyHSD() 함수를 활용합니다.

# aov() 적용 이후 TukeyHSD() 함수 적용
temp %>% 
  aov(streaming_count ~ song_id, data = .) %>% 
  TukeyHSD(., which = "song_id") %>% 
  #tidy() %>% filter(adj.p.value < 0.05)
  tidy()

## # A tibble: 10 x 7
##    term    contrast null.value estimate conf.low conf.high adj.p.value
##    <chr>   <chr>         <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
##  1 song_id b-a               0    1.16     -2.27     4.59        0.884
##  2 song_id c-a               0   -1.96     -5.61     1.69        0.576
##  3 song_id d-a               0   -1.38     -4.65     1.89        0.773
##  4 song_id e-a               0   -0.743    -4.35     2.87        0.979
##  5 song_id c-b               0   -3.12     -6.89     0.657       0.157
##  6 song_id d-b               0   -2.53     -5.94     0.873       0.246
##  7 song_id e-b               0   -1.90     -5.63     1.83        0.624
##  8 song_id d-c               0    0.583    -3.05     4.21        0.992
##  9 song_id e-c               0    1.22     -2.72     5.15        0.913
## 10 song_id e-d               0    0.633    -2.95     4.22        0.988

• 위 결과를 보면 각 변인간 (contrast 컬럼) p-value 값(adj.p.value 컬럼)을 확인해볼 수 있습니다.
• 신뢰구간 conf.low, conf.high 사이에 0이 속하는 것으로 보아 각 변인들간 평균차이는 0이라고 볼 수 있습니다.