ANOVA, или кто комментирует?

В комментариях проскальзывала мысль, что люди мало комментируют статьи на Habrahabr, т.к. боятся потерять карму. Получается, что в основном пишут те, у кого карма побольше. Попробуем исследовать эту гипотезу подробнее и получить результаты, подкрепленные не только интуитивно, но и статистически ^[1].

Нам необходимо проверить, оказывает ли карма пользователя статистически значимое ^[2] влияние на количество комментариев, которое он в среднем оставляет. Т.к. количество сравниваемых групп будет больше двух, то t-тест ^[3] не подойдет, и придется использовать дисперсионный анализ ^[4] — именно так расшифровывается ANOVA (analysis of variance).

Я воспользуюсь данными, которые ранее собрал varagian ^[5] и выложил тут ^[6]:

user_data <- read.csv('user_dataset.csv', stringsAsFactors=F, na.strings=c("", "NA"))

Для однофакторного дисперсионного анализа ^[7] понадобятся две переменные:

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

   И такое распределение — не самое удачное для дисперсионного анализа, т.к. для его проведения должны выполняться некоторые предпосылки, как, например, нормальность зависимой переменной. К счастью, в данном случае переменную можно «сделать» почти нормальной с помощью ее лог-транформации:

   comments_log <- log1p(user_data$comments)
   

   

2. Факторная переменная, влияние которой на зависимую переменную и исследуется. Посмотрим на распределение кармы:

   summary(user_data$karma)
      Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0.00    0.00    5.00   17.92   18.00 1230.00 
   

   Чтобы сгруппировать данные, введем новую переменную, которая будет представлять собой интервальный «срез» кармы и играть роль факторной переменной:

   karma_cut <- cut(user_data$karma, breaks=c(-Inf, 0, 5, 15, 25, 50, 100, Inf))
   table(karma_cut)
     (-Inf,0]      (0,5]     (5,15]    (15,25]    (25,50]   (50,100] (100, Inf] 
         5488       2059       2955       1423       1411        859        480 
   

   Самая многочисленная группа — это пользователи с кармой меньше или равной 0.

Nota bene!

Тут я должен сказать несколько слов о преобразовании непрерывной переменной в категориальную. Такая практика широко распространена в социологических, психологических и медико-биологических исследованиях: например, величину артериального давления можно условно определить как пониженное, нормальное, высокое. Но в целом с точки зрения статистики эта процедура несколько порочна, т.к. любое разбиение весьма условно, что потенциально ведет к потере информации.

Что же касается зависимости «количество комментариев» ~ «карма», то тут есть небольшая положительная корреляция, а линейная регрессия, выполненная на основе этих двух численных показателей (см. выше), являясь значимой, выглядит неубедительно, чтобы на ее основе делать какие-то статистические выводы: например, RESET-тест Рамсея ^[8] сигнализирует о пропущенных переменных, а тест Бройша-Пагана ^[9] — о гетероскедастичности ^[10] случайных ошибок. Кроме того, я заранее ставлю задачу сравнить группы пользователей, у которых карма воспринимается как «маленькая», «средняя» и.т.д.

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

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

Альтернативная гипотеза, соответственно, утверждает, что различия все же не случайны. Чтобы принять или отклонить нулевую гипотезу, нам надо сравнить межгрупповую VAR_b и внутригрупповую VAR_w дисперсии. Обе эти величины по-своему оценивают дисперсию генеральной совокупности, и при верной нулевой гипотезе их отношение находится недалеко от 1, т.е. внутригрупповая и межгрупповая дисперсии не различаются. Формулы для вычисления этих дисперсий приведены ниже:

Тут K — количество групп, N — общий объем выборки. Теперь надо оценить соотношение

В этом контексте принято говорить, что величина F следует F-распределению ^[11] со степенями свободы K-1 и N-K. Небольшой фрагмент кода на R, который вычисляет статистику F_test и критическое значение F_crit при уровне значимости α=5%:

alpha <- 0.05 
K <- length(levels(karma_cut))
N <- nrow(user_data)
df1 <- K - 1
df2 <- N - K
M <- mean(comments_log)

var_b <- sum(tapply(comments_log, karma_cut, function(x){
  length(x) * (mean(x) - M) ^ 2})) / df1

var_w <- sum(tapply(comments_log,karma_cut, function(x){
  sum((x - mean(x))^2)})) / df2

Ftest <- var_b / var_w

Fcrit <- qf(alpha, df1, df2, lower.tail=F)

Ниже на графике показано распределение F(6, 14668), критическое и вычисленное значения статистик:

Очевидно, что F_crit<<F_test, а значит, гипотезу H₀ можно отвергнуть. В R, конечно, все вышеописанное уже давно реализовано, поэтому просто проверим вычисления:

m.aov <- aov(comments_log ~ karma_cut)
summary(m.aov)
               Df Sum Sq Mean Sq F value Pr(>F)    
karma_cut       6   5519   919.9   556.2 <2e-16 ***
Residuals   14668  24260     1.7                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Значение F-статистики 556.2, что совпадает с вычисленным ранее.

Post hoc

Расчеты показывают, что гипотеза H₀ о равенстве средних должна быть отклонена (это единственное, что мы можем сказать, выполнив ANOVA). Что делать с этим выводом дальше? Обычно после ANOVA выполняют post hoc анализ ^[12]: можно применить t-тест с поправкой Бонферрони ^[13] (в R ^[14] это pairwise.t.test ^[15]) или, например, HSD тест Тьюки ^[16], чтобы выяснить, между какими группами существуют различия. Последний тест мне больше нравится:

TukeyHSD(m.aov, "karma_cut", conf.level=0.95)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = comments_log ~ karma_cut)

$karma_cut
                           diff         lwr        upr     p adj
(0,5]-(-Inf,0]      -0.08189943 -0.17990408 0.01610521 0.1725558
(5,15]-(-Inf,0]      0.29578008  0.20925195 0.38230821 0.0000000
(15,25]-(-Inf,0]     0.92766799  0.81485590 1.04048009 0.0000000
(25,50]-(-Inf,0]     1.36762030  1.25442791 1.48081269 0.0000000
(50,100]-(-Inf,0]    1.66955943  1.53041195 1.80870691 0.0000000
(100, Inf]-(-Inf,0]  1.87283311  1.69233134 2.05333488 0.0000000
(5,15]-(0,5]         0.37767952  0.26881660 0.48654243 0.0000000
(15,25]-(0,5]        1.00956743  0.87883646 1.14029840 0.0000000
(25,50]-(0,5]        1.44951973  1.31846045 1.58057902 0.0000000
(50,100]-(0,5]       1.75145887  1.59742628 1.90549145 0.0000000
(100, Inf]-(0,5]     1.95473255  1.76252197 2.14694313 0.0000000
(15,25]-(5,15]       0.63188791  0.50952455 0.75425128 0.0000000
(25,50]-(5,15]       1.07184022  0.94912615 1.19455428 0.0000000
(50,100]-(5,15]      1.37377935  1.22678192 1.52077678 0.0000000
(100, Inf]-(5,15]    1.57705303  1.39043280 1.76367327 0.0000000
(25,50]-(15,25]      0.43995231  0.29748058 0.58242403 0.0000000
(50,100]-(15,25]     0.74189144  0.57803877 0.90574410 0.0000000
(100, Inf]-(15,25]   0.94516512  0.74499878 1.14533146 0.0000000
(50,100]-(25,50]     0.30193913  0.13782440 0.46605386 0.0000012
(100, Inf]-(25,50]   0.50521281  0.30483189 0.70559373 0.0000000
(100, Inf]-(50,100]  0.20327368 -0.01283281 0.41938017 0.0811265

Предварительно можно сделать вывод, что пользователи с кармой (-∞, 0] и (0, 5] оставляют примерно одинаковое log-количество комментариев, как и пользователи с кармой (50,100] и (100, ∞]. Остальные группы существенно отличаются.

Что не так?

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

При более детальном рассмотрении она оказывается весьма далекой от нормального распределения, на что намекает и тест Шапиро-Уилка ^[17]. Другой критерий — тест Бартлетта ^[18] — свидетельствует о том, что дисперсия в группах неоднородна:

bartlett.test(comments_log~karma_cut)

	Bartlett test of homogeneity of variances

data: comments_log by karma_cut
Bartlett's K-squared = 84.811, df = 6, p-value = 3.613e-16

С одной стороны — нарушение предпосылок ANOVA налицо, с другой — есть ряд работ, которые говорят, что дисперсионный анализ весьма устойчив ^[19] и к нарушению условия нормальности, и к неоднородности дисперсии. Если нарушение условий все же кого-то стесняют, то можно воспользоваться непараметрическим тестом Краскела-Уоллиса ^[20]. Тут в качестве нулевой гипотезы выступает предположение о равенстве медиан (а не средних!) во всех группах:

kruskal.test(comments_log ~ karma_cut)

	Kruskal-Wallis rank sum test

data:  comments_log by karma_cut
Kruskal-Wallis chi-squared = 2816.1, df = 6, p-value < 2.2e-16

Тест Краскела-Уоллиса также позволяет отклонить нулевую гипотезу. Т.к. критерий Краскела-Уоллиса непараметрический, то и post-hoc анализ тоже должен основываться на непараметрических тестах, как, например, U-критерий Манна-Уитни ^[21] с последующей коррекцией p-значений, тест Данна. Последний содержится в пакете R dunn.test ^[22]:

dunn.test::dunn.test(comments_log, karma_cut, method="bonferroni")

При 5% уровне значимости, как и тест Тьюки, критерий Данна говорит о том, что пользователи с кармой (-∞, 0] и (0, 5], и (50,100] и (100, ∞]. оставляют примерно одинаковое log-количество комментариев.

ANOVA и линейная регрессия

Говоря о дисперсионном анализе, нельзя не упомянуть о линейной регрессии. Более того, можно утверждать, что это практически одно и то же, т.к. модель дисперсионного анализа — частный случай линейной регрессии. Запишем линейную регрессию с факторной переменной и проверим ее коэффициенты:

m.lm <- lm(comments_log ~ karma_cut-1)
summary(m.lm)
Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
karma_cut(-Inf,0]    3.20288    0.01736  184.50   <2e-16 ***
karma_cut(0,5]       3.12098    0.02834  110.12   <2e-16 ***
karma_cut(5,15]      3.49866    0.02366  147.88   <2e-16 ***
karma_cut(15,25]     4.13055    0.03409  121.16   <2e-16 ***
karma_cut(25,50]     4.57050    0.03424  133.50   <2e-16 ***
karma_cut(50,100]    4.87244    0.04388  111.04   <2e-16 ***
karma_cut(100, Inf]  5.07571    0.05870   86.47   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.286 on 14668 degrees of freedom
Multiple R-squared:  0.8914,	Adjusted R-squared:  0.8913 
F-statistic: 1.719e+04 on 7 and 14668 DF,  p-value: < 2.2e-16

Теперь найдем среднее количество лог-комментариев по каждому кармическому срезу:

tapply(comments_log, karma_cut, mean)
(-Inf,0]      (0,5]     (5,15]    (15,25]    (25,50]   (50,100] (100, Inf] 
3.202879   3.120979   3.498659   4.130547   4.570499   4.872438   5.075712 

Коэффициенты в линейной регрессии — это и есть средние значения, т.е. справедлива такая формула:

где μ_i — среднее соответствующей группы, ε_ij — случайная ошибка, распределенная нормально.

Заключение

Итак, анализ дисперсии выполняется по такой схеме:

1. Определяются зависимая и факторная переменная.
2. Проверяются предпосылки для проведения ANOVA.
3. При нормальной зависимой переменной и однородной дисперсии выполняется ANOVA, иначе — непараметрический тест Краскела-Уоллиса
4. После анализа дисперсии и отвергания нулевой гипотезы, выполняется соответствующий рost-hoc анализ — параметрический или непараметрический.

При многофакторном анализе дисперсий ^[23] так же применяется F-тест ^[24], но формулы для расчета дисперсии усложняются; кроме того, возникают дополнительные термы, описывающие взаимодействие между факторами, и появляется несколько нулевых гипотез.

Что касается комментариев на Хабре, то и ANOVA, и тест Краскела-Уоллиса говорят, что есть зависимость между log-количеством комментариев и кармой пользователя, причем группы с кармой (-∞, 0] и (0, 5], и (50,100] и (100, ∞] оставляют примерно одинаковое количество комментариев. В целом, с ростом самой кармы увеличивается и количество комментариев. Интересно, что результаты и непараметрического, и параметрического тестов совпали, хотя предпосылки для проведения последнего были явно нарушены. Объяснить этот факт можно как робастностью дисперсионного анализа, так и тем, что при большом количестве данных (у нас есть тысячи измерений) нарушение предположения о нормальности данных оказывает гораздо меньшее влияние на результат, чем в случае с малой выборкой.

Автор: kxx

Источник ^[25]