Quality, v 2
А. Созонтов
04 01 2022
Загрузка данных
Убраны такие параметры как Sпней укр с пнями. Квм и Лапник. Квм (неполные данные)
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(tidyverse)
library(vegan)
env <- readxl::read_xlsx("Quality.xlsx")
anim <- env %>%
select(14:17, 2:13) %>%
group_by(year, zone , km , line) %>%
summarise_all(sum) %>%
ungroup()
env <- env %>%
select(14:17, 19:ncol(.)) %>%
group_by(year, zone , km , line) %>%
summarise_all(mean) %>%
ungroup()Визуализация мультиколлинеарности
tmp <- expand_grid(x1 = 5:ncol(env),
x2 = 5:ncol(env)) #%>% filter(x1 > x2)
env2 <- list()
for(i in 1:nrow(tmp)) {
env2[[i]] <- cbind(env[,tmp$x1[i]], env[,tmp$x2[i]])
}
env2 <- env2 %>% lapply(function(a){b <- cor.test(a[,1], a[,2])
paste(c(b$estimate, b$p.value,
colnames(a)[1], colnames(a)[2]),
collapse = "__")
}
) %>% map_chr(c)
env2 <- strsplit(x = env2, "__") %>%
lapply(., function(a){matrix(a, ncol = 4)}) %>%
do.call(rbind, .)
env2 <- tibble(est = as.numeric(env2[,1]),
pval = as.numeric(env2[,2]),
x1 = env2[,3], x2 = env2[,4])
env2 %>%
mutate(est = case_when(x1 != x2 ~ est),
pval2 = case_when(pval < 0.001 & !is.na(est) ~ paste0(round(est, 2), "***"),
pval < 0.01 & !is.na(est) ~ paste0(round(est, 2), "**"),
pval < 0.05 & !is.na(est) ~ paste0(round(est, 2), "*"))) %>%
ggplot(., aes(x=x1, y=x2, fill = est)) +
geom_tile(color = "black") +
geom_text(aes(x = x1, y = x2, label = pval2), size = 2.5,
na.rm = TRUE, angle = 30)+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_bw()+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
coord_fixed() +
labs(x = NULL, y = NULL)
Обилие и кол-во видов
div <- tibble(anim[,1:4],
abu = apply(anim[,5:16], 1, sum) ,
n.spec = apply(anim[,5:16], 1, function(a){length(a[a>0])}),
Clglar = anim$Cl_glareolus,
Clruti = anim$Cl_rutilus
)
div## # A tibble: 17 x 8
## year zone km line abu n.spec Clglar Clruti
## <dbl> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 2005 B 6 8 94 7 47 10
## 2 2005 B 6 9 86 6 46 13
## 3 2005 B 6 10 22 4 4 5
## 4 2005 F 20 5 198 10 144 11
## 5 2005 F 20 6 172 9 127 22
## 6 2005 F 20 7 97 8 43 4
## 7 2005 Imp 2 12 24 6 6 8
## 8 2005 Imp 2 13 17 3 7 8
## 9 2005 Imp 2 14 55 4 4 8
## 10 2021 B 6 8 102 5 47 24
## 11 2021 B 6 9 95 7 47 22
## 12 2021 F 20 5 122 9 100 6
## 13 2021 F 20 6 127 6 108 1
## 14 2021 F 20 7 105 9 80 6
## 15 2021 Imp 2 12 34 4 12 6
## 16 2021 Imp 2 13 17 4 4 7
## 17 2021 Imp 2 14 16 5 3 4Среда
dis <- env[,5:18] %>%
scale %>%
vegan::vegdist(method = "euclidean")
pcoa <- ape::pcoa(dis)
eig <- round(pcoa$values$Eigenvalues/sum(pcoa$values$Eigenvalues)*100, 1)
pcoa <- tibble(env[,1:4],
as.data.frame(pcoa$vectors[,1:2]))
pcoa %>% mutate(year =as.character(year),
km = as.character(km)) %>%
ggplot( aes(x = Axis.1, y = Axis.2, color = year, shape = zone)) +
geom_point(size = 3) +
# facet_wrap(~year) +
# stat_ellipse() +
labs(x = paste0("Axis 1, ", eig[1], "%"),
y = paste0("Axis 2, ", eig[2], "%")) +
theme_bw()
pcoa## # A tibble: 17 x 6
## year zone km line Axis.1 Axis.2
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2005 B 6 8 -1.82 1.06
## 2 2005 B 6 9 -1.97 1.36
## 3 2005 B 6 10 0.952 4.01
## 4 2005 F 20 5 -3.08 -2.03
## 5 2005 F 20 6 -2.93 -0.579
## 6 2005 F 20 7 -2.40 -0.342
## 7 2005 Imp 2 12 -0.0147 1.53
## 8 2005 Imp 2 13 -1.39 2.43
## 9 2005 Imp 2 14 -0.669 3.03
## 10 2021 B 6 8 0.0972 -0.511
## 11 2021 B 6 9 3.02 -0.212
## 12 2021 F 20 5 -0.973 -3.43
## 13 2021 F 20 6 0.885 -2.30
## 14 2021 F 20 7 -0.305 -2.85
## 15 2021 Imp 2 12 3.26 -0.584
## 16 2021 Imp 2 13 3.57 -0.371
## 17 2021 Imp 2 14 3.76 -0.214Модели
Обилие общее
div %>%
left_join(pcoa, by = c("year", "zone", "km", "line")) %>%
lm(abu ~ Axis.1*Axis.2, data = .) %>%
summary##
## Call:
## lm(formula = abu ~ Axis.1 * Axis.2, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.324 -25.251 -1.335 17.416 53.959
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 81.353 7.307 11.134 5.1e-08 ***
## Axis.1 -13.603 3.501 -3.886 0.00188 **
## Axis.2 -14.875 3.888 -3.826 0.00210 **
## Axis.1:Axis.2 3.810 3.250 1.172 0.26212
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.13 on 13 degrees of freedom
## Multiple R-squared: 0.7615, Adjusted R-squared: 0.7064
## F-statistic: 13.83 on 3 and 13 DF, p-value: 0.0002436Количество видов
div %>%
left_join(pcoa, by = c("year", "zone", "km", "line")) %>%
lm(n.spec ~ Axis.1*Axis.2, data = .) %>%
summary##
## Call:
## lm(formula = n.spec ~ Axis.1 * Axis.2, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4898 -0.6070 -0.1571 0.6901 2.0758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.2353 0.2525 24.692 2.62e-12 ***
## Axis.1 -0.4217 0.1210 -3.486 0.00402 **
## Axis.2 -0.6032 0.1344 -4.489 0.00061 ***
## Axis.1:Axis.2 0.2612 0.1123 2.325 0.03687 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.041 on 13 degrees of freedom
## Multiple R-squared: 0.8122, Adjusted R-squared: 0.7689
## F-statistic: 18.75 on 3 and 13 DF, p-value: 5.278e-05Обилие рыжей п.
div %>%
left_join(pcoa, by = c("year", "zone", "km", "line")) %>%
lm(Clglar ~ Axis.1*Axis.2, data = .) %>%
summary##
## Call:
## lm(formula = Clglar ~ Axis.1 * Axis.2, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.371 -10.332 -6.472 10.087 41.278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.765 5.568 8.758 8.19e-07 ***
## Axis.1 -9.865 2.668 -3.698 0.002681 **
## Axis.2 -14.687 2.963 -4.957 0.000262 ***
## Axis.1:Axis.2 3.510 2.477 1.417 0.179964
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.96 on 13 degrees of freedom
## Multiple R-squared: 0.8063, Adjusted R-squared: 0.7617
## F-statistic: 18.04 on 3 and 13 DF, p-value: 6.435e-05Обилие красной п.
div %>%
left_join(pcoa, by = c("year", "zone", "km", "line")) %>%
lm(Clruti ~ Axis.1*Axis.2, data = .) %>%
summary##
## Call:
## lm(formula = Clruti ~ Axis.1 * Axis.2, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.0641 -3.7947 -1.9970 -0.2383 14.3784
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.70588 1.80245 5.385 0.000124 ***
## Axis.1 -0.49934 0.86360 -0.578 0.573006
## Axis.2 0.06792 0.95919 0.071 0.944623
## Axis.1:Axis.2 0.02117 0.80179 0.026 0.979334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.432 on 13 degrees of freedom
## Multiple R-squared: 0.02912, Adjusted R-squared: -0.1949
## F-statistic: 0.13 on 3 and 13 DF, p-value: 0.9406Связь новых осей и исходных параметров
full_join(env, pcoa, by = c("year", "zone", "km", "line")) %>%
select(-year, -zone, -km, -line) %>%
cor() %>%
as.data.frame() %>%
rownames_to_column("v1") %>%
pivot_longer(names_to = "v2", values_to = "val", -v1) %>%
filter(v1 %in% c("Axis.1", "Axis.2"),
!(v2 %in% c("Axis.1", "Axis.2"))) %>%
ggplot(aes(x = v1, y = v2, fill = val)) +
geom_tile() +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_bw() +
labs(x = NULL, y = NULL)