Сравнение трех и более групп

Для сравнения трех и более групп не подходят уже известные нам t-критерий Стьюдента и критерий Вилкоксона. Мы будем пользоваться параметрическим методом, который называется однофакторный дисперсионный анализ – ANOVA (analysis of variance) и непараметрическим ранговым критерием Краскелла-Уоллиса. Выбор между критериями осуществляется на основании знания о виде распределения. Если в большинстве групп данные имеют нормальное распределение, то мы можем использовать дисперсионный анализ. Если нормальности в данных нет, мы выбираем критерий Краскелла-Уоллиса.

ANOVA

Однофакторный дисперсионный анализ позволяет нам проверить гипотезу о равенстве средних значений некоторой переменной между группами. Нулевая гипотеза предполагает, что все средние одинаковы:

\[H_{0}: \mu_{1}=\mu_{2}=...=\mu_{k}\]

Альтернативная гипотеза предполагает, что существует хотя бы одно отличающееся среднее:

\[H_{1}: хотя~бы~одно~среднее~отличается\]

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

  • Внутригрупповая дисперсия показывает, насколько сильно различаются наблюдения внутри одной группы. Иными словами, насколько сильно каждое из наблюдений отклоняется от внутригруппового среднего. Закономерно предположить, что сам факт принадлежности к той или иной группе не будет оказывать влияния на внутригрупповую дисперсию.

\[ \displaystyle MSW = \frac{\sum_{j=1}^{k} \sum_{i=1}^{N} (y_{ij}-\bar{y}_{j})^2}{N-k}~,\] где \(MSW\) – внутригрупповой средний квадрат (отклонений), \(N\) – количество наблюдений во всех группах вместе, \(k\) – количество групп.

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

\[ \displaystyle MSB = \frac{\sum_{j=1}^k n_{j} (\bar{y}_{j}-\bar{y})^2}{k-1}~,\] где \(MSB\) – межгрупповой средний квадрат (отклонений), \(n_j\) – количество наблюдений в группе \(j\), \(k\) – общее количество групп.

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

\[(y_{ij}-\bar{y_{j}})+(\bar{y_{j}} - \bar{y})=y_{ij}-\bar{y}\]

Статистический критерий, который используется для проверки нулевой гипотезы – F-критерий – расчитывается как отношение межгрупповой дисперсии к внутригрупповой дисперсии. Чем больше межгрупповая дисперсия будет по отношению к внутригрупповой дисперсии, тем более правдоподобно, что группы различаются между собой. Статистика критерия имеет распределение Фишера с \(k-1\) и \(N-k\) степенями свободы.

Критерий Краскелла-Уоллиса

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

\[H_{0}: F_1(x)=F_2(x)=...=F_k(x)\]

\[H_{1}: хотя~бы~одно~распределение~отличается\] Для расчета статистики критерия строится объединенный вариационный ряд по всем наблюдениям всех групп и вычисляются ранги наблюдений в этом объединенном вариационном ряду. После этого вычисляется средний ранг по группам. Формула статистики критерия выглядит следующим образом и имеет распределение хи-квадрат с \(k-1\) степенями свободы:

\[H=\frac{12}{N(N+1)}\sum_{j=1}^k \frac{R_j^2}{n_j} - 3(N+1)~,\] где \(R_{j}^2\) – квадрат суммы рангов в группе \(j\), \(n_{j}\) – количество наблюдений в группе \(j\), \(k\) – общее количество групп, \(N\) – общее количество наблюдений.

Рост студентов в зависимости от любимого предмета

Подгрузим необходимые пакеты. Сегодня нам понадобятся пакеты psych, ggplot2 и car.

library(psych) # описательные статистики по группам
library(ggplot2) # графики по группам
library(car) # проверка однородности дисперсий

Откроем базу данных.

df <- read.csv2("http://math-info.hse.ru/f/2018-19/psych-ms/data/psych_survey.csv")
head(df)
##   height math bio subject gender  residence length angle soft   subject2
## 1  169.0   76  96       1   male Not Moscow     NA    NA    R   Sciences
## 2  168.0   62  79       2 female     Moscow     NA    NA    R   Sciences
## 3  164.0   80  94       2 female Not Moscow     18    11    R   Sciences
## 4  158.0   75  77       4 female     Moscow     43    12    R Humanities
## 5  186.0   62  74       2   male Not Moscow     20    15    R   Sciences
## 6  168.5   70  65       1 female     Moscow     25    15    R   Sciences

Нам понадобятся переменные height – рост студента – и subject2 – любимый предмет в школе. Посмотрим на уровни фактора переменной любимого предмета. Она принимает три значения и содержит два пропущенных значения:

summary(df$subject2)
## Humanities    Neither   Sciences       NA's 
##         36         19         66          2

Избавимся от пропущенных значений в переменной любимого предмета:

df <- subset(df, subject2 != "NA")

Выведем описательные статистики для роста по всей выборке:

summary(df$height)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     145     163     169     169     175     190

Мы видим, что в нашей выборке нет пропущенных наблюдений.

Выведем описательные статистики по группам:

describeBy(df$height, group = df$subject2)
## 
##  Descriptive statistics by group 
## group: Humanities
##    vars  n   mean   sd median trimmed mad min max range  skew kurtosis
## X1    1 36 166.03 8.69  164.5   166.2 8.9 145 183    38 -0.12    -0.48
##      se
## X1 1.45
## -------------------------------------------------------- 
## group: Neither
##    vars  n   mean   sd median trimmed  mad min max range skew kurtosis
## X1    1 19 168.34 8.44    168  168.03 5.93 156 186    30 0.28    -0.75
##      se
## X1 1.94
## -------------------------------------------------------- 
## group: Sciences
##    vars  n   mean   sd median trimmed  mad min max range skew kurtosis
## X1    1 66 170.77 8.46    170  170.46 8.15 150 190    40 0.29    -0.35
##      se
## X1 1.04

Мы видим, что средний рост студентов, указавших в качестве любимого гуманитарный предмет, равен \(166\) см. Средний рост студентов, указавших в качестве любимого естественно-научный предмет, равен \(170.7\) см. Средний рост студентов, которые указали, что им не нравится ни то, ни другое, составил \(168.34\). Напомню, что значения, которые мы получили – это выборочные оценки средних, которые имеют свою изменчивость (об изменчивости оценки среднего свидетельствует его стандартная ошибка, se, которая расчитывается как \(\frac{sd}{\sqrt{n}}\)). Но можем ли мы утверждать, что существуют значимые различия между группами относительно роста?

Дополнительно: красивые ящики с усами

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

ggplot(data = df, aes(y = height, x = subject2, fill = subject2)) + 
  geom_boxplot() + 
  xlab("Favorite Subject") + ylab("Height")

Мы можем отметить на ящиках с усами групповые средние, используя функцию stat_summary(fun.y=mean, geom="point"):

ggplot(data = df, aes(y = height, x = subject2, fill = subject2)) + 
  geom_boxplot() + 
  xlab("Favorite Subject") + ylab("Height")+ 
  stat_summary(fun.y=mean, geom="point", size=5, color="red", fill="red")

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

ggplot(data = df, aes(y = height, x = subject2, fill = subject2)) + 
  geom_boxplot() + 
  xlab("Favorite Subject") + ylab("Height") +
  stat_summary(fun.y=mean, geom="point", shape=20, size=10, color="red", fill="red") +
  geom_jitter()

Проверка на нормальность

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

hum <- subset(df, subject2 == "Humanities")
sci <- subset(df, subject2 == "Sciences")
nei <- subset(df, subject2 == "Neither")

Проверим нормальность с помощью теста Колмогорова-Смирнова. Напомню, что тест сравнивает теоретическую и эмпирическую функцию распределения. Чтобы сравнение выполнялось правильно, надо указать R, какие параметры имеет теоретическое распределение: точно так же, как мы накладываем кривую нормального распределения на гистограмму.

ks.test(hum$height, "pnorm",
        mean = mean(hum$height, na.rm = T),
        sd = sd(hum$height, na.rm = T))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  hum$height
## D = 0.10263, p-value = 0.8427
## alternative hypothesis: two-sided
ks.test(sci$height, "pnorm",
        mean = mean(sci$height, na.rm = T),
        sd = sd(sci$height, na.rm = T))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sci$height
## D = 0.12729, p-value = 0.2352
## alternative hypothesis: two-sided
ks.test(nei$height, "pnorm",
        mean = mean(nei$height, na.rm = T),
        sd = sd(nei$height, na.rm = T))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  nei$height
## D = 0.10631, p-value = 0.9827
## alternative hypothesis: two-sided

Проверим то же самое с помощью критерия Шапиро-Уилка:

shapiro.test(hum$height)
## 
##  Shapiro-Wilk normality test
## 
## data:  hum$height
## W = 0.98039, p-value = 0.7583
shapiro.test(sci$height)
## 
##  Shapiro-Wilk normality test
## 
## data:  sci$height
## W = 0.97645, p-value = 0.2415
shapiro.test(nei$height)
## 
##  Shapiro-Wilk normality test
## 
## data:  nei$height
## W = 0.95783, p-value = 0.5304

Согласно всем тестам, во всех выборках нулевая гипотеза о нормальном распределении не отвергается. Следовательно, мы выбираем однофакторный дисперсионный анализ для сравнения средних в нескольких группах.

ANOVA

Чтобы провести однофакторный дисперсионный анализ в R в предположении о равных дисперсиях, нужно воспользоваться командой aov(), которая имеет такое же устройство, как и команда для линейной регрессии. Нужно сохранить результаты в отдельный объект, а потом вывести для них описание.

anova <- aov(height ~ subject2, data = df)
summary(anova)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## subject2      2    534   266.8    3.67 0.0284 *
## Residuals   118   8578    72.7                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Перед нами появилась стандартная ANOVA-таблица. Давайте обсудим ее устройство. В первой строчке находится межгрупповая дисперсия, во второй – внутригрупповая дисперсия.

  • В первой колонке указаны степени свободы df;
  • во второй колонке – сумма квадратов отклонений, Sum Sq;
  • в третьей колонке – средний квадрат, Mean Sq, который расчитывается как отношение суммы квадратов к степеням свободы;
  • в четвертой колонке – значение F-статистики, отношение междгруппового среднего квадрата к внутригрупповому;
  • в пятой колонке – p-value для F-статистики, Pr(>F).

Если мы сложим межгрупповую сумму квадратов и внутригрупповую сумму квадратов, то получим общую сумму квадратов: \(\sum(y_{i}-\bar{y})^2\). Общая сумма квадратов – это числитель выборочной дисперсии для объединенной выборки. Давайте это проверим.

ssb<-534
ssw<-8578
ssb+ssw
## [1] 9112
var(df$height)
## [1] 75.93082
var(df$height)*120 # умножаем на знаменатель, n-1, чтобы получить числитель
## [1] 9111.698

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

Апостериорные сравнения

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

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

# multiple pairwise comparison 
TukeyHSD(anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = height ~ subject2, data = df)
## 
## $subject2
##                         diff        lwr      upr     p adj
## Neither-Humanities  2.314327 -3.4245208 8.053176 0.6052814
## Sciences-Humanities 4.744949  0.5517231 8.938176 0.0223743
## Sciences-Neither    2.430622 -2.8384258 7.699670 0.5190756

Мы видим, что существует значимое различие между ростом студентов, кто выбрал в качестве любимого предмета точные науки, и студентов, кто выбрал в качестве любимого предмета гуманитарные науки. Первые выше.

Можно визуализировать эти различия:

plot(TukeyHSD(anova))

Случай неравных дисперсий

Как и в t-критерии Стьюдента, процедура дисперсионного анализа предполагает, что дисперсии между группами равны. Мы можем проверить это предположение с помощью критерия Ливиня (нулевая гипотеза: все дисперсии равны; альтернативная гипотеза: хотя бы одна отличается). Для реализации команды нужен пакет car:

# car package required
leveneTest(height ~ subject2, data = df)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.0875 0.9163
##       118

В нашем случае дисперсии дейстительно оказались равными (\(H_{0}\) не отвергается). Если же нулевая гипотеза была отвергнута, то следует реализовывать дисперсионный анализ с поправкой на неравные дисперсии. Это делается следующей командой:

oneway.test(height ~ subject2, data = df, var.equal = F)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  height and subject2
## F = 3.5608, num df = 2.000, denom df = 46.952, p-value = 0.03632

Критерий Краскелла-Уоллиса

Реализуем критерий Краскелла-Уоллиса:

kruskal.test(height ~ subject2, data = df)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  height by subject2
## Kruskal-Wallis chi-squared = 5.7817, df = 2, p-value = 0.05553

Мы видим, что значение p-value оказалось пограничным. В зависимости от выбранного уровня значимости, мы можем как отвергнуть, так и не отвергать нулевую гипотезу о равенстве распределний между группами.