Дисперсионный анализ (ANOVA)
Андрей Ловаков
03.02.2017
Установка пакетов и рабочей директории
Установим необходимые для занятия пакеты (это действие необходимо, только если пакеты ещё не установлены на ваш компьютер) (для установки выполните следующую строку без символа # вначале).
# install.packages(c('ggplot2', 'reshape2', 'dplyr', 'gplots', 'compute.es', 'lsr', 'pwr'))Подгружаем загруженные пакеты для их использования в текущей сессии работы R.
library(ggplot2)
library(reshape2)
library(dplyr)
library(gplots)
library(compute.es)
library(lsr)
library(pwr)Устанавливаем рабочую директорию: через меню или командой: Session -> Set working directory -> Choose directory
#setwd('...')
#(вместо '...' укажите путь к папке, в которой хранится файл)Открываем файл experiment.csv. В нём три переменные: pre - результаты теста интеллекта до приёма чудо-таблетки, post - результаты теста интеллекта после приёма чудо-таблетки, group - делит выборку на экспериментальную и контрольную (плацебо) группы.
df <- read.table("experiment.csv")
str(df)## 'data.frame': 200 obs. of 3 variables:
## $ group: Factor w/ 2 levels "control","experiment": 1 1 1 1 1 1 1 1 1 1 ...
## $ pre : int 102 78 105 103 111 96 91 129 110 122 ...
## $ post : int 83 100 95 90 126 96 108 57 121 95 ...Однофакторный дисперсионный анализ (One-way ANOVA)
Сначала нужно трансформировать данные из “широкого” формата в “длинный” формат. Подробнее про разницу в этих форматах см. http://seananderson.ca/2013/10/19/reshape.html. Воспользуемся функцией melt из пакета reshape2.
df_long <- melt(df)
df_long <- rename(df_long, measure=variable, IQ=value) # переименуем получившиеся переменныеПосмотрим на структуру новых данных.
str(df_long) ## 'data.frame': 400 obs. of 3 variables:
## $ group : Factor w/ 2 levels "control","experiment": 1 1 1 1 1 1 1 1 1 1 ...
## $ measure: Factor w/ 2 levels "pre","post": 1 1 1 1 1 1 1 1 1 1 ...
## $ IQ : int 102 78 105 103 111 96 91 129 110 122 ...С помощью однофакторного дисперсионного анализа можно сравнить средние значения зависимой переменной для двух и более групп, заданных категориальной группирующей переменной. Проверим, влияет ли номера группы на рост.
Запустим однофакторный дисперсионный анализ (One-way ANOVA), чтобы проверить влияние чудо-таблетки на уровень IQ.
fit <- aov(IQ ~ group, data=subset(df_long, measure=="post"))
summary(fit)## Df Sum Sq Mean Sq F value Pr(>F)
## group 1 10310 10310 44.86 2.14e-10 ***
## Residuals 198 45508 230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Дисперсионный анализ показал, что средний уровень IQ в повторном измерении оказался различным в экспериментальной и контрольной группах: F(1, 35) = 44.86, p < 0.001.
Посмотрим на средние в этих двух группах.
aggregate(df_long$IQ, by=list(df_long$group), FUN=mean)## Group.1 x
## 1 control 100.220
## 2 experiment 107.275Визуализируем результаты
Один из самых простых и быстрых вариантов визуализации - функция plotmeans() из пакета gplots&
plotmeans(IQ ~ group, data=subset(df_long, measure=="post"), connect=FALSE) # если connect=FALSE изменить на connect=TRUE, точки соединяться, это стоит делать, если фактор по оси X является ранговой или интервальной переменной (в нашем случае это не так, т.к. группа является качественной переменной)
Проверим, выполняется ли предположение о гомогенности дисперсии IQ в каждой из групп.
bartlett.test(IQ ~ group, data=subset(df_long, measure=="post"))##
## Bartlett test of homogeneity of variances
##
## data: IQ by group
## Bartlett's K-squared = 0.037649, df = 1, p-value = 0.8461Однофакторный дисперсионный анализ (One-way ANOVA) с более чем 2 уровнями фактора
Возьмём другой набор данных. В R есть встроенный набор данных ToothGrowth, содержащий данный о влиянии витамина C на рост зубов свинок. Переменная len - длина зубов, переменная dose - доза витамина C (было протестировано три варианта дозы: 0.5, 1, 2). Глянем на структуру данных.
str(ToothGrowth)## 'data.frame': 60 obs. of 3 variables:
## $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
## $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
## $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...Попробуем применить ANOVA, чтобы проверить, влияет ли дозировака витамина C на длину зубов.
fit2 <- aov(len ~ as.factor(dose), data=ToothGrowth) # Обратите внимание на то, что мы одновременно изменили тип переменной dose с количественной (num, см. выше) на номинальную, или качественную as.factor(dose). Это необходимо для того, чтобы с помощью этой переменной образовать три группы.
summary(fit2)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(dose) 2 2426 1213 67.42 9.53e-16 ***
## Residuals 57 1026 18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Дисперсионный анализ показал, что дозировака витамина C влияет на длину зубов F(2, 57) = 67.42, p-value < 0.001. Вспомним, что p-value = 9.53e-16 означает, что если нулевая гипотеза (H0: независимая переменная не влияет на зависимую) верна, то вероятность получить вычесленное значение критерия F* (в данном случае F = 67.42) равна 9.53e-16, что сильно меньше конвенционального значения 0.05. Следовательно, мы должны отклонить нулевую гипотезу об отсутствии эффекта и пинять альтернативную, т.е. эффект дозироввки на длину зубов есть. Однако у нас три уровня фактора дозировки (три группы). Различия есть между всеми тремя? Полученные результаты не дают возможности это утверждать, т.к. различия только между двуми из них могут показать значимый эффект. Чтобы понять, между какими именно группами (уровнчми фактора) есть различия, необходимо их попарно сравнить. Для этого существуют несколько тестов. В данном случае можно воспользоваться тестом Тьюки (Tukey’s HSD test).
TukeyHSD(fit2)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = len ~ as.factor(dose), data = ToothGrowth)
##
## $`as.factor(dose)`
## diff lwr upr p adj
## 1-0.5 9.130 5.901805 12.358195 0.00e+00
## 2-0.5 15.495 12.266805 18.723195 0.00e+00
## 2-1 6.365 3.136805 9.593195 4.25e-05Результатом применения теста является три попарных сравнения, трёх групп между собой. В колонке p adj приводятся значения p-value. В данном случае нулевая гипотеза о том, что между двумя группами нет различий. Поскольку во всех трёх случаях p-value меньше конвенционального значения 0.05, то различия в длине зубок есть между всеми тремя группами свинок.
Визуализируем результат с помощью более красивого рисунка. Но за красоту надо платить, поэтому придётся немного поработать руками.
Сначала рассчитаем средние значения длины зубов и их стандартные отклонения для каждой группы. Сделаем это с помощью очень полезной функции tapply. Почитайте про неё, например, здесь
means <- tapply(ToothGrowth$len, as.factor(ToothGrowth$dose), mean)
sd <- tapply(ToothGrowth$len, as.factor(ToothGrowth$dose), sd)
dose <- c("0.5", "1", "2")
plot_data <- data.frame(means, sd, dose) # объединим их в отдельный data.frame, добавив к нему групирующую переменную Посмотрим на получившийся data.frame
plot_data## means sd dose
## 0.5 10.605 4.499763 0.5
## 1 19.735 4.415436 1
## 2 26.100 3.774150 2Теперь на его основе построим стобликовую диаграмму с доверительными интервалами.
ggplot(plot_data, aes(x=dose, y=means)) +
geom_bar(stat="identity", width=.4, fill="gold") +
geom_errorbar(aes(ymin=means-sd, ymax=means+sd), width=.2) +
theme_bw()
Двухфакторный дисперсионный анализ (Two-way ANOVA)
С помощью двухфакторного дисперсионного анализа можно проверить влияние двух факторов на зависимую переменную. Продолжим работать с данными про свинок и их зубы. Попробуем проверить, влияет ли на длину зубов и дозировака витамина C, и тип введения его в организм свинки (supp: апельсиновый сок (OJ) или аскорбиновая кислота (VC)).
fit3 <- aov(len ~ as.factor(dose) + supp, data=ToothGrowth) # Обратите внимание на то, что мы одновременно изменили тип переменной dose с количественной (num, см. выше) на номинальную, или качественную as.factor(dose). Это необходимо для того, чтобы с помощью этой переменной образовать три группы.
summary(fit3)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(dose) 2 2426.4 1213.2 82.81 < 2e-16 ***
## supp 1 205.4 205.4 14.02 0.000429 ***
## Residuals 56 820.4 14.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Дисперсионный анализ (3x2 ANOVA) показал, что и дозировака витамина C влияет на длину зубов, F(2, 56) = 82.81, p < 0.001, и тип введения его в организм свинки, F(1, 56) = 14.02, p < 0.001.
Независимые переменные могут взаимодействовать дуг с другом, и их взаимодействие также может оказывать влияние на зависимую переменную. В предыдущий анализ кроме двух основных факторов можно добавить их взаимодействие. Давайте это сделаем.
fit4 <- aov(len ~ as.factor(dose) + supp + as.factor(dose):supp, data=ToothGrowth) # взаимодействие факторов добавляется через с помощью добавление в формула нового члена, состоящего из имен двух взаимодействующих факторов, разделённых двоеточием - as.factor(dose):supp
summary(fit4)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(dose) 2 2426.4 1213.2 92.000 < 2e-16 ***
## supp 1 205.4 205.4 15.572 0.000231 ***
## as.factor(dose):supp 2 108.3 54.2 4.107 0.021860 *
## Residuals 54 712.1 13.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Можно было сделать короче: **len ~ as.factor(dose)*supp**. Эта формула добавляет сразу и две переменные по-отдельности, и их взаимодействие. Результат получится тот же самый.
fit5 <- aov(len ~ as.factor(dose)*supp, data=ToothGrowth)
summary(fit5)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(dose) 2 2426.4 1213.2 92.000 < 2e-16 ***
## supp 1 205.4 205.4 15.572 0.000231 ***
## as.factor(dose):supp 2 108.3 54.2 4.107 0.021860 *
## Residuals 54 712.1 13.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Дисперсионный анализ (3x2 ANOVA) показал, что есть как два основных эффекта дозы и типа введения по-отдельности, F(2, 54) = 92.00, p < 0.001, F(1, 54) = 15.57, p < 0.001, так и эффект их взаимодействия, F(2, 54) = 4.107, p = 0.022.
Посмотрим, какие именно группу отличаются друг от друга с помощью теста Тьюки.
TukeyHSD(fit5)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = len ~ as.factor(dose) * supp, data = ToothGrowth)
##
## $`as.factor(dose)`
## diff lwr upr p adj
## 1-0.5 9.130 6.362488 11.897512 0.0e+00
## 2-0.5 15.495 12.727488 18.262512 0.0e+00
## 2-1 6.365 3.597488 9.132512 2.7e-06
##
## $supp
## diff lwr upr p adj
## VC-OJ -3.7 -5.579828 -1.820172 0.0002312
##
## $`as.factor(dose):supp`
## diff lwr upr p adj
## 1:OJ-0.5:OJ 9.47 4.671876 14.2681238 0.0000046
## 2:OJ-0.5:OJ 12.83 8.031876 17.6281238 0.0000000
## 0.5:VC-0.5:OJ -5.25 -10.048124 -0.4518762 0.0242521
## 1:VC-0.5:OJ 3.54 -1.258124 8.3381238 0.2640208
## 2:VC-0.5:OJ 12.91 8.111876 17.7081238 0.0000000
## 2:OJ-1:OJ 3.36 -1.438124 8.1581238 0.3187361
## 0.5:VC-1:OJ -14.72 -19.518124 -9.9218762 0.0000000
## 1:VC-1:OJ -5.93 -10.728124 -1.1318762 0.0073930
## 2:VC-1:OJ 3.44 -1.358124 8.2381238 0.2936430
## 0.5:VC-2:OJ -18.08 -22.878124 -13.2818762 0.0000000
## 1:VC-2:OJ -9.29 -14.088124 -4.4918762 0.0000069
## 2:VC-2:OJ 0.08 -4.718124 4.8781238 1.0000000
## 1:VC-0.5:VC 8.79 3.991876 13.5881238 0.0000210
## 2:VC-0.5:VC 18.16 13.361876 22.9581238 0.0000000
## 2:VC-1:VC 9.37 4.571876 14.1681238 0.0000058Визуализируем результаты
pd = position_dodge(0.4)
ggplot(ToothGrowth, aes(as.factor(dose), len, color = supp)) +
stat_summary(fun.data = mean_cl_boot, geom = 'errorbar', width = 0.2, lwd = 0.8, position = pd)+
stat_summary(fun.data = mean_cl_boot, geom = 'line', size = 1.5, position = pd) +
stat_summary(fun.data = mean_cl_boot, geom = 'point', size = 5, position = pd, pch=15) +
theme_bw() +
xlab('Доза витамина C')+
ylab('Длина')
Size effect
Рассчитаем размер эффекта типа введение Витамина C на длину зубов. Воспользуемся пакетом lsr, в котором есть функция etaSquared, рассчитывающая размер эффекта eta-Squared.
etaSquared(fit2, anova = TRUE)## eta.sq eta.sq.part SS df MS F
## as.factor(dose) 0.7028642 0.7028642 2426.434 2 1213.21717 67.41574
## Residuals 0.2971358 NA 1025.775 57 17.99605 NA
## p
## as.factor(dose) 8.881784e-16
## Residuals NAPазмер эффекта eta-Squared = 0.703. Это эффект cредней величины.
Анализ мощности
Теперь проверим, сколько необходимо людей, чтобы зафиксировать эффект такого размера. Воспользуется пакетом pwr, в котором есть функция pwr.anova.test (есть и для других дизайнов).
pwr.anova.test(k = 2, f= 0.703, sig.level = 0.05, power = 0.80)##
## Balanced one-way analysis of variance power calculation
##
## k = 2
## n = 9.010349
## f = 0.703
## sig.level = 0.05
## power = 0.8
##
## NOTE: n is number in each groupАнализ мощности показал, что для того, чтобы с 80% мощностью зафиксировать эффект (eta-Squared = 0.703), необходимо всего по 9 человек на группу.