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

Узнайте, как выполнить двусторонний ANOVA в R. Вы также узнаете его цель, гипотезы, предположения и способы интерпретации результатов.

Введение

Двусторонний ANOVA (дисперсионный анализ) — это статистический метод, который позволяет оценить одновременное влияние двух категориальных переменных на количественно непрерывную переменную.

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

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

• Корреляция измеряет взаимосвязь между двумя количественными переменными. Множественная линейная регрессия также измеряет взаимосвязь между двумя переменными, но на этот раз с учетом потенциального влияния других ковариат.
• Односторонний ANOVA проверяет, различаются ли количественные переменные между группами. Двусторонний дисперсионный анализ также проверяет, различается ли количественная переменная между группами, но на этот раз с учетом влияния другой качественной переменной.

Ранее мы обсуждали однофакторный дисперсионный анализ в R. Теперь мы покажем, когда, почему и как выполнять двусторонний дисперсионный анализ в R.

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

Критерий Стьюдента используется для оценки влияния одной категориальной переменной на количественную непрерывную переменную, когда категориальная переменная имеет ровно 2 уровня:

• Критерий Стьюдента для независимых выборок, если наблюдения независимы (например, если мы сравниваем возраст женщин и мужчин).
• критерий Стьюдента для парных выборок, если наблюдения зависимы, то есть когда они идут парами (это случай, когда одни и те же испытуемые измеряются дважды, в два разные моменты времени, например, до и после лечения)

Чтобы оценить влияние одной категориальной переменной на количественную переменную, когда категориальная переменная имеет 3 или более уровней: 1

• однофакторный дисперсионный анализ (часто называемый просто дисперсионным анализом), если группы независимы (например, группа пациентов, получавших лечение А, другая группа пациентов, получавших лечение B и последняя группа пациентов, не получавших лечения или получавших плацебо)
• повторные измерения ANOVA, если группы зависимы (когда одни и те же субъекты измеряются три раза, в три разных момента времени, например, до, во время и после лечения)

Двухфакторный дисперсионный анализ используется для оценки влияния двух категориальных переменных (и их потенциального взаимодействия) на количественную непрерывную переменную. Это тема поста.

Линейная регрессия используется для оценки взаимосвязи между количественной непрерывной зависимой переменной и одной или несколькими независимыми переменными:

• простая линейная регрессия, если есть только одна независимая переменная (которая может быть количественной или качественной)
• множественная линейная регрессия, если есть как минимум две независимые переменные (которые могут быть количественными, качественными или их сочетанием)

ANCOVA (анализ ковариации) используется для оценки влияния категориальной переменной на количественную переменную с учетом влияния другой количественной переменной (известной как ковариация). ANCOVA на самом деле является частным случаем множественной линейной регрессии с сочетанием одной качественной и одной количественной независимых переменных.

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

Данные

Чтобы проиллюстрировать, как выполнить двусторонний ANOVA в R, мы используем набор данных penguins, доступный в пакете {palmerpenguins}.

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

# install.packages("palmerpenguins")
library(palmerpenguins)
dat <- penguins # rename dataset
str(dat) # structure of dataset
## tibble [344 × 8] (S3: tbl_df/tbl/data.frame)
##  $ species          : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ island           : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ bill_length_mm   : num [1:344] 39.1 39.5 40.3 NA 36.7 39.3 38.9 39.2 34.1 42 ...
##  $ bill_depth_mm    : num [1:344] 18.7 17.4 18 NA 19.3 20.6 17.8 19.6 18.1 20.2 ...
##  $ flipper_length_mm: int [1:344] 181 186 195 NA 193 190 181 195 193 190 ...
##  $ body_mass_g      : int [1:344] 3750 3800 3250 NA 3450 3650 3625 4675 3475 4250 ...
##  $ sex              : Factor w/ 2 levels "female","male": 2 1 1 NA 1 2 1 2 NA NA ...
##  $ year             : int [1:344] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...

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

summary(dat)
##       species          island    bill_length_mm  bill_depth_mm  
##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
##  Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  
##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
##                                  Mean   :43.92   Mean   :17.15  
##                                  3rd Qu.:48.50   3rd Qu.:18.70  
##                                  Max.   :59.60   Max.   :21.50  
##                                  NA's   :2       NA's   :2      
##  flipper_length_mm  body_mass_g       sex           year     
##  Min.   :172.0     Min.   :2700   female:165   Min.   :2007  
##  1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  
##  Median :197.0     Median :4050   NA's  : 11   Median :2008  
##  Mean   :200.9     Mean   :4202                Mean   :2008  
##  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  
##  Max.   :231.0     Max.   :6300                Max.   :2009  
##  NA's   :2         NA's   :2

В этом посте мы сосредоточимся на следующих трех переменных:

• species: вид пингвина (Адели, Антарктический или Генту)
• sex: пол пингвина (самка и самец)
• body_mass_g: масса тела пингвина (в граммах).

При необходимости дополнительную информацию об этом наборе данных можно найти, запустив ?penguins в R.

body_mass_g является количественной непрерывной переменной и будет зависимой переменной, тогда как species и sex являются качественными переменными.

Эти две последние переменные будут нашими независимыми переменными, также называемыми факторами. Убедитесь, что R читает их как факторы. Если это не так, их нужно будет преобразовать в факторы.

Цель и гипотезы двустороннего дисперсионного анализа

Как упоминалось выше, двусторонний дисперсионный анализ используется для одновременной оценки влияния двух категориальных переменных на одну количественную непрерывную переменную.

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

Здесь мы хотели бы знать, зависит ли масса тела от вида и/или пола. В частности, нас интересует:

1. измерение и тестирование взаимосвязи между видами и массой тела,
2. измерение и тестирование взаимосвязи между полом и массой тела, а также
3. потенциально проверить, отличается ли взаимосвязь между видом и массой тела для самок и самцов (что эквивалентно проверке того, зависит ли взаимосвязь между полом и массой тела от вида)

Первые две взаимосвязи называются основными эффектами, а третья точка известна как эффект взаимодействия.

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

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

Если вернуться к нашему примеру, то мы имеем следующие проверки гипотез:

Основное влияние секса на массу тела:

• H0: средняя масса тела одинакова у самок и самцов.
• H1: средняя масса тела у самок и самцов разная.

Основное влияние видов на массу тела:

• H0: средняя масса тела одинакова для всех 3 видов.
• H1: средняя масса тела отличается по крайней мере для одного вида.

Взаимодействие между полом и видом:

• H0: взаимодействие между полом и видом отсутствует, что означает, что взаимосвязь между видом и массой тела одинакова для самок и самцов (аналогично, взаимосвязь между полом и массой тела одинакова для всех 3 видов)
• H1: существует взаимодействие между полом и видом, что означает, что связь между видом и массой тела у самок иная, чем у мужчин (аналогично, связь между полом и массой тела зависит от вида)

Предположения двустороннего дисперсионного анализа

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

Допущения двухфакторного дисперсионного анализа такие же, как и для однофакторного дисперсионного анализа. Обобщить:

• Тип переменной: зависимая переменная должна быть количественно непрерывной, а две независимые переменные должны быть категориальными (как минимум с двумя уровнями).
• Независимость: наблюдения должны быть независимыми между группами и внутри каждой группы.
• Нормальность:
• Для небольших выборок данные должны следовать примерно нормальному распределению.
• Для больших выборок (обычно n ≥ 30 в каждой группе/выборке) нормальность не требуется (благодаря центральной предельной теореме).
• Равенство различий: различия должны быть одинаковыми для всех групп.
• Выбросы: ни в одной группе не должно быть значительных выбросов.

Подробнее об этих предположениях можно узнать в разделе Предположения однофакторного дисперсионного анализа.

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

Тип переменной

Зависимая переменная масса тела является количественно непрерывной, тогда как независимые переменные пола и вида являются качественными переменными (имеющими не менее 2 уровней).

Следовательно, это предположение выполняется.

Независимость

Независимость обычно проверяют на основе плана эксперимента и способа сбора данных.

Для простоты наблюдения обычно таковы:

• независимый, если каждая экспериментальная единица (здесь пингвин) была измерена только один раз и наблюдения получены от репрезентативной и случайно выбранной части популяции, или
• зависит, если каждая экспериментальная единица была измерена как минимум дважды (как это часто бывает в области медицины, например, с двумя измерениями на одних и тех же субъектах; одно до и одно после лечения).

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

Нормальность

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

table(dat$species, dat$sex)
##            
##             female male
##   Adelie        73   73
##   Chinstrap     34   34
##   Gentoo        58   61

так что нормальность проверять не надо.

Для полноты мы еще покажем, как проверить нормальность, как если бы у нас были небольшие выборки.

Существует несколько методов проверки предположения о нормальности. Наиболее распространены следующие методы:

• QQ-график по группам или по остаткам, и/или
• гистограмма по группам или по остаткам, и/или
• тест на нормальность (например, тест Шапиро-Уилка) по группам или по остаткам.

Самый простой/кратчайший способ — проверить нормальность с помощью QQ-графика остатков. Чтобы нарисовать этот график, нам сначала нужно сохранить модель:

# save model
mod <- aov(body_mass_g ~ sex * species,
  data = dat
)

Этот фрагмент кода будет объяснен далее.

Теперь мы можем нарисовать QQ-график по остаткам. Мы покажем два способа сделать это, первый с помощью функции plot(), а второй с помощью функции qqPlot() из пакета {car}:

# method 1
plot(mod, which = 2)

# method 2
library(car)
qqPlot(mod$residuals,
  id = FALSE # remove point identification
)

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

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

Если вы предпочитаете проверять нормальность на основе гистограммы остатков, вот код:

# histogram
hist(mod$residuals)

Гистограмма остатков показывает гауссово распределение, которое соответствует выводу из графика QQ.

Хотя QQ-графика и гистограммы в значительной степени достаточно для проверки нормальности, если вы хотите проверить это более формально с помощью статистического теста, тест Шапиро-Уилка также может быть применен к остаткам:

# normality test
shapiro.test(mod$residuals)
## 
## 	Shapiro-Wilk normality test
## 
## data:  mod$residuals
## W = 0.99776, p-value = 0.9367

⇒ Мы не отвергаем нулевую гипотезу о том, что остатки следуют нормальному распределению (p-значение = 0,937).

Из QQ-графика, гистограммы и критерия Шапиро-Уилка мы заключаем, что не отвергаем нулевую гипотезу о нормальности остатков.

Таким образом, предположение о нормальности проверено, теперь мы можем проверить равенство дисперсий. 2

Однородность отклонений

Равенство дисперсий, также называемое однородностью дисперсий или гомоскедастичностью, можно проверить визуально с помощью функции plot():

plot(mod, which = 3)

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

Приведенного выше диагностического графика достаточно, но при желании его также можно протестировать более формально с помощью теста Левена (также из пакета {car}): 3

leveneTest(mod)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   5  1.3908 0.2272
##       327

⇒ Мы не отвергаем нулевую гипотезу о том, что дисперсии равны (p-значение = 0,227).

И визуальный, и формальный подходы дают один и тот же вывод; мы не отвергаем гипотезу об однородности дисперсий.

Выбросы

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

Для женщин и мужчин:

library(ggplot2)
# boxplots by sex
ggplot(dat) +
  aes(x = sex, y = body_mass_g) +
  geom_boxplot()

Для трех видов:

# boxplots by species
ggplot(dat) +
  aes(x = species, y = body_mass_g) +
  geom_boxplot()

Согласно критерию межквартильного диапазона, у вида Chinstrap есть два выброса. Эти точки, тем не менее, не являются достаточно экстремальными, чтобы исказить результаты.

Поэтому мы считаем, что предположение об отсутствии значительных выбросов выполняется.

Двусторонний дисперсионный анализ

Мы показали, что все предположения выполняются, поэтому теперь мы можем перейти к реализации двухфакторного дисперсионного анализа в R.

Это позволит нам ответить на следующие вопросы исследования:

• С поправкой на вид, значительно ли различается масса тела у представителей двух полов?
• С учетом пола, значительно ли различается масса тела хотя бы у одного вида?
• Отличается ли взаимосвязь между видом и массой тела у самок и самцов пингвинов?

Предварительные анализы

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

Это можно сделать с помощью описательной статистики или графиков.

Описательная статистика

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

# mean by group
aggregate(body_mass_g ~ species + sex,
  data = dat,
  FUN = mean
)
##     species    sex body_mass_g
## 1    Adelie female    3368.836
## 2 Chinstrap female    3527.206
## 3    Gentoo female    4679.741
## 4    Adelie   male    4043.493
## 5 Chinstrap   male    3938.971
## 6    Gentoo   male    5484.836

Или, наконец, среднее значение и стандартное отклонение для каждой подгруппы с использованием пакета {dplyr}:

# mean and sd by group
library(dplyr)
group_by(dat, sex, species) %>%
  summarise(
    mean = round(mean(body_mass_g, na.rm = TRUE)),
    sd = round(sd(body_mass_g, na.rm = TRUE))
  )
## # A tibble: 8 × 4
## # Groups:   sex [3]
##   sex    species    mean    sd
##   <fct>  <fct>     <dbl> <dbl>
## 1 female Adelie     3369   269
## 2 female Chinstrap  3527   285
## 3 female Gentoo     4680   282
## 4 male   Adelie     4043   347
## 5 male   Chinstrap  3939   362
## 6 male   Gentoo     5485   313
## 7 <NA>   Adelie     3540   477
## 8 <NA>   Gentoo     4588   338

Сюжеты

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

Наиболее подходящим графиком, когда у нас есть одна количественная и две качественные переменные, является блочная диаграмма по группам. Это легко сделать с помощью пакета {ggplot2}:

# boxplot by group
library(ggplot2)
ggplot(dat) +
  aes(x = species, y = body_mass_g, fill = sex) +
  geom_boxplot()

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

dat %>%
  filter(!is.na(sex)) %>%
  ggplot() +
  aes(x = species, y = body_mass_g, fill = sex) +
  geom_boxplot()

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

dat %>%
  filter(!is.na(sex)) %>%
  ggplot() +
  aes(x = sex, y = body_mass_g, fill = species) +
  geom_boxplot()

Но для более удобочитаемого графика я предпочитаю помещать в качестве цвета переменную с наименьшим количеством уровней (которая на самом деле является аргументом fill в слое aes()) и переменную с наибольшим количеством категорий на оси x ( т. е. аргумент x в слое aes()).

Из средних значений и диаграмм по подгруппам мы уже можем видеть, что в нашем образце:

• самки пингвинов, как правило, имеют меньшую массу тела, чем самцы, и это характерно для всех рассматриваемых видов, и
• масса тела папуасских пингвинов выше, чем у двух других видов.

Имейте в виду, что эти выводы справедливы только в рамках нашей выборки! Чтобы обобщить эти выводы на популяцию, нам нужно выполнить двусторонний дисперсионный анализ и проверить значимость объясняющих переменных. Это цель следующего раздела.

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

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

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

Таким образом, мы начинаем с модели, которая включает два основных эффекта (т. е. пол и вид) и взаимодействие:

# Two-way ANOVA with interaction
# save model
mod <- aov(body_mass_g ~ sex * species,
  data = dat
)
# print results
summary(mod)
##              Df    Sum Sq  Mean Sq F value   Pr(>F)    
## sex           1  38878897 38878897 406.145  < 2e-16 ***
## species       2 143401584 71700792 749.016  < 2e-16 ***
## sex:species   2   1676557   838278   8.757 0.000197 ***
## Residuals   327  31302628    95727                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 11 observations deleted due to missingness

Сумма квадратов (столбец Sum Sq) показывает, что виды объясняют большую часть изменчивости массы тела. Это самый важный фактор в объяснении этой изменчивости.

P-значения отображаются в последнем столбце выходных данных выше (Pr(>F)). Из этих p-значений мы делаем вывод, что на уровне значимости 5%:

• с учетом вида масса тела значительно различается между двумя полами,
• с учетом пола, масса тела значительно различается по крайней мере для одного вида, и
• взаимодействие между полом и видом (отображаемое в строке sex:species в выходных данных выше) имеет большое значение.

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

Если бы, наоборот, взаимодействие не было значимым (т. е. если значение p ≥ 0,05), мы бы убрали этот эффект взаимодействия из модели. В иллюстративных целях ниже приведен код для двустороннего дисперсионного анализа без взаимодействия, называемого аддитивной моделью:

# Two-way ANOVA without interaction
aov(body_mass_g ~ sex + species,
  data = dat
)

Для читателей, которые привыкли выполнять линейные регрессии в R, вы заметите, что структура кода для двустороннего дисперсионного анализа на самом деле аналогична:

• формула dependent variable ~ independent variables
• знак + используется для включения независимых переменных без взаимодействия4
• знак * используется для включения независимых переменных с взаимодействием

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

Обратите внимание, что следующий код также работает и дает те же результаты:

# method 2
mod2 <- lm(body_mass_g ~ sex * species,
  data = dat
)
Anova(mod2)
## Anova Table (Type II tests)
## 
## Response: body_mass_g
##                Sum Sq  Df F value    Pr(>F)    
## sex          37090262   1 387.460 < 2.2e-16 ***
## species     143401584   2 749.016 < 2.2e-16 ***
## sex:species   1676557   2   8.757 0.0001973 ***
## Residuals    31302628 327                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

• ANOVA типа II, когда нет отсутствия значительного взаимодействия, что можно сделать в R с Anova(mod, type = "II"), 5 и
• ANOVA типа III, когда есть значительное взаимодействие, которое можно сделать в R с Anova(mod, type = "III").

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

Попарные сравнения

Учитывая значимость двух основных эффектов, мы пришли к выводу, что:

• с учетом вида, масса тела у самок и самцов различается, и
• с учетом пола масса тела различается по крайней мере для одного вида.

Если масса тела различается между двумя полами, учитывая, что существует ровно два пола, это должно быть потому, что масса тела значительно различается между женщинами и мужчинами.

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

Однако не все так просто для вида. Позвольте мне объяснить, почему это не так просто, как для полов.

Существует три вида (Adelie, Chinstrap и Gentoo), поэтому существует 3 пары видов:

1. Адели и Чинстрап
2. Адели и Генту
3. Подбородочный ремень и Gentoo

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

• масса тела значительно отличается у Адели и антарктического ремешка, но существенно не отличается между Адели и генту и незначительно отличается между антаркальным ремешком и генту, или
• масса тела существенно различается между Адели и Генту, но существенно не отличается между Адели и антарктическим ремешком и незначительно отличается между антарктическим ремешком и генту, или
• масса тела существенно различается у антарктических и генту, но существенно не отличается между Адели и антарктических и не имеет существенных различий между Адели и генту.

Или может быть и так:

• масса тела значительно различается у Адели и антарктического ремешка, а также между Адели и генту, но существенно не отличается между антарктическим ремешком и генту, или
• масса тела значительно различается у Адели и антарктического ремешка, а также между антарктического ремешка и генту, но существенно не различается между антарктически и антарктически, или
• масса тела значительно различается у антарктических и генту, а также между Адели и генту, но существенно не различается между антарктических и антарктических.

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

Что касается одностороннего ANOVA, то на данном этапе мы не можем точно знать, какой вид отличается от какого по массе тела. Чтобы узнать это, нам нужно сравнить каждый вид по два с помощью апостериорных тестов (также известных как попарные сравнения).

Существует несколько апостериорных тестов, наиболее распространенным из которых является HSD Тьюки, который проверяет все возможные пары групп. Как упоминалось ранее, этот тест необходимо проводить только для переменной вида, поскольку для пола существует только два уровня.

Что касается однофакторного дисперсионного анализа, HSD Тьюки можно выполнить в R следующим образом:

# method 1
TukeyHSD(mod,
  which = "species"
)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = body_mass_g ~ sex * species, data = dat)
## 
## $species
##                        diff       lwr       upr     p adj
## Chinstrap-Adelie   26.92385  -80.0258  133.8735 0.8241288
## Gentoo-Adelie    1377.65816 1287.6926 1467.6237 0.0000000
## Gentoo-Chinstrap 1350.73431 1239.9964 1461.4722 0.0000000

или используя пакет {multcomp}:

# method 2
library(multcomp)
summary(glht(
  aov(body_mass_g ~ sex + species,
    data = dat
  ),
  linfct = mcp(species = "Tukey")
))
## 
## 	 Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = body_mass_g ~ sex + species, data = dat)
## 
## Linear Hypotheses:
##                         Estimate Std. Error t value Pr(>|t|)    
## Chinstrap - Adelie == 0    26.92      46.48   0.579     0.83    
## Gentoo - Adelie == 0     1377.86      39.10  35.236   <1e-05 ***
## Gentoo - Chinstrap == 0  1350.93      48.13  28.067   <1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

или используя функцию pairwise.t.test(), используя метод корректировки значения p по вашему выбору: 6

# method 3
pairwise.t.test(dat$body_mass_g, dat$species,
  p.adjust.method = "BH"
)
## 
## 	Pairwise comparisons using t tests with pooled SD 
## 
## data:  dat$body_mass_g and dat$species 
## 
##           Adelie Chinstrap
## Chinstrap 0.63   -        
## Gentoo    <2e-16 <2e-16   
## 
## P value adjustment method: BH

Обратите внимание, что при использовании второго метода именно модель без взаимодействия должна быть указана в функции glht(), даже если взаимодействие является значительным. Более того, не забудьте заменить mod и species в моем коде на имя вашей модели и имя вашей независимой переменной.

Оба метода дают одинаковые результаты, а именно:

• масса тела незначительно различается между подбородочным ремнем и Адели (скорректированное значение p = 0,83),
• масса тела значительно различается между Gentoo и Adelie (скорректированное значение p ‹ 0,001), и
• масса тела значительно различается между Gentoo и Chinstrap (скорректированное значение p ‹ 0,001).

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

Если вы хотите сравнить все комбинации групп, это можно сделать с помощью функции TukeyHSD() и указать взаимодействие в аргументе which:

# all combinations of sex and species
TukeyHSD(mod,
  which = "sex:species"
)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = body_mass_g ~ sex * species, data = dat)
## 
## $`sex:species`
##                                      diff       lwr       upr     p adj
## male:Adelie-female:Adelie        674.6575  527.8486  821.4664 0.0000000
## female:Chinstrap-female:Adelie   158.3703  -25.7874  342.5279 0.1376213
## male:Chinstrap-female:Adelie     570.1350  385.9773  754.2926 0.0000000
## female:Gentoo-female:Adelie     1310.9058 1154.8934 1466.9181 0.0000000
## male:Gentoo-female:Adelie       2116.0004 1962.1408 2269.8601 0.0000000
## female:Chinstrap-male:Adelie    -516.2873 -700.4449 -332.1296 0.0000000
## male:Chinstrap-male:Adelie      -104.5226 -288.6802   79.6351 0.5812048
## female:Gentoo-male:Adelie        636.2482  480.2359  792.2606 0.0000000
## male:Gentoo-male:Adelie         1441.3429 1287.4832 1595.2026 0.0000000
## male:Chinstrap-female:Chinstrap  411.7647  196.6479  626.8815 0.0000012
## female:Gentoo-female:Chinstrap  1152.5355  960.9603 1344.1107 0.0000000
## male:Gentoo-female:Chinstrap    1957.6302 1767.8040 2147.4564 0.0000000
## female:Gentoo-male:Chinstrap     740.7708  549.1956  932.3460 0.0000000
## male:Gentoo-male:Chinstrap      1545.8655 1356.0392 1735.6917 0.0000000
## male:Gentoo-female:Gentoo        805.0947  642.4300  967.7594 0.0000000

Или с помощью функции HSD.test() из пакета {agricolae}, которая обозначает подгруппы, существенно не отличающиеся друг от друга одной буквой:

library(agricolae)
HSD.test(mod,
  trt = c("sex", "species"),
  console = TRUE # print results
)
## 
## Study: mod ~ c("sex", "species")
## 
## HSD Test for body_mass_g 
## 
## Mean Square Error:  95726.69 
## 
## sex:species,  means
## 
##                  body_mass_g      std  r  Min  Max
## female:Adelie       3368.836 269.3801 73 2850 3900
## female:Chinstrap    3527.206 285.3339 34 2700 4150
## female:Gentoo       4679.741 281.5783 58 3950 5200
## male:Adelie         4043.493 346.8116 73 3325 4775
## male:Chinstrap      3938.971 362.1376 34 3250 4800
## male:Gentoo         5484.836 313.1586 61 4750 6300
## 
## Alpha: 0.05 ; DF Error: 327 
## Critical Value of Studentized Range: 4.054126 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##                  body_mass_g groups
## male:Gentoo         5484.836      a
## female:Gentoo       4679.741      b
## male:Adelie         4043.493      c
## male:Chinstrap      3938.971      c
## female:Chinstrap    3527.206      d
## female:Adelie       3368.836      d

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

# set axis margins so labels do not get cut off
par(mar = c(4.1, 13.5, 4.1, 2.1))
# create confidence interval for each comparison
plot(TukeyHSD(mod, which = "sex:species"),
  las = 2 # rotate x-axis ticks
)

Из результатов и графика выше мы делаем вывод, что все комбинации пола и вида значительно различаются, за исключением между самкой антарктического ремешка и самкой Адели (значение p = 0,138) и самцом антарктического ремешка и самцом Адели (значение p = 0,581).

Эти результаты, которые, кстати, соответствуют диаграммам, показанным выше, и которые будут подтверждены визуализациями ниже, завершают двухсторонний ANOVA в R.

Визуализации

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

Во-первых, со средним значением и стандартной ошибкой среднего по подгруппе с помощью функции allEffects() из пакета {effects}:

# method 1
library(effects)
plot(allEffects(mod))

Или с помощью пакета {ggpubr}:

# method 2
library(ggpubr)
ggline(subset(dat, !is.na(sex)), # remove NA level for sex
  x = "species",
  y = "body_mass_g",
  color = "sex",
  add = c("mean_se") # add mean and standard error
) +
  labs(y = "Mean of body mass (g)")

В качестве альтернативы, используя {Rmisc} и {ggplot2}:

library(Rmisc)
# compute mean and standard error of the mean by subgroup
summary_stat <- summarySE(dat,
  measurevar = "body_mass_g",
  groupvars = c("species", "sex")
)
# plot mean and standard error of the mean
ggplot(
  subset(summary_stat, !is.na(sex)), # remove NA level for sex
  aes(x = species, y = body_mass_g, colour = sex)
) +
  geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), # add error bars
    width = 0.1 # width of error bars
  ) +
  geom_point() +
  labs(y = "Mean of body mass (g)")

Во-вторых, если вы предпочитаете рисовать только среднее значение по подгруппе:

with(
  dat,
  interaction.plot(species, sex, body_mass_g)
)

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

# plot mean and standard error of the mean as barplots
ggplot(
  subset(summary_stat, !is.na(sex)), # remove NA level for sex
  aes(x = species, y = body_mass_g, fill = sex)
) +
  geom_bar(position = position_dodge(), stat = "identity") +
  geom_errorbar(aes(ymin = body_mass_g - se, ymax = body_mass_g + se), # add error bars
    width = 0.25, # width of error bars
    position = position_dodge(.9)
  ) +
  labs(y = "Mean of body mass (g)")

Заключение

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

Все это было проиллюстрировано с помощью набора данных penguins, доступного в пакете {palmerpenguins}.

Спасибо за прочтение.

Я надеюсь, что эта статья поможет вам провести двусторонний ANOVA с вашими данными.

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

1. Теоретически однофакторный дисперсионный анализ также можно использовать для сравнения 2 групп, а не только 3 и более. Тем не менее, на практике часто бывает так, что критерий Стьюдента используется для сравнения двух групп, а однофакторный дисперсионный анализ — для сравнения трех или более групп. Выводы, полученные с помощью t-критерия Стьюдента для независимых выборок и однофакторного дисперсионного анализа с 2 группами, будут схожими. ↩︎
2. Обратите внимание, что если предположение о нормальности не выполняется, для его улучшения можно применить множество преобразований, наиболее распространенным из которых является логарифмическое преобразование (функция log() в R).↩︎
3. Обратите внимание, что тест Бартлетта также подходит для проверки предположения о равных дисперсиях. ↩︎
4. В аддитивной модели предполагается, что две объясняющие переменные независимы; они не взаимодействуют друг с другом. ↩︎
5. Где mod — название сохраненной модели. ↩︎
6. Здесь мы используем поправку Benjamini & Hochberg (1995), но вы можете выбрать один из нескольких методов. Подробнее см. ?p.adjust. ↩︎

Статьи по Теме

• АНОВА в R
• Тест Стьюдента в R и вручную: как сравнить две группы при разных сценариях?
• Одновыборочный критерий Уилкоксона в R
• Тест Краскела-Уоллиса, или непараметрический вариант дисперсионного анализа
• Проверка гипотез вручную

Первоначально опубликовано на https://statsandr.com 19 июня 2023 г.