• Еще раз про множественные сравнения
• Пример: Поведение песчанок в тесте открытое поле
• perMANOVA
• Условия применимости perMANOVA
• Более подробная интерпретация результатов perMANOVA
• Более сложные дизайны в perMANOVA

Вы сможете

• Вспомнить основные дизайны для тестирования гипотез в рамках дисперсионного анализа
• Применить функции R для тестирования гипотез с помощью perMANOVA

В чем опасность множественных сравнений?

$\alpha$ — это вероятность совершить ошибку первого рода при тестировании гипотезы (= вероятность отвергнуть истинную нулевую гипотезу, = вероятность найти различия там, где их нет).

Обычно принимается, что $H_0$ отвергают на уровне значимости $\alpha = 0.05$.

Когда у нас два средних — все просто, сравнение всего одно.

Естественно, вероятность совершить ошибку I рода для группы сравнений $\alpha _{familywise}$ равна уровню значимости для единственного сравнения $\alpha _{per\,comparison}$.

$\alpha _{familywise} = \alpha _{per\,comparison}$

Но если сравнений много, то растет вероятность совершить хотя бы одну ошибку I рода (найти различия там, где их нет).

Если сравнений много...

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

Пусть мы решили, что в каждом из сравнений будем использовать уровень значимости $0.05$.

Тогда в каждом из сравнений вероятность совершить ошибку первого рода будет $\alpha_{per\,comparison} = 0.05$.

Если сделать $N$ тестов, то вероятность совершить хотя бы одну ошибку I рода в группе тестов (family-wise error rate, FWER) значительно возрастает.

$$\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N$$

Чем больше сравнений, тем больше вероятность обнаружить различия там, где их на самом деле нет.

$$\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N$$

В таблице даны значения $\alpha _{familywise}$ для разного числа сравнений, если $\alpha _{per\,comparison} = 0.05$:

Для решения проблемы есть два подхода

1. Взять более жесткий порог уровня значимости. Например,
   • поправка Бонферрони,
   • поправка Хольма-Бонферрони,
   • процедура Бенъямини-Хохберга.
2. Изменить схему тестирования гипотезы — тестировать не три независимых гипотезы, а одну сложную (так это, например, происходит в ANOVA).

Поправка Бонферрони

$\alpha^*_{per\,comparison} = \alpha_{familywise} /{N}$

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

Ниже $\alpha _{per\,comparison}$ после поправки Бонферрони, сохраняющие $\alpha _{familywise} = 0.05$:

Метод Хольма-Бонферрони

Метод Хольма-Бонферрони — это пошаговая процедура.

Чтобы зафиксировать $FWER \le \alpha_{familywise}$:

1. Сортируем в порядке возрастания $N$ значений $p$, полученные в тестах, и присваиваем каждой ранг $j$ от 1 до $N$:
   $$p_{1} \le p_{2} \le \cdots \le p_{N - 1} \le p_{N}$$
2. Вводим поправку для значений $p$
   $${p^*_{j}} = min{\{(N - j + 1) \cdot p_{j},\;1\}}$$
3. Сравниваем с $\alpha$. Отвергаем $H_0$, если $p^*_{j} < \alpha_{familywise}$

Пример поправки Хольма-Бонферрони

В таблице ниже даны результаты нескольких сравнений ($N = 5$). С помощью поправки Хольма-Бонферрони для каждого $p_j$ мы получим свое скорректированное значение $p^*_{j}$.

Классический дисперсионный анализ

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

$SS_{total} = SS_{between} + SS_{within}$

Тестирование гипотезы о влиянии фактора

Чтобы сравнить межгрупповую изменчивость ( $SS_{between}$) и внутригрупповую ( $SS_{within}$), Фишер предложил статистику

$$F = \frac{SS_{between} / (a - 1)}{SS_{within} / (N - a)}$$ где $a$ — число групп.

F-распределение

Если межгрупповая изменчивость равна внутригрупповой ( $H_0$), то $F$ принадлежит $F$-распределению с двумя параметрами $df_{between} = a - 1$ и $df_{within} = N - a$, где $a$ — число классов, $N$ — общее количество объектов в анализе.

Однофакторный дизайн

Выявляется влияние фактора А.

Многофакторный ортогональный дизайн

Выявляется влияние фактора А, B и их взаимодействия A $\times$ B.

Иерархический дизайн

Некоторые дискретные предикторы могут быть иерархически соподчинены.

Выявляется влияние вложенного фактора B внутри градаций фактора A.

Влияние вложенного фактора B можно оценить, но чаще всего оно не представляет интереса для исследователя.

Рандомизированный полный блочный дизайн

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

Тест “открытое поле”

Задание

• Какую меру различия можно использовать с этими данными?
• Как можно преобразовать данные?
• Постройте ординацию объектов в осях MDS и раскрасьте точки в соответствии с видами

## Данные наблюдений
pesch <- read.csv("data/pesch.csv", header = TRUE, sep = ";")
head(pesch)

Решение

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

options(digits = 4)
log_pesch <- pesch
log_pesch[, 3:ncol(pesch)] <- log(pesch[, 3:ncol(pesch)] + 1)
head(log_pesch)

##   Gender Species Time_to_entrance Urination Defecation Quadr_Number
## 1      f    karl            0.000    0.6931      0.000        3.871
## 2      f    karl            3.045    0.0000      1.386        5.762
## 3      m    karl            5.204    0.0000      0.000        5.182
## 4      f    karl            0.000    0.0000      0.000        3.497
## 5      f    karl            4.942    0.0000      0.000        5.328
## 6      m    karl            0.000    0.0000      1.099        3.664
##   Vert_Number Displ_Act Time_in_Centre
## 1       2.485     1.609          0.000
## 2       4.078     1.386          1.946
## 3       4.025     1.946          1.099
## 4       1.792     1.386          0.000
## 5       3.401     2.398          1.386
## 6       2.303     2.197          0.000

Решение

library(ggplot2)
library(vegan)
theme_set(theme_bw(base_size = 14))
mds_pesch <- metaMDS(log_pesch[, 3:ncol(pesch)], distance = "euclidean")
mds_pesch <- as.data.frame(mds_pesch$points)
mds_pesch$Species <- pesch$Species
ggplot(mds_pesch, aes(x = MDS1, y = MDS2, colour = Species)) + 
  geom_point(size = 5)

Различаются ли виды песчанок по поведению?

Каким способом можно ответить на этот вопрос?

• Мы могли бы взять каждый из поведенческих признаков отдельно и провести одномерный однофакторный дисперсионный анализ.

• Но нас интересует поведение в целом — нужен многомерный анализ.

ANOVA разработан для одномерных данных

Что делать если мы хотим оценивать объект по многим признакам сразу?

Примеры:

• Сообщество как целое (много видов)
• Поведение как целое (много элементов)
• Метаболом как целое (много метаболитов)
• Ответы респондентов на множество вопросов в анкетах

Варианты решений:

• MANOVA (Fisher, 1925, Wilks, 1932)
• perMANOVA (Anderson, 2001; McArdle, Anderson, 2001)
• distance-based Redundancy Analysis (db-RDA) (Legendre, Anderson, 1999)

Многомерный дисперсионный анализ (MANOVA)

Давно разработан параметрический метод MANOVA (Multivariate Analysis Of Variance). Он дает возможность проводить анализ аналогичный ANOVA. В основе MANOVA лежат представление о многомерном нормальном распределении и расстояниях между центроидами.

В MANOVA сравниваются:

• отклонения точек от групповых центроидов (аналог $SS_{within}$)
• отклонения групповых центроидов от общего центроида (аналог $SS_{between}$).

Permutational Multivariate Analysis of Variance (perMANOVA)

Теорема

Сумма квадратов отклонений от центроидов равна сумме квадратов взаимных расстояний, деленной на число объектов

Для Евклидовых расстояний эта закономерность была известна давно (например, Kendall, Stuart 1963).

Для других мер различия центроиды найти гораздо сложнее. Например, для коэффициента Брея-Куртиса (не метрика), среднее значение не будет соответствовать центроиду.

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

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

Anderson, 2001

Разложение дисперсии становится очень простым

$SS_{total} = \frac{1}{N}\sum \sum {d_{ij}^2}$
Сумма квадратов взаимных расстояний — это сумма квадратов субдиагональных элементов, деленная на число объектов N.

$SS_{within} = \frac{1}{n}\sum \sum {d_{ij}^2 \cdot \epsilon_{ij}}$
Внутригрупповая сумма квадратов — это сумма всех сумм квадратов расстояний между элементами для каждой группы, деленная на n число объектов в группе.

Тогда межгрупповая сумма квадратов $SS_{between} = SS_{total} - SS_{within}$

Псевдо-F статистика

$$F = \frac{SS_{between} / (a - 1)}{SS_{within}/(N - a)}$$

Для оценки значимости псевдо-F используется пермутационная процедура:

• Случайным образом перетасовываются строки исходной матрицы
• После каждого акта пермутации вычисляется $F_{perm}$
• Уровень значимости
• Внимание! Для иерархического дизайна процедура пермутации имеет свои особенности (обсудим позднее).

$$p = \frac {N_{F_{perm} \ge F}} {N_{permutations}}$$

Применим метод perMANOVA (функция adonis())

library(vegan)
permanova_pesch <- adonis(log_pesch[3:9] ~ Species, data = log_pesch,
                          method = "euclidean")
permanova_pesch

## 
## Call:
## adonis(formula = log_pesch[3:9] ~ Species, data = log_pesch,      method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)  
## Species    2        60   30.02    3.58 0.237  0.012 *
## Residuals 23       193    8.39         0.763         
## Total     25       253                 1.000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Мы видим традиционную для ANOVA таблицу результатов. Что здесь что?

Результаты perMANOVA

## 
## Call:
## adonis(formula = log_pesch[3:9] ~ Species, data = log_pesch,      method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)  
## Species    2        60   30.02    3.58 0.237  0.012 *
## Residuals 23       193    8.39         0.763         
## Total     25       253                 1.000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Поведение разных видов песчанок в тесте открытое поле значимо различалось (perMANOVA, p = 0.012).

Условия применимости perMANOVA

1. Требуется равенство внутригрупповых дисперсий (гомоскедастичность).

2. Желательно использование сбалансированных данных (с равным объемом групп), т.к. этом случае perMANOVA устойчив к гетерогенности дисперсии (Anderson, Walsh, 2013).

Проверка равенства внутригрупповых дисперсий

Для проверки можно использовать функцию betadisper(), изначально предназначенную для сравнения $\beta$-разнообразия сообществ.

Функция betadisper() вычисляет внутригрупповые центроиды и координаты точек в пространстве главных координат (Principal Coordinate Analysis = PCoA = metric MDS).

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

dist_pesch <- vegdist(log_pesch[,3:ncol(pesch)], method  = "euclidean")
PCO_pesch <- betadisper(dist_pesch, log_pesch$Species)

График ординации PCoA

Объект, возвращаемый betadisper(), позволяет также нарисовать наши объекты в пространстве главных координат (PCoA).

plot(PCO_pesch, main = "PCoA ordination")

Процедура PERMDISP2 для проверки равенства внутригрупповых дисперсий

Процедура PERMDISP2 реализована в пакете vegan в функции anova.betadisper(). Это многомерный аналог теста Левина на гомогенность дисперсий в группах, который иногда используется для проверки условий применимости дисперсионного анализа.

anova(PCO_pesch)

## Analysis of Variance Table
## 
## Response: Distances
##           Df Sum Sq Mean Sq F value Pr(>F)
## Groups     2    9.7    4.86    2.05   0.15
## Residuals 23   54.5    2.37

• Не выявлено значимых различий разброса внутригрупповых расстояний.

Для визуализации можно нарисовать боксплот

boxplot(PCO_pesch)

Post hoc тесты в perMANOVA

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

Для выявления попарных различий нужны попарные сравнения.

Функция для попарных сравнений perMANOVA

pairwise_permanova <- function(dat, group, strata = NULL, ...){
  pair <- combn(unique(as.character(group)), 2)
  ncomb <- ncol(pair)
  res <- rep(NA, ncomb)
  for (i in 1:ncomb) {
    filter <- group %in% pair[, i]
    if(is.null(strata)){
      posthoc <- adonis(dat[filter, ] ~ group[filter], ...)$aov.tab$Pr[1]
    } else {
      posthoc <- adonis(dat[filter, ] ~ group[filter], 
                        strata = strata[filter], ...)$aov.tab$Pr[1]
    }
    res[i] <- posthoc
    names(res)[i] <- paste(pair[, i], collapse = " vs. ")
  }
  return(res)
}

Результаты попарных сравнений

p_vals <- pairwise_permanova(
  dat = log_pesch[, -c(1:2)], group = log_pesch$Species, 
  method = "euclidean", permutations=9999)
p_vals

##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. zhirnokhvost 
##                  0.0011                  0.4029                  0.0031

Это все? Пишем статью?

Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

Поправка на множественные сравнения

Мы делали три пары сравнений — нужно ввести поправку для уровня значимости. Будем считать значимыми различия в тех сравнениях, где после введения поправки $p < 0.05$.

Можно посчитать скорректированные значения уровня значимости $p$ с учетом поправки Хольма-Бонферрони.

p.adjust(p_vals, method = "holm")

##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. zhirnokhvost 
##                  0.0033                  0.4029                  0.0062

Многофакторный дизайн в perMANOVA

Выясним, влияет ли пол и вид песчанок на поведение.

Отфильтруем исходные данные (в случае с жирнохвостыми песчанками были изучены только самки)

log_pesch2 <- log_pesch[log_pesch$Species != "zhirnokhvost", ]

Проведем двухфакторный анализ perMANOVA

Различается ли поведение песчанок в зависимости от видовой принадлежности и пола?

twofact_pesch <- adonis(log_pesch2[,3:ncol(pesch)] ~ Gender * Species,
                        data = log_pesch2, method = "euclidean")
twofact_pesch

## 
## Call:
## adonis(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender * Species,      data = log_pesch2, method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)   
## Gender          1       7.2     7.2    1.13 0.049  0.292   
## Species         1      45.5    45.5    7.16 0.311  0.004 **
## Gender:Species  1       4.5     4.5    0.72 0.031  0.529   
## Residuals      14      88.9     6.4         0.608          
## Total          17     146.1                 1.000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Здесь возможен иерархический дизайн

Различается ли поведение самцов и самок у этих видов песчанок?

nested_pesch <- adonis(log_pesch2[, 3:ncol(pesch)] ~ Gender, 
                       data = log_pesch2, strata = log_pesch2$Species, 
                       method = "euclidean")
nested_pesch

## 
## Call:
## adonis(formula = log_pesch2[, 3:ncol(pesch)] ~ Gender, data = log_pesch2,      method = "euclidean", strata = log_pesch2$Species) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)
## Gender     1       7.2    7.21    0.83 0.049    0.3
## Residuals 16     138.9    8.68         0.951       
## Total     17     146.1                 1.000

Внимание! Пермутации производятся только в пределах группирующего фактора (аргумент strata)

Задание

• Создайте датафрейм из файла simulated_data.csv (Это данные симулированные по алгоритму, приведенному в справке по функции adonis())
• В этом датафрейме записано обилие двух видов на экспериментальных площадках двух типов: без добавления и с добавлением NO3, по 6 повторностей в каждом эксперименте.
  Эксперименты были независимо проведены на 3 полях.
• Оцените, зависит ли структура совместного поселения этих двух видов от концентрации NO3.

Решение

com <- read.csv("data/simulated_data.csv", sep = ',', header = T)
# Ошибочный дизай
com_permanova <- adonis(com[,1:2] ~ com$NO3)
# Правильный дизайн
com_permanova2 <- adonis(com[,1:2] ~ com$NO3, strata = com$field)

Summary

• При одновременном тестировании нескольких гипотез растет вероятность ошибки I рода.
• Чтобы контролировать вероятность ошибки I рода, либо используют более жесткий уровень значимости, либо тестируют сложную гипотезу вместо нескольких простых.
• perMANOVA дает возможность тестировать сложные гипотезы в отношении явлений, описанных по многим переменным (т.е. на многомерных данных).
• В perMANOVA можно использовать любые коэффициенты различия.
• Для применения perMANOVA требуется равенство разбросов точек между центроидами их групп, но при равных объемах групп анализ устойчив к отклонениям от этого условия.
• При использовании perMANOVA важно не запутаться в дизайне.

Другие программы

• Primer 6.0 + PERMANOVA Коммерческий продукт.
• PAST Некоммерческая программа. Здесь метод называется NPMANOVA.
• Оригинальная программа М. Андерсон PERMANOVA и PERMDISP.

Что почитать

• Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26: 32–46.
• Anderson, M.J. 2005. PERMANOVA: a FORTRAN computer program for permutational multivariate analysis of variance. Department of Statistics, University of Auckland, New Zealand.
• Anderson, M.J. (2004). PERMDISP: a FORTRAN computer program for permutational analysis of multivariate dispersions (for any two-factor ANOVA design) using permutation tests. Department of Statistics, University of Auckland, New Zealand.
• Legendre P., Legendre L. (2012) Numerical ecology. Second english edition. Elsevier, Amsterdam.