Двухфакторный дисперсионный анализ («дисперсионный анализ») используется для определения того, существует ли статистически значимое различие между средними значениями трех или более независимых групп, разделенных на два фактора.
В этом руководстве объясняется, как выполнить двухфакторный дисперсионный анализ в 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) %>%
summarise (mean = mean(weight_loss),
sd = sd(weight_loss))
# A tibble: 6 x 4
# Groups: gender [2]
# gender exercise mean 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
par(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 выглядит следующим образом:
aov (переменная ответа ~ переменная_предиктора1 * переменная_предиктора2, данные = набор данных)
Обратите внимание, что * между двумя переменными-предикторами указывает, что мы также хотим проверить эффект взаимодействия между двумя переменными-предикторами.
В нашем примере мы можем использовать следующий код, чтобы соответствовать двухфакторной модели ANOVA, используя weight_loss в качестве переменной ответа и пол и упражнения в качестве двух наших переменных-предикторов.
Затем мы можем использовать функцию 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
#---
#Signif. codes: 0 '\*\*\*' 0.001 '\*\*' 0.01 '\*' 0.05 '.' 0.1 ' ' 1
Из результатов модели мы видим, что пол , физические упражнения и взаимодействие между двумя переменными статистически значимы на уровне значимости 0,05.
Проверка предположений модели
Прежде чем идти дальше, мы должны проверить, выполняются ли предположения нашей модели, чтобы наши результаты из модели были надежными. В частности, двусторонний дисперсионный анализ предполагает:
1. Независимость – наблюдения в каждой группе должны быть независимыми друг от друга. Поскольку мы использовали рандомизированный дизайн , это предположение должно выполняться, поэтому нам не нужно слишком беспокоиться об этом.
2. Нормальность – зависимая переменная должна быть примерно нормально распределена для каждой комбинации групп двух факторов.
Один из способов проверить это предположение — построить гистограмму остатков модели. Если остатки примерно нормально распределены, это предположение должно выполняться.
#define model residuals
resid <- 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: 0,990364).
Мы также можем визуализировать 95% доверительные интервалы, полученные в результате теста Тьюки, с помощью функции plot() в R:
#set axis margins so labels don't get cut off
par(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 были выполнены.