R에서 일원 분산 분석을 수행하는 방법

일원 분산 분석은 3개 이상의 독립 그룹 평균 간에 통계적으로 유의한 차이가 있는지 여부를 확인하는 데 사용됩니다.

이러한 유형의 테스트는 예측 변수가 반응 변수에 미치는 영향을 분석하기 때문에 일원 분산 분석이라고 합니다.

참고 : 반응 변수에 대한 두 예측 변수의 영향에 관심이 있는 경우 양방향 분산 분석을 수행할 수 있습니다.

R에서 일원 분산 분석을 수행하는 방법

다음 예에서는 R에서 일원 분산 분석을 수행하는 방법을 보여줍니다.

배경

세 가지 다른 운동 프로그램이 체중 감량에 서로 다른 영향을 미치는지 확인하고 싶다고 가정해 보겠습니다. 우리가 연구하는 예측 변수는 운동 프로그램 이고 반응 변수는 파운드 단위로 측정되는 체중 감소 입니다.

일원 분산 분석을 수행하여 세 가지 프로그램으로 인한 체중 감량 간에 통계적으로 유의미한 차이가 있는지 확인할 수 있습니다.

우리는 한 달 동안 프로그램 A, 프로그램 B 또는 프로그램 C를 따르도록 무작위로 30명을 할당하는 실험에 참여할 90명을 모집합니다.

다음 코드는 우리가 작업할 데이터 프레임을 생성합니다.

 #make this example reproducible
set.seed(0)

#create data frame
data <- data.frame(program = rep(c("A", "B", "C"), each = 30),
                   weight_loss = c(runif(30, 0, 3),
                                   runif(30, 0, 5),
                                   runif(30, 1, 7)))

#view first six rows of data frame
head(data)

# program weight_loss
#1 A 2.6900916
#2 A 0.7965260
#3 A 1.1163717
#4 A 1.7185601
#5 A 2.7246234
#6 A 0.6050458

데이터 프레임의 첫 번째 열은 그 사람이 한 달 동안 참여한 프로그램을 보여주고 두 번째 열은 프로그램이 끝날 때 그 사람이 경험한 총 체중 감소를 파운드 단위로 보여줍니다.

데이터 탐색

일원 분산 분석 모델을 적용하기 전에 dplyr 패키지를 사용하여 세 가지 프로그램 각각에 대한 체중 감량의 평균 및 표준 편차를 찾아 데이터를 더 잘 이해할 수 있습니다.

 #load dplyr package
library (dplyr)

#find mean and standard deviation of weight loss for each treatment group
data %>%
  group_by (program) %>%
  summarize (mean = mean(weight_loss),
            sd = sd(weight_loss))

# A tibble: 3 x 3
# program mean sd
#      
#1 A 1.58 0.905
#2 B 2.56 1.24 
#3 C 4.13 1.57  

또한 세 가지 프로그램 각각에 대한 상자 그림을 만들어 각 프로그램의 체중 감량 분포를 시각화할 수도 있습니다.

 #create boxplots
boxplot(weight_loss ~ program,
data = data,
main = "Weight Loss Distribution by Program",
xlab = "Program",
ylab = "Weight Loss",
col = "steelblue",
border = "black")

이 상자 그림에서 프로그램 C 참가자의 평균 체중 감소가 가장 높고 프로그램 A 참가자의 평균 체중 감소가 가장 낮다는 것을 알 수 있습니다.

또한 체중 감량에 대한 표준 편차(상자 그림의 “길이”)가 다른 두 프로그램에 비해 프로그램 C에서 약간 더 높다는 것을 알 수 있습니다.

다음으로, 이러한 시각적 차이가 실제로 통계적으로 유의한지 확인하기 위해 일원 분산 분석 모델을 데이터에 적용하겠습니다.

일원 분산 분석 모델 피팅

R에서 일원 분산 분석 모델을 피팅하는 일반적인 구문은 다음과 같습니다.

aov(응답 변수 ~ 예측_변수, 데이터 = 데이터 세트)

이 예에서는 다음 코드를 사용하여 체중 손실을 응답 변수로 사용하고 프로그램을 예측 변수로 사용하여 일원 분산 분석 모델을 맞출 수 있습니다. 그런 다음 summary() 함수를 사용하여 모델 결과를 표시할 수 있습니다.

 #fit the one-way ANOVA model
model <- aov(weight_loss ~ program, data = data)

#view the model output
summary(model)

# Df Sum Sq Mean Sq F value Pr(>F)    
#program 2 98.93 49.46 30.83 7.55e-11 ***
#Residuals 87 139.57 1.60                     
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

모델 결과에서 예측변수 프로그램이 0.05 유의수준에서 통계적으로 유의미하다는 것을 알 수 있습니다.

즉, 세 가지 프로그램의 평균 체중 감량 간에는 통계적으로 유의미한 차이가 있는 것입니다.

모델 가정 확인

더 나아가기 전에 모델 결과가 신뢰할 수 있도록 모델의 가정이 충족되는지 확인해야 합니다. 특히 일원 분산 분석에서는 다음을 가정합니다.

1. 독립성 – 각 그룹의 관찰은 서로 독립적이어야 합니다. 무작위 설계(즉, 참가자를 운동 프로그램에 무작위로 할당)를 사용했기 때문에 이 가정이 충족되므로 너무 걱정할 필요가 없습니다.

2. 정규성 – 종속변수는 예측변수의 각 수준에 대해 대략 정규분포를 가져야 합니다.

3. 등분산 – 각 그룹의 분산이 동일하거나 거의 동일합니다.

정규성 과 등분산 의 가정을 확인하는 한 가지 방법은 4개의 모델 확인 플롯을 생성하는 plot() 함수를 사용하는 것입니다. 특히 우리는 다음 두 가지 플롯에 특히 관심이 있습니다.

• 잔차 대 적합 – 이 그래프는 잔차와 적합치 사이의 관계를 보여줍니다. 이 그래프를 사용하여 그룹 간의 분산이 대략 같은지 여부를 대략적으로 평가할 수 있습니다.
• QQ 플롯 – 이 플롯은 이론적 분위수에 대한 표준화된 잔차를 표시합니다. 이 그래프를 사용하여 정규성 가정이 충족되는지 여부를 대략적으로 평가할 수 있습니다.

다음 코드를 사용하여 이러한 모델 검사 플롯을 생성할 수 있습니다.

 plot(model)

위의 QQ 그래프를 통해 정규성 가정을 확인할 수 있습니다. 이상적으로는 표준화된 잔차가 그림의 직선 대각선을 따라 놓이게 됩니다. 그러나 위의 그래프에서는 잔차가 선에서 시작과 끝으로 약간 벗어나는 것을 볼 수 있습니다. 이는 정규성 가정이 위반될 수 있음을 나타냅니다.

잔차 대 위 의 조정된 그래프를 통해 등분산 가정을 확인할 수 있습니다. 이상적으로는 잔차가 적합치의 각 수준에 균등하게 분포되기를 원합니다.

적합치가 높을수록 잔차가 훨씬 더 많이 퍼져 있는 것을 볼 수 있는데, 이는 등분산 가정이 위반될 수 있음을 나타냅니다.

등분산을 공식적으로 테스트하기 위해 car 패키지를 사용하여 Levene 테스트를 수행할 수 있습니다.

 #load car package
library (car)

#conduct Levene's Test for equality of variances
leveneTest(weight_loss ~ program, data = data)

#Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)  
#group 2 4.1716 0.01862 *
#87                  
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

검정의 p-값은 0.01862 입니다. 유의 수준 0.05를 사용하면 세 프로그램에서 분산이 동일하다는 귀무 가설을 기각합니다. 그러나 유의수준 0.01을 사용하면 귀무가설을 기각하지 않습니다.

정규성과 등분산에 대한 가정이 충족되도록 데이터를 변환하려고 시도할 수 있지만 지금은 이에 대해 너무 걱정하지 않겠습니다.

치료 차이점 분석

모델 가정이 충족되는지(또는 합리적으로 충족되는지) 확인되면 사후 테스트를 수행하여 정확히 어떤 치료 그룹이 서로 다른지 확인할 수 있습니다.

사후 테스트에서는 TukeyHSD() 함수를 사용하여 다중 비교를 위한 Tukey 테스트를 수행합니다.

 #perform Tukey's Test for multiple comparisons
TukeyHSD(model, conf.level=.95) 

#Tukey multiple comparisons of means
# 95% family-wise confidence level
#
#Fit: aov(formula = weight_loss ~ program, data = data)
#
#$program
# diff lwr upr p adj
#BA 0.9777414 0.1979466 1.757536 0.0100545
#CA 2.5454024 1.7656076 3.325197 0.0000000
#CB 1.5676610 0.7878662 2.347456 0.0000199

p-값은 각 프로그램 간에 통계적으로 유의미한 차이가 있는지 여부를 나타냅니다. 그 결과, 각 프로그램의 평균 체중 감량은 유의수준 0.05에서 통계적으로 유의미한 차이가 있는 것으로 나타났습니다.

또한 R의 plot(TukeyHSD()) 함수를 사용하여 Tukey 테스트의 결과인 95% 신뢰 구간을 시각화할 수 있습니다.

 #create confidence interval for each comparison
plot(TukeyHSD(model, conf.level=.95), las = 2)

신뢰 구간의 결과는 가설 검정 결과와 일치합니다.

특히, 프로그램 간의 평균 체중 감량에 대한 신뢰구간 중 어느 것도 0 값을 포함하고 있지 않음을 알 수 있는데, 이는 세 프로그램 간의 평균 체중 감량에 통계적으로 유의미한 차이가 있음을 나타냅니다.

이는 가설 검정 에 대한 모든 p-값이 0.05 미만인 것과 일치합니다.

일원 분산 분석 결과 보고

마지막으로 결과를 요약하는 방식으로 일원 분산 분석의 결과를 보고할 수 있습니다.

운동 프로그램의 효과를 조사하기 위해 일원 분산 분석을 수행했습니다.   체중 감량 (파운드로 측정). 체중 감량에 대한 세 가지 프로그램의 효과 간에는 통계적으로 유의미한 차이가 있었습니다(F(2, 87) = 30.83, p = 7.55e-11). 사후 Tukey의 HSD 테스트가 수행되었습니다.

프로그램 C 참가자의 평균 체중 감소는 프로그램 B 참가자의 평균 체중 감소보다 훨씬 큽니다(p < 0.0001).

프로그램 C 참가자의 평균 체중 감소는 프로그램 A 참가자의 평균 체중 감소보다 훨씬 큽니다(p < 0.0001).

또한 프로그램 B 참가자의 평균 체중 감소는 프로그램 A 참가자의 평균 체중 감소보다 훨씬 컸습니다(p = 0.01).

추가 리소스

다음 자습서에서는 일원 분산 분석에 대한 추가 정보를 제공합니다.

일원 분산 분석 소개
ANOVA를 통한 사후 테스트 사용 가이드
전체 가이드: ANOVA 결과 보고 방법