- Протестировать гипотезу о различиях между дискретными группами многомерных данных с помощью теста ANOSIM
- Выявить переменные, вносящие наибольший вклад в формирование различий между группами, применив процедуру SIMPER
Вы сможете
Вспомним основы
Этапы работы с гипотезами (Underwood, 1997)
Формулировка биологической гипотезы
Численное выражение биологической гипотезы (\(H\))
Формулировка альтернативной гипотезы (\(H_0\) - нулевой гипотезы)
Тестирование нулевой гипотезы
Если \(H_0\) ложна, то признаем, что верна \(H\) и, следовательно, биологическую гипотезу можно рассматривать как истинное утверждение
Если \(H_0\) истинна, то…
ANOSIM: Analysis Of Similarity
Пример: Динамика сообществ мидиевых банок
Существуют ли различия в сообществах мидиевых банок, сформированных старыми и молодыми мидиями (Khaitov, 2013)?
com <- read.csv("data/mussel_beds.csv", sep=';', header = T)
com$Mussel_size <- factor(com$Mussel_size)
ascam <- read.csv("data/ASCAM.csv", sep=';', header = T)
com — усредненные данные по обилию 12 видов для 3 мидиевых банок (Vor2, Vor4, Vor5).ascam — (averaged size class abundance matrix) средние плотности поселения мидий разных размеров (6] размерных классов)
Задание
- Постройте ординацию всех описаний датасета
com(прологарифмированные данные) в осях nMDS на основе матрицы Брея-Куртиса - Раскрасьте разными цветами точки, относящиеся к двум разным группам: “Large-dominated” и “Small-dominated”
Решение
library(vegan) library(ggplot2) log_com <- log(com[,-c(1:3)] + 1) ord_log_com <- metaMDS(log_com, distance = "bray", k=2, autotransform = F, trace = 0) MDS <- data.frame(ord_log_com$points)
Решение
Обратите внимание, здесь есть две группы расстояний между точками
ggplot(MDS, aes(x = MDS1, y = MDS2, fill = com$Mussel_size)) +
geom_point(shape = 21, size = 4) +
scale_fill_manual(values = c("red", "blue")) +
labs(fill = "Mussel size structure type") + ggtitle(paste("Stress = ", round(ord_log_com$stress, 3), sep = " "))
Расстояния между объектами
1. Внутригрупповые расстояния
Within group distance
2. Межгрупповые расстояния
between group distance
Ранги расстояний
Для работы удобно (но не обязательно!) перейти от исходных расстояний между объектами, к их рангам.
Обозначим внутригрупповые расстояния (ранги), как \(r_w\), а межгрупповые, как \(r_b\).
Вычислим
- средние значения внутригрупповых рангов расстояний \(R_w\)
- средние значения межгрупповых рангов расстояний \(R_b\)
R - статистика
На основе полученных значений можно построить статистику (Clarke, 1988, 1993)
\[R_{global} = \frac{R_b-R_w}{n(n-1)/4}\]
Эта статистика распределена в интервале -1 < \(R_{global}\) < 1
\(R_{global}\) - оценивает степень разделенности групп в пространстве признаков.
\(R_{global} > 0\) – как минимум есть тенденция к разделению групп
\(R_{global} \rightarrow 1\) – группы разделяются хорошо
\(R_{global} \leq 0\) – нет разделения на группы
ANOSIM вручную
Ранги расстояний в пространстве признаков
Задание:
- Вычислите матрицу коэффициентов Брея-Куртиса на основе матрицы
log_com - Разверните полученную матрицу в вектор.
- На основе полученного вектора создайте вектор, содержащий ранги расстояний.
Решение
dist_com <- vegdist(log_com, method = "bray") rank_dist_com <- rank(dist_com)
Внутри- и межгрупповые расстояния
Задание:
- Создайте треугольную матрицу
dummy_dist, той же размерности, что и матрицаdist_com, в которой0будет с стоять в тех ячейках, которые соответствуют межгрупповым расстояниями, а1– внутригрупповым.
Решение
dummy_dist <- dist(as.numeric(factor(com$Mussel_size))) dummy_dist <- ifelse(dummy_dist == 0, 0, 1)
R-статистка
Задание:
- Вычислите средние значения рангов внутри- и межгрупповых расстояний
- Вычислите R-статистику
Решение
dists <- data.frame(rank_dist = rank_dist_com, dummy = as.vector(dummy_dist)) library(dplyr) mean_dists <- dists %>% group_by(dummy) %>% summarize(rank_mean = mean(rank_dist)) n <- nrow(log_com) R_glob <- (mean_dists$rank_mean[2] - mean_dists$rank_mean[1])/(n * (n - 1)/4)
Вопрос: Каков дальнейший ход процедуры
ANOSIM?
Пермутации
Задание:
- Напишите пользовательскую функцию (пусть она будет называться
R_perm), которая пермутировала бы принадлежность каждого объекта к той или иной группе и вычисляла бы значение R-статистики для новой комбинации. - Используя функцию
for()вычислите 10000 значений R-статистики и запишите их в вектор.
Hint. Весь код для пользовательской функции почти такой же, как написали ранее.
Решение
# функция. Весь код почти такой же, как написали ранее.
R_perm <- function(comm, group){
require(vegan)
dist_com <- vegdist(comm)
rank_dist_com <- rank(dist_com)
dummy_dist <- dist(sample(as.numeric(group))) #Перемешиваем группы
dummy_dist <- ifelse(dummy_dist == 0, 0, 1)
dists <- data.frame(rank = rank_dist_com, dummy = as.vector(dummy_dist))
require(dplyr)
mean_dists <- dists %>% group_by(dummy) %>% summarize(rank_type = mean(rank))
n <- nrow(log_com)
R_p <- (mean_dists$rank_type[2] - mean_dists$rank_type[1])/(n * (n - 1)/4)
R_p
}
R_perms <- rep(NA, 10000)
for(i in 1:9999) R_perms[i] <- R_perm(comm = log_com, group = com$Mussel_size)
Вопрос: Что надо добавить в вектор R_perms?
Распределение пермутационных оценок R-статистики
Задание:
- Постройте частотную гистограмму, характеризующую распределение пермутационных оценок.
- Нанесите на нее полученное значение \(R_{global}\).
- Вычислите уровень значимости.
Решение
R_perms[10000] <- R_glob Pl_manual <- ggplot(data.frame(R_perms), aes(x = R_perms)) + geom_histogram(binwidth = 0.01) + geom_vline(xintercept = R_glob, linetype = 2) + xlim(-0.4, 0.4) Pl_manual
Решение
Вычисляем уровень значимости
p_level <- mean(R_perms >= R_glob) p_level
## [1] 0.0009
Процедура ANOSIM в пакете vegan
Специализированная функция anosim()
com_anosim <- anosim(log_com,
grouping = com$Mussel_size,
permutations = 9999,
distance = "bray")
Задание
Изучите структуру объекта com_anosim и постройте частотное распределение значений \(R_{global}\), полученных при каждом акте пермутации
Решение
anosim_perm <- data.frame(perm = com_anosim$perm) anosim_perm[(com_anosim$permutations + 1), 1] <- com_anosim$statistic Pl_prof <- ggplot(anosim_perm, aes(x = perm)) + geom_histogram(binwidth = 0.01, color = "black", fill = "blue") + geom_vline(xintercept = com_anosim$statistic, linetype = 2) + xlim(-0.4, 0.4) Pl_prof
Сравним результаты ручной и автоматической обработки
Результаты процедуры ANOSIM
summary(com_anosim)
## ## Call: ## anosim(x = log_com, grouping = com$Mussel_size, permutations = 9999, distance = "bray") ## Dissimilarity: bray ## ## ANOSIM statistic R: 0.219 ## Significance: 0.001 ## ## Permutation: free ## Number of permutations: 9999 ## ## Upper quantiles of permutations (null model): ## 90% 95% 97.5% 99% ## 0.0567 0.0791 0.1062 0.1363 ## ## Dissimilarity ranks between and within classes: ## 0% 25% 50% 75% 100% N ## Between 1 302 549 766 946 459 ## Large_dominated 4 193 360 645 945 351 ## Small_dominated 3 275 478 683 938 136
Ограничения (Assumptions) для применения ANOSIM
1) Внутригрупповые расстояния (ранги) для разных групп должны иметь приблизительно равные медианы и пределы варьирования.
Для проверки этого допущения надо сравнить распределения внутригрупповых и межгрупповых расстояний (рангов)
Распределение расстояний имеет следующий вид
plot(com_anosim)
ANOSIM позволяет сравнивать одновременно и несколько групп
НО! Есть одно очень важное ограничение:
2) Попарные сравненения групп можно осуществлять только если было показано, что \(R_{global}\) статистически значимо.
Если это условие выполнено, то можно проводить попарные сравнения
Пример
Пусть у нас есть три группы объектов: A, B, C.
Можно вычислить \(R_{A vs B}\), \(R_{A vs C}\), \(R_{B vs C}\).
Но при больших объемах выборки даже незначительные различия будут достоверны. Важно обращать внимание не только на оценку статистической значимости, но и на значения R!
Важно! При множественных сравнениях придется вводить поправку Бонферрони.
NB! Для сравнения нескольких групп многомерных объектов, есть более мощное средство - PERMANOVA
Задание
- Постройте ординацию в осях nMDS, раскрасив точки в разные цвета в зависимости от номера мидиевой банки
- Проверьте гипотезу о различиях в структуре сообществ на разных банках
- Проверьте условия применимости ANOSIM
- Проведите попарное сравнение всех банок
Решение
График ординации
ggplot(MDS, aes(x = MDS1, y = MDS2, fill = com$Bank)) +
geom_point(shape = 21, size = 4) +
scale_fill_manual(values = c("red", "blue", "green")) +
labs(fill = "Mussel beds") + ggtitle(paste("Stress = ", round(ord_log_com$stress, 3), sep = " "))
Решение
Проверка гипотезы о различиях в структуре сообществ на разных банках
bed_anosim <- anosim(log_com, grouping = com$Bank,permutations = 999, distance = "bray") bed_anosim
## ## Call: ## anosim(x = log_com, grouping = com$Bank, permutations = 999, distance = "bray") ## Dissimilarity: bray ## ## ANOSIM statistic R: 0.44 ## Significance: 0.001 ## ## Permutation: free ## Number of permutations: 999
Решение
Условия применимости
plot(bed_anosim)
Решение
Попарные сравнения Vor2 vs Vor4
# Vor2 vs Vor4
anosim(log_com [com$Bank %in% c("Vor2", "Vor4"),],
grouping = com$Bank[com$Bank %in% c("Vor2", "Vor4")])
##
## Call:
## anosim(x = log_com[com$Bank %in% c("Vor2", "Vor4"), ], grouping = com$Bank[com$Bank %in% c("Vor2", "Vor4")])
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.588
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
Попарные сравнения Vor2 vs Vor5
# Vor2 vs Vor5
anosim(log_com[com$Bank %in% c("Vor2", "Vor5"),], grouping = com$Bank[com$Bank %in% c("Vor2", "Vor5")])
##
## Call:
## anosim(x = log_com[com$Bank %in% c("Vor2", "Vor5"), ], grouping = com$Bank[com$Bank %in% c("Vor2", "Vor5")])
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.309
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
Попарные сравнения Vor4 vs Vor5
# Vor4 vs Vor5
anosim(log_com [com$Bank %in% c("Vor4", "Vor5"),], grouping = com$Bank[com$Bank %in% c("Vor4", "Vor5")])
##
## Call:
## anosim(x = log_com[com$Bank %in% c("Vor4", "Vor5"), ], grouping = com$Bank[com$Bank %in% c("Vor4", "Vor5")])
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.426
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
Проблема малых выборок
Мощность ANOSIM невелика.
При малых выборках пермутационная оценка уровня значимости может не выявить достоверных различий, даже при очень высоком значении R.
SIMPER: Similarity Percentages
Какие признаки зависимой матрицы вносят наибольший вклад в формирование различий между группами?
log_com_simper <- simper(log_com, group = com$Mussel_size, permutations = 999) summary (log_com_simper)
## ## Contrast: Large_dominated_Small_dominated ## ## average sd ratio ava avb cumsum p ## Polydora_quadrilobata 0.04185 0.02941 1.42 0.862 2.425 0.177 0.003 ** ## Hydrobia_ulvae 0.02598 0.01812 1.43 3.259 2.826 0.287 0.903 ## Skeneopsis_planorbis 0.02343 0.01722 1.36 0.380 1.282 0.387 0.001 *** ## Fabricia_sabella 0.02234 0.01799 1.24 0.921 1.312 0.481 0.348 ## Cricotopus_vitripennis 0.02097 0.01512 1.39 1.687 2.221 0.570 0.130 ## Onoba_aculeus 0.01879 0.01319 1.42 1.066 1.575 0.650 0.026 * ## Nemertini 0.01718 0.01417 1.21 2.159 2.511 0.722 0.145 ## Macoma_balthica 0.01544 0.01299 1.19 2.605 2.678 0.788 0.654 ## Littorina_saxatilis 0.01542 0.01003 1.54 2.045 1.689 0.853 0.005 ** ## Filamentous_algae 0.01505 0.01145 1.31 0.828 0.629 0.917 0.958 ## Gammarus_sp. 0.01145 0.00892 1.28 1.766 1.916 0.965 0.889 ## Tubificoides_benedeni 0.00821 0.00613 1.34 4.859 4.755 1.000 0.924 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Permutation: free ## Number of permutations: 999
Оценка вклада в формирование различий
\[Contr_i = \frac{|{y_{i,j}-y_{i,k}}|}{\sum{(y_{i,j}-y_{i,k})}}\]
\(y_{i,j}-y_{i,k}\) - разность значений \(i\)-той переменной (вида) в объектах (пробах) \(j\) и \(k\)
Величины \(Contr_i\) далее усредняются для всех межгрупповых пар \(average\)
\(sd\) - среднеквадратичное отклонение
Отношение \(average/sd\) характеризует дискриминирующую силу переменной
Задание
Выявите виды, отвечающие за различия в сообществах разных банок
Решение
summary(simper(log_com, group = com$Bank, permutations = 999 ))
## ## Contrast: Vor2_Vor4 ## ## average sd ratio ava avb cumsum p ## Polydora_quadrilobata 0.04661 0.02795 1.67 0.6680 2.804 0.180 0.002 ** ## Hydrobia_ulvae 0.03246 0.01559 2.08 3.6840 2.456 0.305 0.004 ** ## Skeneopsis_planorbis 0.02820 0.01816 1.55 0.0486 1.460 0.414 0.001 *** ## Fabricia_sabella 0.02229 0.01699 1.31 1.5530 1.095 0.500 0.471 ## Nemertini 0.02103 0.01389 1.51 1.9237 2.531 0.582 0.001 *** ## Filamentous_algae 0.01947 0.01314 1.48 1.1917 0.419 0.657 0.004 ** ## Onoba_aculeus 0.01845 0.01375 1.34 0.8647 1.416 0.728 0.214 ## Littorina_saxatilis 0.01741 0.01013 1.72 2.1832 1.361 0.795 0.001 *** ## Cricotopus_vitripennis 0.01703 0.01182 1.44 2.0702 2.343 0.861 0.986 ## Macoma_balthica 0.01426 0.01224 1.16 2.7410 2.767 0.916 0.831 ## Gammarus_sp. 0.01218 0.00936 1.30 1.8253 1.994 0.963 0.389 ## Tubificoides_benedeni 0.00949 0.00664 1.43 5.0329 4.712 1.000 0.066 . ## ## Contrast: Vor2_Vor5 ## ## average sd ratio ava avb cumsum p ## Hydrobia_ulvae 0.0289 0.02376 1.216 3.6840 3.139 0.133 0.182 ## Fabricia_sabella 0.0275 0.01847 1.491 1.5530 0.531 0.260 0.001 *** ## Polydora_quadrilobata 0.0220 0.02301 0.955 0.6680 0.889 0.361 1.000 ## Cricotopus_vitripennis 0.0217 0.01407 1.541 2.0702 1.222 0.461 0.190 ## Onoba_aculeus 0.0192 0.01373 1.396 0.8647 1.525 0.550 0.112 ## Filamentous_algae 0.0180 0.01220 1.479 1.1917 0.634 0.633 0.026 * ## Macoma_balthica 0.0168 0.01317 1.274 2.7410 2.375 0.710 0.288 ## Nemertini 0.0165 0.01208 1.368 1.9237 2.439 0.786 0.565 ## Skeneopsis_planorbis 0.0138 0.01094 1.261 0.0486 0.673 0.850 0.995 ## Gammarus_sp. 0.0123 0.00893 1.374 1.8253 1.641 0.907 0.335 ## Littorina_saxatilis 0.0102 0.00712 1.430 2.1832 2.196 0.953 1.000 ## Tubificoides_benedeni 0.0101 0.00690 1.461 5.0329 4.703 1.000 0.014 * ## ## Contrast: Vor4_Vor5 ## ## average sd ratio ava avb cumsum p ## Polydora_quadrilobata 0.04702 0.02844 1.65 2.804 0.889 0.196 0.001 *** ## Cricotopus_vitripennis 0.02797 0.01679 1.67 2.343 1.222 0.313 0.001 *** ## Hydrobia_ulvae 0.02465 0.01588 1.55 2.456 3.139 0.416 0.835 ## Skeneopsis_planorbis 0.02372 0.01545 1.54 1.460 0.673 0.514 0.002 ** ## Fabricia_sabella 0.02009 0.01698 1.18 1.095 0.531 0.598 0.856 ## Littorina_saxatilis 0.01874 0.01087 1.72 1.361 2.196 0.676 0.001 *** ## Onoba_aculeus 0.01641 0.01182 1.39 1.416 1.525 0.745 0.807 ## Macoma_balthica 0.01628 0.01375 1.18 2.767 2.375 0.813 0.395 ## Nemertini 0.01431 0.01378 1.04 2.531 2.439 0.872 0.893 ## Filamentous_algae 0.01229 0.00855 1.44 0.419 0.634 0.924 0.998 ## Gammarus_sp. 0.01154 0.00931 1.24 1.994 1.641 0.972 0.624 ## Tubificoides_benedeni 0.00676 0.00494 1.37 4.712 4.703 1.000 0.994 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Permutation: free ## Number of permutations: 999
Summary
- ANOSIM — простейший вариант сравнения нескольких групп объектов, охарактеризованных по многим признакам.
- С помощью процедуры SIMPER можно оценить вклад отдельных переменных в формирование различий между группами.
Что почитать
- Clarke, K. R., Gorley R. N. (2006) PRIMER v6: User Manual/Tutorial. PRIMER-E, Plymouth.
- Legendre P., Legendre L. (2012) Numerical ecology. Second english edition. Elsevier, Amsterdam.