Статья 1. Широкими мазками
Sozontov
21 07 2020
Загрузка данных
library(tidyverse)
library(vegan)
library(foreach)
library(doSNOW)
long <- readxl::read_xlsx("data4.xlsx", sheet = 1)
labs <- readxl::read_xlsx("data4.xlsx", sheet = 7)
wide <- long %>%
select(id,spec, freqlog) %>%
pivot_wider(names_from = spec, values_from = "freqlog") %>%
left_join(., select(labs, id, variable, n.species, n.bones)) %>%
filter(variable == "numeric", n.species > 1, n.bones > 4) %>% # Фильтр
select(-variable, -n.species, -n.bones)Матрица сходства
data.perm <- wide %>%
left_join(., select(labs, id, fs, foreststeppe, region, period2), by = "id") %>%
filter(fs != "islands", period2 != "mix") #filters!
dist.perm <- vegan::vegdist(data.perm[,2:12], method="bray", binary=FALSE)
dist.ord <- vegan::vegdist(wide[,-1], method="bray", binary=FALSE)
m <- as.matrix(dist.ord)
colnames(m) <- paste0("id_", wide$id)
m <- m %>%
as.data.frame %>%
mutate_all(function(a){return(1-round(a, 3))}) %>%
cbind(ID = paste0("id_", wide$id), .)
formattable::as.datatable(formattable::formattable(m))Матрица сходства пусть будет.
Permanova
vegan::adonis(dist.perm ~ fs * period2, permutations = 999, data = data.perm)##
## Call:
## vegan::adonis(formula = dist.perm ~ fs * period2, data = data.perm, permutations = 999)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## fs 6 18.721 3.12018 36.281 0.42311 0.001 ***
## period2 2 1.200 0.59977 6.974 0.02711 0.001 ***
## fs:period2 11 2.051 0.18649 2.168 0.04636 0.001 ***
## Residuals 259 22.274 0.08600 0.50342
## Total 278 44.246 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ординация
pcoa <- ape::pcoa(dist.ord)
eig <- data.frame(eig = pcoa$values$Eigenvalues) %>%
filter(eig > 0) %>%
mutate(eig = round(eig/sum(eig)*100, 1)) %>%
pull(eig)
pcoa <- cbind(id = wide$id, pcoa$vectors[,1:2]) %>%
as.data.frame %>%
left_join(., select(labs, id, fs, period2), by = "id")
ggplot(filter(pcoa, period2 != "mix"), aes(x = Axis.1, y = Axis.2, color = period2)) +
geom_point() +
stat_ellipse(level = 0.8) +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "PCoA, bray-curtis", x = paste0("Axis 1 (", eig[1], "%)"), y = paste0("Axis 2 (", eig[2], "%)"))
ggplot(filter(pcoa, period2 != "mix", fs != "islands"), aes(x = Axis.1, y = Axis.2, color = period2)) +
geom_point() +
stat_ellipse(level = 0.8) +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "PCoA, bray-curtis", x = paste0("Axis 1 (", eig[1], "%)"), y = paste0("Axis 2 (", eig[2], "%)"),
subtitle = "Какую-то одну лесостепь убрать бы, чтоб лучше смотрелось") +
facet_wrap(~fs)
ggplot(filter(pcoa, fs != "islands"), aes(x = Axis.1, y = Axis.2, color = fs)) +
geom_point() +
stat_ellipse(level = 0.8) +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "PCoA, bray-curtis", x = paste0("Axis 1 (", eig[1], "%)"), y = paste0("Axis 2 (", eig[2], "%)"))
ggplot(filter(pcoa, period2 != "mix", fs != "islands"), aes(x = Axis.1, y = Axis.2, color = fs)) +
geom_point() +
stat_ellipse(level = 0.8) +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "PCoA, bray-curtis", x = paste0("Axis 1 (", eig[1], "%)"), y = paste0("Axis 2 (", eig[2], "%)")) + facet_wrap(~period2) +
coord_fixed(ratio = 0.6)
Рецентные 1
rec.nolog <- rbind(
cbind(select(data.perm, id, Благородный, Кабан, Косуля, Лось, fs, period2), year = NA, region = NA, zone = NA),
readxl::read_xlsx("data4.xlsx", sheet = "recent") #load additional data from excel
) %>% mutate(total = Благородный + Кабан + Косуля + Лось) %>%
pivot_longer(names_to = "spec", values_to = "val",
-c("id", "fs", "period2", "year", "region", "zone", "total")) %>%
mutate(val = val / total, period2 = case_when(
period2 == "AT" ~ "1.AT",
period2 == "SB" ~ "2.SB",
period2 == "SA" ~ "3.SA",
period2 == "recent" ~ "4.Rec"))
myplot <- function(data, sp) {
data %>%
filter(spec == sp, val > 0, val != 1) %>%
ggplot(aes(x = period2, y = val)) +
geom_boxplot() +
theme_bw() +
labs(x = "", y = "", title = sp)
}myplot(rec.nolog, sp = "Благородный")
myplot(rec.nolog, sp = "Благородный") + facet_wrap(~fs)
myplot(rec.nolog, sp = "Кабан")
myplot(rec.nolog, sp = "Кабан") + facet_wrap(~fs)
myplot(rec.nolog, sp = "Косуля")
myplot(rec.nolog, sp = "Косуля") + facet_wrap(~fs)
myplot(rec.nolog, sp = "Лось")
myplot(rec.nolog, sp = "Лось") + facet_wrap(~fs)
Рецентные 2
myplot2 <- function(data, reg) {
data %>%
filter(fs == reg, val > 0, val != 1) %>%
ggplot(aes(x = period2, y = val)) +
geom_boxplot() +
theme_bw() +
labs(x = "", y = "", title = reg) +
facet_wrap(~spec)
}
myplot2(rec.nolog, "1.Pannonskaia") + coord_fixed(ratio = 2)
myplot2(rec.nolog, "2.Dneprovskaia")
myplot2(rec.nolog, "3.Srednerusskaia")
myplot2(rec.nolog, "4.Zavolzhskaia")
myplot2(rec.nolog, "5.Tobolskaia")
myplot2(rec.nolog, "6.Ishimskaia")
myplot2(rec.nolog, "7.Barab.Ob")
Рецентные 3
rec.nolog <- rec.nolog %>% pivot_wider(values_from = val, names_from = spec)
rec.nolog$nsp <- apply(rec.nolog[,8:11], 1, FUN = function(a) {length(a[a>0])})
rec.nolog <- filter(rec.nolog, nsp > 1)
dist.ord <- vegdist(rec.nolog[,8:11], method = "bray", binary = FALSE)
pcoa <- ape::pcoa(dist.ord)
eig <- data.frame(eig = pcoa$values$Eigenvalues) %>%
filter(eig > 0) %>%
mutate(eig = round(eig/sum(eig)*100, 1)) %>%
pull(eig)
pcoa <- cbind(rec.nolog[,1:3], pcoa$vectors[,1:2])
ggplot(pcoa, aes(x = Axis.1, y = Axis.2, color = period2, shape = period2)) +
geom_point() +
stat_ellipse(level = 0.8) +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "PCoA, bray-curtis", x = paste0("Axis 1 (", eig[1], "%)"), y = paste0("Axis 2 (", eig[2], "%)"))
ggplot(pcoa, aes(x = Axis.1, y = Axis.2, color = period2, shape = period2)) +
geom_point() +
stat_ellipse(level = 0.8) +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "PCoA, bray-curtis", x = paste0("Axis 1 (", eig[1], "%)"), y = paste0("Axis 2 (", eig[2], "%)"),
subtitle = "хорошо бы разбить на 2 части") + facet_wrap(~fs)
ggplot(pcoa, aes(x = Axis.1, y = Axis.2, color = fs, shape = fs)) +
geom_point() +
stat_ellipse(level = 0.8) +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "PCoA, bray-curtis", x = paste0("Axis 1 (", eig[1], "%)"), y = paste0("Axis 2 (", eig[2], "%)")) + facet_wrap(~period2)
Indicator value
indval <- function(X, class, method = "APCF", perm = 999, setseed = 1) { #APCF, ACF, ASF, PCF, PSF
library(dplyr); library(tidyr)
set.seed(setseed)
if(method == "APCF") {
APCF <- function(X, class) {
indval <- cbind(X, class) %>%
pivot_longer(names_to = "spec", values_to = "x", -class) %>%
mutate(y = case_when(x > 0 ~ 1, x == 0 ~ 0)) %>%
group_by(class, spec) %>%
summarise(x.mean = mean(x), sumy = sum(y), leny = length(y), .groups = 'drop') %>%
group_by(spec) %>%
mutate(A = x.mean/sum(x.mean), B = sumy / leny, indval = A*B*100) %>%
ungroup() %>%
select(class, spec, A, B, indval)
return(indval)
}
res <- APCF(X, class)
res <- cbind(res, matrix(NA, ncol = perm, nrow = nrow(res)))
for(i in 1:perm) {
res[,5+i] <- APCF(X, sample(class, length(class), replace = F))$indval
}
s <- apply(res[,5:ncol(res)], 1, FUN = function(a){
a1 <- a[1]
b <- a[2:length(a)]
return(length(b[b>a1])/length(b))
})
res <- res %>% select(class, spec, A, B, indval) %>%
mutate(indval.s = s, method = method, n.perm = perm)
}
if(method == "ACF") {
ACF <- function(X, class) { #
indval <-cbind(X, class) %>%
pivot_longer(names_to = "spec", values_to = "x", -class) %>%
group_by(class, spec) %>%
mutate(xih = sum(x), nj = length(x), a = case_when(sum(x) == 0 ~ 0,
sum(x) != 0 ~ abs(x/mean(xih) - 1/mean(nj))
)) %>%
summarise(x.mean = mean(x), b = 1 - 0.5*sum(a), .groups = 'drop') %>%
group_by(spec) %>%
mutate(A = x.mean/sum(x.mean), B = b, indval = A*B*100, indval.s = 0) %>%
ungroup() %>%
select(class, spec, A, B, indval, indval.s)
return(indval)
}
res <- ACF(X, class) #
res <- cbind(res, matrix(NA, ncol = perm, nrow = nrow(res)))
for(i in 1:perm) {
res[,5+i] <- ACF(X, sample(class, length(class), replace = F))$indval #
}
s <- apply(res[,5:ncol(res)], 1, FUN = function(a){
a1 <- a[1]
b <- a[2:length(a)]
return(length(b[b>a1])/length(b))
})
res <- res %>% select(class, spec, A, B, indval) %>%
mutate(indval.s = s, method = method, n.perm = perm)
}
if(method == "ASF") {res <- "there is stil not implemented"}
if(method == "PCF ") {
PCF <- function(X, class) {
indval <- cbind(X, class) %>%
pivot_longer(names_to = "spec", values_to = "x", -class) %>%
mutate(y = case_when(x > 0 ~ 1, x == 0 ~ 0)) %>%
group_by(spec) %>%
mutate(ally = sum(y)) %>%
ungroup %>%
group_by(class, spec) %>%
summarise(A = sum(y)/mean(ally), B = sum(y)/length(y), indval = A*B*100, indval.s = 0, .groups = 'drop')
return(indval)
}
res <- PCF (X, class) #
res <- cbind(res, matrix(NA, ncol = perm, nrow = nrow(res)))
for(i in 1:perm) {
res[,5+i] <- PCF (X, sample(class, length(class), replace = F))$indval #
}
s <- apply(res[,5:ncol(res)], 1, FUN = function(a){
a1 <- a[1]
b <- a[2:length(a)]
return(length(b[b>a1])/length(b))
})
res <- res %>% select(class, spec, A, B, indval) %>%
mutate(indval.s = s, method = method, n.perm = perm)
}
if(method == "PSF") {
PSF <- function(X, class) {
indval <- cbind(X, class) %>%
pivot_longer(names_to = "spec", values_to = "x", -class) %>%
mutate(y = case_when(x > 0 ~ 1, x == 0 ~ 0)) %>%
group_by(class, spec) %>%
summarise(ymean = mean(y), B = sum(y)/length(y), .groups = "drop") %>%
group_by(spec) %>%
mutate(A = (ymean - (sum(ymean)-ymean)/(length(ymean) - 1)) / max(ymean), indval = A*B*100, indval.s = 0) %>%
#there can be negative indvals (A, particularly), and it's okay
ungroup() %>%
select(class, spec, A, B, indval, indval.s)
return(indval)
}
res <- PSF(X, class)
res <- cbind(res, matrix(NA, ncol = perm, nrow = nrow(res)))
for(i in 1:perm) {
res[,5+i] <- PSF(X, sample(class, length(class), replace = F))$indval
}
s <- apply(res[,5:ncol(res)], 1, FUN = function(a){
a1 <- a[1]
b <- a[2:length(a)]
return(length(b[b>a1])/length(b))
})
res <- res %>% select(class, spec, A, B, indval) %>%
mutate(indval.s = s, method = method, n.perm = perm)
}
return(as_tibble(res))
}
data.iv <- left_join(wide, select(labs, id, fs, region, period1, n.species, n.bones), by = "id") %>%
rename(period = period1) %>%
left_join(., readxl::read_excel("hier.clust.xlsx", sheet = "data"), by = "period") %>%
filter(n.species > 1, n.bones > 4, n.bones < 1000)cl <- makeCluster(13)
registerDoSNOW(cl)
res1 <- foreach(i=18:24, .combine='rbind') %dopar% {
indval(X = data.iv[,2:12], class = data.iv[[i]], method = "APCF", perm = 999) #
}
res1 <- res1 %>% filter(!is.na(class) & class != "NA")
res1a <- res1 %>% mutate(astr = case_when(indval.s <= 0.001 ~ "***", indval.s <= 0.01 ~ "** ",
indval.s <= 0.05 ~ "* ", indval.s > 0.05 ~ " "),
ind2 = paste0(round(indval, 1), astr, "(", round(indval.s, 3), ")")) %>%
select(class, spec, ind2) %>%
pivot_wider(names_from = `class`, values_from = `ind2`, values_fn = function(a){as.character(a[[1]])})
registerDoSNOW(cl)
res2 <- foreach(i=18:24, .combine='rbind') %dopar% {
indval(X = data.iv[,2:12], class = data.iv[[i]], method = "ACF", perm = 999) #
}
res2 <- res2 %>% filter(!is.na(class) & class != "NA")
res2a <- res2 %>% mutate(astr = case_when(indval.s <= 0.001 ~ "***", indval.s <= 0.01 ~ "** ",
indval.s <= 0.05 ~ "* ", indval.s > 0.05 ~ " "),
ind2 = paste0(round(indval, 1), astr, "(", round(indval.s, 3), ")")) %>%
select(class, spec, ind2) %>%
pivot_wider(names_from = `class`, values_from = `ind2`, values_fn = function(a){as.character(a[[1]])})
data.iv$period <- "all"
registerDoSNOW(cl)
res3 <- foreach(i=13:15, .combine='rbind') %dopar% {
indval(X = data.iv[,2:12], class = data.iv[[i]], method = "APCF", perm = 999) #
}
res3 <- res3 %>% filter(!is.na(class) & class != "NA")
res3a <- res3 %>% mutate(astr = case_when(indval.s <= 0.001 ~ "***", indval.s <= 0.01 ~ "** ",
indval.s <= 0.05 ~ "* ", indval.s > 0.05 ~ " "),
ind2 = paste0(round(indval, 1), astr, "(", round(indval.s, 3), ")")) %>%
select(class, spec, ind2) %>%
pivot_wider(names_from = `class`, values_from = `ind2`, values_fn = function(a){as.character(a[[1]])})
registerDoSNOW(cl)
res4 <- foreach(i=13:15, .combine='rbind') %dopar% {
indval(X = data.iv[,2:12], class = data.iv[[i]], method = "ACF", perm = 999) #
}
res4 <- res4 %>% filter(!is.na(class) & class != "NA")
res4a <- res4 %>% mutate(astr = case_when(indval.s <= 0.001 ~ "***", indval.s <= 0.01 ~ "** ",
indval.s <= 0.05 ~ "* ", indval.s > 0.05 ~ " "),
ind2 = paste0(round(indval, 1), astr, "(", round(indval.s, 3), ")")) %>%
select(class, spec, ind2) %>%
pivot_wider(names_from = `class`, values_from = `ind2`, values_fn = function(a){as.character(a[[1]])})
writexl::write_xlsx( list(
APCF.time = res1a, APCF.time.long = res1,
ACF.time = res2a, ACF.time.long = res2,
APCF.spatial = res3a, APCF.spatial.long = res3,
ACF.spatial = res4a, ACF.spatial.long = res4),
path = "Indicator.Value (22.07.20).xlsx")