Как выполнить двусторонний дисперсионный анализ в r

Двусторонний дисперсионный анализ («дисперсионный анализ») используется для определения того, существует ли статистически значимая разница между средними значениями трех или более независимых групп, которые были разделены по двум факторам.

В этом руководстве объясняется, как выполнить двусторонний дисперсионный анализ в R.

Пример: двусторонний дисперсионный анализ в R

Допустим, мы хотим определить, влияют ли интенсивность тренировок и пол на потерю веса. В данном случае мы рассматриваем два фактора — физические упражнения и пол , а переменной ответа — потерю веса, измеряемую в фунтах.

Мы можем выполнить двусторонний дисперсионный анализ, чтобы определить, влияют ли физические упражнения и пол на потерю веса, а также определить, существует ли взаимодействие между физическими упражнениями и полом на потерю веса.

Мы набираем 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 :

 #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 

Мы также можем создать коробчатую диаграмму для каждой из шести групп лечения, чтобы визуализировать распределение потери веса для каждой группы:

 #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

Общий синтаксис для подбора двуфакторной модели ANOVA в R:

aov(переменная ответа ~predictor_variable1 *predictor_variable2, data = набор данных)

Обратите внимание, что знак * между двумя переменными-предикторами указывает на то, что мы также хотим проверить эффект взаимодействия между двумя переменными-предикторами.

В нашем примере мы можем использовать следующий код для соответствия модели двустороннего дисперсионного анализа, используя вес_потерю в качестве переменной ответа, а пол и упражнения в качестве двух переменных-предсказателей.

Затем мы можем использовать функцию 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.

Проверка предположений модели

Прежде чем идти дальше, нам необходимо убедиться, что предположения нашей модели выполняются, чтобы результаты нашей модели были надежными. В частности, двусторонний дисперсионный анализ предполагает:

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() для выполнения теста Тьюки для множественных сравнений:

 #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).

Мы также можем визуализировать 95% доверительные интервалы, полученные в результате теста Тьюки, с помощью функцииplot() в R:

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

Отчет о результатах двустороннего дисперсионного анализа

Наконец, мы можем сообщить о результатах двустороннего дисперсионного анализа таким образом, чтобы суммировать результаты:

Двусторонний дисперсионный анализ был проведен для изучения влияния пола ( мужской, женский) и программы упражнений (нет, легкие, интенсивные) на потерю веса (измеряется в фунтах). Было статистически значимое взаимодействие между влиянием пола и физических упражнений на потерю веса (F(2, 54) = 4,615, p = 0,0141). Апостериорно были проведены тесты Тьюки HSD.

У мужчин интенсивная программа упражнений привела к значительно большей потере веса, чем легкая программа (p <0,0001) или отсутствие программы упражнений (p <0,0001). Кроме того, у мужчин легкая диета привела к значительно большей потере веса, чем отсутствие режима физических упражнений (p <0,0001).

У женщин интенсивная программа упражнений привела к значительно большей потере веса, чем легкая программа (p <0,0001) или отсутствие программы упражнений (p <0,0001).

Проверки нормальности и тест Левена были выполнены для проверки того, что предположения ANOVA выполнены.