일원 분산분석 1 2

데이터와 탐색적 데이터 분석

타이어 데이터셋Bidyut Ghosh(2017), “One-way ANOVA in R” 일원배치 분산분석(One-way Anova) 블로그 포스트에 제공되고 있다. 인도에서 수집된 데이터로 타이어 브랜드별(Brands)로 연비(Mileage) 차이가 있는가라는 관심을 두고 있다.

물론 귀무가설 첫번째 가정은 타이어 브랜드별로 연비가 평균적으로 동일하다는 것이다.

\[\mu_{Apollo} = \mu_{Bridgestone} = \mu_{CEAT} = \mu_{Falken}\]

데이터를 바로 다운로드 받아 자료형에 맞춰 변환을 한후 skim 팩키지로 탐색적으로 기술 통계량을 뽑아보고, 상자그림으로 시각화한다.

# 0. 환경설정 -----
library(tidyverse)
library(skimr)
library(broom)
library(cowplot)

# 1. 데이터 -----
download.file(url="https://datascienceplus.com/wp-content/uploads/2017/08/tyre.csv", destfile = "data/tyre.csv")
tire_df <- read_csv("data/tyre.csv")

tire_df <- tire_df %>% 
  mutate(Brands = fct_infreq(Brands))

# 2. 탐색적 데이터 분석 -----
## skim
tire_df %>% 
  group_by(Brands) %>% 
  skim()
Skim summary statistics
 n obs: 60 
 n variables: 2 
 group variables: Brands 

-- Variable type:numeric ---------------------------------------------------------------------------------------------------------------------
      Brands variable missing complete  n  mean   sd    p0   p25   p50
      Apollo  Mileage       0       15 15 34.8  2.22 30.62 32.89 34.84
 Bridgestone  Mileage       0       15 15 31.78 2.2  27.88 30.9  32   
        CEAT  Mileage       0       15 15 34.76 2.53 30.43 33.11 34.78
      Falken  Mileage       0       15 15 37.62 1.7  34.31 36.66 37.38
   p75  p100     hist
 36.42 38.33 <U+2582><U+2582><U+2587><U+2582><U+2585><U+2585><U+2587><U+2585>
 33.3  35.01 <U+2583><U+2582><U+2581><U+2587><U+2582><U+2586><U+2583><U+2583>
 36.13 41.05 <U+2582><U+2586><U+2583><U+2587><U+2587><U+2581><U+2581><U+2582>
 38.47 40.66 <U+2582><U+2581><U+2586><U+2587><U+2586><U+2582><U+2582><U+2583>
## 상자그림
tire_df %>% 
  ggplot(aes(x=Brands, y=Mileage, fill=Brands)) +
    geom_boxplot() +
    coord_flip() +
    geom_jitter(colour = "dark gray", width=.1) +
    stat_boxplot(geom ='errorbar',width = 0.4) +
    stat_summary(geom="point", fun.y=mean, color="blue") +
    guides(fill = FALSE) +
    theme_minimal()

분산분석

일원배치 분산분석을 aov() 함수로 수행하면 연비는 브랜드별로 차이가 있다고 나타난다. 하지만, 문제는 제조사가 4개로 조합의 가지수가 \(\choose_{3}{2} = 6\) 만큼 나오게 되는데, 어디서 연비차이가 나는지 확인하기 위해서 TukeyHSD() 검정을 도입하고 이를 시각적으로 표현한다.

# 3. 분산분석 -----
## 전체 비교
tire_aov <- aov(Mileage ~ Brands, data=tire_df)
tidy(tire_aov)
       term df    sumsq    meansq statistic      p.value
1    Brands  3 256.2908 85.430262  17.94151 2.780989e-08
2 Residuals 56 266.6494  4.761597        NA           NA
## 쌍체 비교
TukeyHSD(tire_aov, conf.level = 0.95)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Mileage ~ Brands, data = tire_df)

$Brands
                          diff        lwr       upr     p adj
Bridgestone-Apollo -3.01900000 -5.1288190 -0.909181 0.0020527
CEAT-Apollo        -0.03792661 -2.1477456  2.071892 0.9999608
Falken-Apollo       2.82553333  0.7157143  4.935352 0.0043198
CEAT-Bridgestone    2.98107339  0.8712544  5.090892 0.0023806
Falken-Bridgestone  5.84453333  3.7347143  7.954352 0.0000000
Falken-CEAT         2.86345994  0.7536409  4.973279 0.0037424
TukeyHSD(tire_aov, conf.level = 0.95) %>% 
  tidy() %>% 
  ggplot(aes(x=fct_reorder(comparison, -estimate), y=estimate, color=comparison)) +
    geom_point() +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high)) +
    geom_hline(yintercept = 0, color="darkgray", linetype="dashed") +
    coord_flip() +
    labs(x="", y="브랜드별 차이", title="신뢰구간 95% 쌍체 비교") +
    theme_minimal() +
    guides(color=FALSE)

브랜드 비교 6개 조합수 중에서 CEAT-Apollo 두 브랜드를 제외하고 모두 연비차이가 나는 것으로 파악된다.

모형 가정 검정

일원배치 분산분석모형은 독립성(independence), 등분산성(homogeneity of variance, homoscedasticity), 그리고 정규성(normality)이 충족되어야 모형에서 나온 결과가 유의미하다.

독립성은 브랜드가 다르기 때문에 서로 다른 회사에서 제조되어 특이한 경우 동일공장에서 브랜드를 달리해서 타이어를 생산하거나, 지분관계로 연관성이 없다고 본다면 충족되는 것으로 일단 판정한다. 그리고 나서, 정규성 검정을 해야하는데 정규성은 Shapiro-Wilk, Kolmogorov-Smirnov 검정등이 있지만, qqplot으로 갈음한다. 마지막으로 등분산성은 bartlett.test() 검정을 통해 등분산성이 충족되는 것이 확인된다.

## 가정 검정
tire_glm <- glm(Mileage ~ Brands, data=tire_df)

glm_qqplot_g <- augment(tire_glm) %>% 
  ggplot(aes(sample = .resid)) +
    stat_qq() +
    stat_qq_line(distribution = stats::qnorm, color="lightblue",size=1) +
    theme_minimal() +
    labs(x="이론값", y="관측표본", title="GLM 모형") 

aov_qqplot_g <- residuals(tire_aov) %>% 
  as_tibble %>% 
  ggplot(aes(sample = value)) +
  stat_qq() +
  stat_qq_line(distribution = stats::qnorm, color="lightblue",size=1) +
  theme_minimal() +
  labs(x="이론값", y="관측표본", title="분산분석 모형") 

plot_grid(glm_qqplot_g, aov_qqplot_g)

### 등분산
bartlett.test(Mileage~Brands, tire_df)

    Bartlett test of homogeneity of variances

data:  Mileage by Brands
Bartlett's K-squared = 2.1496, df = 3, p-value = 0.5419