R에서 양방향 anova를 수행하는 방법

양방향 ANOVA (“분산 분석”)는 두 요인에 걸쳐 분할된 3개 이상의 독립 그룹의 평균 간에 통계적으로 유의미한 차이가 있는지 여부를 확인하는 데 사용됩니다.

이 튜토리얼에서는 R에서 양방향 ANOVA를 수행하는 방법을 설명합니다.

예: R의 양방향 ANOVA

운동 강도와 성별이 체중 감량에 영향을 미치는지 확인하고 싶다고 가정해 보겠습니다. 이 경우 우리가 보고 있는 두 가지 요소는 운동 과 성별 이며, 응답 변수는 파운드 단위로 측정되는 체중 감소 입니다.

운동과 성별이 체중 감량에 영향을 미치는지 확인하고 운동과 성별이 체중 감량에 상호 작용이 있는지 확인하기 위해 양방향 분산 분석을 수행할 수 있습니다.

우리는 남성 30명, 여성 30명을 무작위로 10명씩 무작위로 할당하여 한 달 동안 운동을 하지 않거나 가벼운 운동 또는 격렬한 운동 프로그램을 따르는 실험에 참여하도록 모집하고 있습니다.

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

 #make this example reproducible
set.seed(10)

#create data frame
data <- data.frame(gender = rep(c("Male", "Female"), each = 30),
                   exercise = rep(c("None", "Light", "Intense"), each = 10, times = 2),
                   weight_loss = c(runif(10, -3, 3), runif(10, 0, 5), runif(10, 5, 9),
                                   runif(10, -4, 2), runif(10, 0, 3), runif(10, 3, 8)))

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

# gender exercise weight_loss
#1 Male None 0.04486922
#2 Male None -1.15938896
#3 Male None -0.43855400
#4 Male None 1.15861249
#5 Male None -2.48918419
#6 Male None -1.64738030

#see how many participants are in each group
table(data$gender, data$exercise)

# Intense Light None
# Female 10 10 10
# Male 10 10 10

데이터 탐색

양방향 ANOVA 모델을 적용하기 전에 dplyr 패키지를 사용하여 6개 치료 그룹 각각에 대한 체중 감소의 평균 및 표준 편차를 찾아 데이터를 더 잘 이해할 수 있습니다.

 #load dplyr package
library(dplyr)

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

# A tibble: 6 x 4
# Groups: gender [2]
# gender exercise means sd
#          
#1 Female Intense 5.31 1.02 
#2 Female Light 0.920 0.835
#3 Female None -0.501 1.77 
#4 Male Intense 7.37 0.928
#5 Male Light 2.13 1.22 
#6 Male None -0.698 1.12 

또한 6개 치료 그룹 각각에 대한 상자 그림을 만들어 각 그룹의 체중 감소 분포를 시각화할 수도 있습니다.

 #set margins so that axis labels on boxplot don't get cut off
by(mar=c(8, 4.1, 4.1, 2.1))

#create boxplots
boxplot(weight_loss ~ gender:exercise,
data = data,
main = "Weight Loss Distribution by Group",
xlab = "Group",
ylab = "Weight Loss",
col = "steelblue",
border = "black", 
las = 2 #make x-axis labels perpendicular
)

강렬한 운동에 참여한 두 그룹의 체중 감량 값이 더 높은 것으로 즉시 알 수 있습니다. 또한 강렬한 운동 그룹과 가벼운 운동 그룹 모두에서 남성이 여성보다 체중 감량 값이 더 높은 경향이 있음을 알 수 있습니다.

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

양방향 ANOVA 모델 피팅

R에서 양방향 ANOVA 모델을 피팅하는 일반적인 구문은 다음과 같습니다.

aov(응답 변수 ~predictor_variable1 *predictor_variable2, 데이터 = 데이터 세트)

두 예측 변수 사이의 * 는 두 예측 변수 간의 상호 작용 효과도 테스트하고 싶다는 것을 나타냅니다.

이 예에서는 체중 감량을 응답 변수로 사용하고 성별 과 운동을 두 예측 변수로 사용하여 양방향 ANOVA 모델을 맞추기 위해 다음 코드를 사용할 수 있습니다.

그런 다음 summary() 함수를 사용하여 모델 결과를 표시할 수 있습니다.

 #fit the two-way ANOVA model
model <- aov(weight_loss ~ gender * exercise, data = data)

#view the model output
summary(model)

# Df Sum Sq Mean Sq F value Pr(>F)    
#gender 1 15.8 15.80 11.197 0.0015 ** 
#exercise 2 505.6 252.78 179.087 <2e-16 ***
#gender:exercise 2 13.0 6.51 4.615 0.0141 *  
#Residuals 54 76.2 1.41                   
#---
#Significant. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

모델 결과를 보면 성별 , 운동 , 두 변수 간의 상호작용이 모두 유의수준 0.05에서 통계적으로 유의미한 것으로 나타났습니다.

모델 가정 확인

더 나아가기 전에 모델 결과가 신뢰할 수 있도록 모델의 가정이 충족되는지 확인해야 합니다. 특히 양방향 ANOVA는 다음을 가정합니다.

1. 독립성 – 각 그룹의 관찰은 서로 독립적이어야 합니다. 무작위 설계를 사용했기 때문에 이 가정이 충족되어야 하므로 크게 걱정할 필요가 없습니다.

2. 정규성 – 종속 변수는 두 요인 그룹의 각 조합에 대해 대략적인 정규 분포를 가져야 합니다.

이 가정을 테스트하는 한 가지 방법은 모델 잔차의 히스토그램을 만드는 것입니다. 잔차가 대략적으로 정규 분포를 따르는 경우 이 가정이 충족되어야 합니다.

 #define model residuals
reside <- model$residuals

#create histogram of residuals
hist(resid, main = "Histogram of Residuals", xlab = "Residuals", col = "steelblue")

잔차는 대략 정규 분포를 따르므로 정규성 가정이 충족된다고 가정할 수 있습니다.

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

이 가정을 확인하는 한 가지 방법은 car 패키지를 사용하여 등분산 테스트를 수행하는 것입니다.

 #load car package
library(car)

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

#Levene's Test for Homogeneity of Variance (center = median)
# Df F value Pr(>F)
#group 5 1.8547 0.1177
#54  

검정의 p-값이 유의 수준 0.05보다 크므로 그룹 간 분산 동일성 가정이 충족된다고 가정할 수 있습니다.

치료 차이점 분석

모델 가정이 충족되는지 확인한 후에는 사후 테스트를 수행하여 정확히 어떤 치료 그룹이 서로 다른지 확인할 수 있습니다.

사후 테스트에서는 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 ~ gender * exercise, data = data)
#
#$gender
# diff lwr upr p adj
#Male-Female 1.026456 0.4114451 1.641467 0.0014967
#
#$exercise
# diff lwr upr p adj
#Light-Intense -4.813064 -5.718493 -3.907635 0.0e+00
#None-Intense -6.938966 -7.844395 -6.033537 0.0e+00
#None-Light -2.125902 -3.031331 -1.220473 1.8e-06
#
#$`gender:exercise`
# diff lwr upr p adj
#Male:Intense-Female:Intense 2.0628297 0.4930588 3.63260067 0.0036746
#Female:Light-Female:Intense -4.3883563 -5.9581272 -2.81858535 0.0000000
#Male:Light-Female:Intense -3.1749419 -4.7447128 -1.60517092 0.0000027
#Female:None-Female:Intense -5.8091131 -7.3788841 -4.23934219 0.0000000
#Male:None-Female:Intense -6.0059891 -7.5757600 -4.43621813 0.0000000
#Female:Light-Male:Intense -6.4511860 -8.0209570 -4.88141508 0.0000000
#Male:Light-Male:Intense -5.2377716 -6.8075425 -3.66800066 0.0000000
#Female:None-Male:Intense -7.8719429 -9.4417138 -6.30217192 0.0000000
#Male:None-Male:Intense -8.0688188 -9.6385897 -6.49904786 0.0000000
#Male:Light-Female:Light 1.2134144 -0.3563565 2.78318536 0.2185439
#Female:None-Female:Light -1.4207568 -2.9905278 0.14901410 0.0974193
#Male:None-Female:Light -1.6176328 -3.1874037 -0.04786184 0.0398106
#Female:None-Male:Light -2.6341713 -4.2039422 -1.06440032 0.0001050
#Male:None-Male:Light -2.8310472 -4.4008181 -1.26127627 0.0000284
#Male:None-Female:None -0.1968759 -1.7666469 1.37289500 0.9990364

p-값은 각 그룹 간에 통계적으로 유의한 차이가 있는지 여부를 나타냅니다.

예를 들어, 위의 마지막 행에서 운동을 하지 않은 남성 그룹은 운동을 하지 않은 여성 그룹에 비해 체중 감량에 통계적으로 유의미한 차이가 없음을 알 수 있습니다(p-값: 0.990364).

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

 #set axis margins so labels don't get cut off
by(mar=c(4.1, 13, 4.1, 2.1))

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

양방향 ANOVA 결과 보고

마지막으로 결과를 요약하는 방식으로 양방향 ANOVA의 결과를 보고할 수 있습니다.

성별( 남성, 여성) 및 운동 프로그램 (없음, 가벼운 운동, 강렬한 운동)이 체중 감량 (파운드로 측정)에 미치는 영향을 조사하기 위해 양방향 ANOVA를 수행했습니다. 성별과 운동이 체중 감량에 미치는 영향 사이에는 통계적으로 유의미한 상호작용이 있었습니다(F(2, 54) = 4.615, p = 0.0141). 사후 Tukey의 HSD 테스트가 수행되었습니다.

남성의 경우, 강렬한 운동 프로그램은 가벼운 프로그램(p < 0.0001)이나 운동 프로그램을 하지 않는 것 (p < 0.0001)에 비해 체중 감량 효과가 훨씬 컸습니다. 또한 남성의 경우 가벼운 식사는 운동을 하지 않는 것 보다 훨씬 더 큰 체중 감소를 가져왔습니다(p < 0.0001).

여성의 경우, 강렬한 운동 프로그램은 가벼운 프로그램(p < 0.0001)이나 운동 프로그램을 하지 않는 것 (p < 0.0001)에 비해 훨씬 더 큰 체중 감소를 가져왔습니다.

ANOVA의 가정이 충족되는지 확인하기 위해 정규성 검사와 Levene 테스트를 수행했습니다.