Загрузка

library(vegan)
library(tidyverse)
library(plotly)
theme_set(theme_bw() + theme(legend.position = "bottom"))
S <- list(RD   = c("C_g", "C_r", "A_o", "M_ag", "M_ar", "S_u", "M_m"), 
          IV   = c("S_a","S_c", "S_m", "N_f"), 
          b.lv = c("PF", "BF", "BFM", "FF", "RS", "MM", "D", "W"))
df <-readxl::read_xlsx("Karabash_23.07.2022.xlsx", 
                       sheet = "main",
                       range = "A2:W246") %>% 
    select(-biotop_old) %>% 
    mutate(biotop = factor(biotop, levels = S$b.lv))
# Recompute ---------------------------------------------------------------
df.l <- df %>% 
    select(-method, -effort) %>% # -biotop_exact, 
    pivot_longer(names_to = "sp", values_to = "abu", -c(1:7)) %>% 
    group_by(zone, line, year, biotop, biotop_exact, sp) %>% 
    summarise(abu = sum(abu), days = mean(days), traps = sum(traps), .groups = "drop") %>% 
    mutate(abu = abu/traps/days*100) %>% 
    select(-traps, -days) %>% 
    filter(!(sp %in% c("T_s", "M_n")))
df.w <- df.l %>% pivot_wider(names_from = sp, values_from = abu)

Обозначения

readxl::read_xlsx("Karabash_23.07.2022.xlsx", 
                       sheet = "labs",
                       range = "A2:C10") %>% 
    select(1, 3) %>% 
    formattable::format_table()
названия биотопов abb2
сосновый лес PF
березовый лес BF
березовые редины BFM
пойменный лес FF
болото тростниковое RS
луг MM
свалка бытовых отходов D
техногенная пустошь W

Базовые расчеты

Весь градиент

df %>% 
    select(-zone:-effort) %>% 
    sum %>% 
    cat("Всего особей отловлено:", .)
## Всего особей отловлено: 443
df %>% 
    select(-zone:-effort) %>% 
    colSums() %>% 
    (function(x){length(x[x>0])}) %>% 
    cat("Всего видов обнаружено:", ., " (Вместе с лаской, - одна особь)")
## Всего видов обнаружено: 13  (Вместе с лаской, - одна особь)
df %>% 
    pull(effort) %>% 
    sum %>% 
    cat("Отработано", ., "ловушко-суток")
## Отработано 5466 ловушко-суток
df %>% 
    select(-biotop:-effort) %>% 
    pivot_longer(names_to = "spec", values_to = "abu", -zone) %>% 
    mutate(taxa = case_when(spec %in% S$RD ~ "Rodents", 
                            spec %in% S$IV ~ "Insectivores")) %>% 
    group_by(taxa) %>% 
    summarise(abu = sum(abu)) %>% 
    mutate(`part, %` = round(abu/sum(abu)*100))
## # A tibble: 3 × 3
##   taxa           abu `part, %`
##   <chr>        <dbl>     <dbl>
## 1 Insectivores    91        21
## 2 Rodents        349        79
## 3 <NA>             3         1
df %>% 
    select(-biotop:-effort) %>% 
    pivot_longer(names_to = "spec", values_to = "abu", -zone) %>% 
    mutate(taxa = case_when(spec %in% S$RD ~ "Rodents", 
                            spec %in% S$IV ~ "Insectivores")) %>% 
    group_by(zone, taxa) %>% 
    summarise(abu = sum(abu), .groups = "drop_last") %>% 
    mutate(`part, %` = round(abu/sum(abu)*100)) %>% 
    ungroup() 
## # A tibble: 6 × 4
##   zone  taxa           abu `part, %`
##   <chr> <chr>        <dbl>     <dbl>
## 1 F     Insectivores    67        25
## 2 F     Rodents        203        75
## 3 F     <NA>             2         1
## 4 I     Insectivores    24        14
## 5 I     Rodents        146        85
## 6 I     <NA>             1         1
df %>% 
    transmute(zone, biotop, abu = round(apply(df[,10:ncol(df)], 1, sum)/df$effort*100, 1)) %>% 
    arrange(desc(abu)) %>% 
    group_by(zone) %>% 
    filter(abu == max(abu)) %>% 
    ungroup %>% 
    mutate(where = "Максимальное обилие: ", .before = 1)
## # A tibble: 2 × 4
##   where                   zone  biotop   abu
##   <chr>                   <chr> <fct>  <dbl>
## 1 "Максимальное обилие: " F     D       80  
## 2 "Максимальное обилие: " I     FF      73.3
df %>% 
    select(line, year, 10:ncol(df)) %>% 
    pivot_longer(names_to = "sp", values_to = "abu", -c("line", "year")) %>% 
    group_by(line, year) %>% 
    summarise(abu = sum(abu), .groups = "drop") %>%  
    filter(abu == 0) %>% 
    nrow() %>% 
    cat("`Нулевых` линий всего: ", ., "; Это ", round(./123*100, 1), "%")
## `Нулевых` линий всего:  43 ; Это  35 %

По зонам

cat("Imp Zone:")
## Imp Zone:
df %>% 
    filter(zone == "I") %>% 
    select(-zone:-effort) %>% 
    sum %>% 
    cat("Импактная зона. Всего особей отловлено:", .)
## Импактная зона. Всего особей отловлено: 171
df %>% 
    filter(zone == "I") %>% 
    select(-zone:-effort) %>% 
    colSums() %>% 
    (function(x){length(x[x>0])}) %>% 
    cat("Импактная зона. Всего видов обнаружено:", ., " (Вместе с лаской, - 1 особь)")
## Импактная зона. Всего видов обнаружено: 11  (Вместе с лаской, - 1 особь)
df %>% 
    filter(zone == "I") %>% 
    pull(effort) %>% 
    sum %>% 
    cat("Импактная зона. Отработано", ., "ловушко-суток")
## Импактная зона. Отработано 2925 ловушко-суток
## FON
cat("Fon Zone:")
## Fon Zone:
df %>% 
    filter(zone == "F") %>% 
    select(-zone:-effort) %>% 
    sum %>% 
    cat("Фоновая зона. Всего особей отловлено:", .)
## Фоновая зона. Всего особей отловлено: 272
df %>% 
    filter(zone == "F") %>% 
    select(-zone:-effort) %>% 
    colSums() %>% 
    (function(x){length(x[x>0])}) %>% 
    cat("Фоновая зона. Всего видов обнаружено:", ., "")
## Фоновая зона. Всего видов обнаружено: 12
df %>% 
    filter(zone == "F") %>% 
    pull(effort) %>% 
    sum %>% 
    cat("Фоновая зона. Отработано", ., "ловушко-суток")
## Фоновая зона. Отработано 2541 ловушко-суток

Метод ловушек

df2 <- tibble(
    select(df, zone, biotop, year, method), 
    (df[,10:ncol(df)] / df$effort * 100)
    ) %>% 
    filter(apply(df[,10:ncol(df)], 1, sum) != 0)

dis <- vegan::vegdist(select(df2, -zone:-method), method = "bray", binary = FALSE)

pcoa <- ape::pcoa(dis)
p <- tibble(df2[,1:4], as_tibble(pcoa$vectors[,1:2])) %>% 
    mutate(year = as.character(year)) %>% 
    # unite("type", zone, biotop) %>% 
    ggplot(aes(x = Axis.1, y = Axis.2, 
               shape = year, color = method)) + 
    geom_point() + 
    facet_wrap(vars(zone, biotop))
p

plotly::ggplotly(p)
# визуально разницы нет

vegan::adonis2(dis ~ zone + year + method, data = df2)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = dis ~ zone + year + method, data = df2)
##           Df SumOfSqs      R2      F Pr(>F)    
## zone       1    1.437 0.03052 4.0000  0.001 ***
## year       1    1.564 0.03324 4.3559  0.001 ***
## method     1    1.331 0.02827 3.7053  0.001 ***
## Residual 119   42.738 0.90797                  
## Total    122   47.070 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
vegan::adonis2(dis ~ zone + year + method + biotop, data = df2)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = dis ~ zone + year + method + biotop, data = df2)
##           Df SumOfSqs      R2      F Pr(>F)    
## zone       1    1.437 0.03052 5.2790  0.001 ***
## year       1    1.564 0.03324 5.7487  0.001 ***
## method     1    1.331 0.02827 4.8901  0.001 ***
## biotop     6   11.988 0.25468 7.3419  0.001 ***
## Residual 113   30.751 0.65330                  
## Total    122   47.070 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# стат.значимо, но r2 очень маленький все равно

tibble(df2[,1:4], abu = apply(df2[,5:ncol(df2)], 1, sum)) %>% 
    ggplot(aes(x = method, y = abu, fill = method)) + 
    geom_boxplot(color = "black")  +
    geom_jitter(shape = 21, color = "black", fill = "black", alpha = 0.3, size = 2)

tibble(df2[,1:4], abu = apply(df2[,5:ncol(df2)], 1, sum)) %>% 
    t.test(abu ~ method, data = .)
## 
##  Welch Two Sample t-test
## 
## data:  abu by method
## t = -5.1617, df = 109.58, p-value = 1.101e-06
## alternative hypothesis: true difference in means between group давилки and group живоловки is not equal to 0
## 95 percent confidence interval:
##  -19.400644  -8.635949
## sample estimates:
##   mean in group давилки mean in group живоловки 
##                11.31504                25.33333
# уловистость живоловок значимо выше

df %>% 
    select(-biotop_exact, -effort) %>% 
    pivot_longer(names_to = "sp", values_to = "abu", -c(1:7)) %>% 
    group_by(zone, biotop, line, year, method, sp) %>% 
    summarise(abu = sum(abu), days = mean(days), traps = sum(traps), .groups = "drop") %>% 
    mutate(abu = abu/traps/days*100) %>% 
    select(-traps, -days) %>% 
    mutate(taxa = case_when(sp %in% S$RD ~ "Rodents", sp %in% S$IV ~ "Insectivores")) %>% 
    filter(!is.na(taxa)) %>% 
    group_by(zone, biotop, line, year, method, taxa) %>% 
    summarise(abu = sum(abu), .groups = "drop") %>% 
    ggplot(aes(x = "", y = abu, color = method)) + 
    geom_boxplot() + 
    geom_jitter(alpha = 0.3) + 
    facet_wrap(vars(zone, taxa)) + 
    labs(x = NULL, y = "Обилие, экз. (особей)")

df %>% 
    select(-biotop_exact, -effort) %>% 
    pivot_longer(names_to = "sp", values_to = "abu", -c(1:7)) %>% 
    group_by(zone, biotop, line, year, method, sp) %>% 
    summarise(abu = sum(abu), days = mean(days), traps = sum(traps), .groups = "drop") %>% 
    mutate(abu = abu/traps/days*100) %>% 
    select(-traps, -days) %>% 
    mutate(taxa = case_when(sp %in% S$RD ~ "Rodents", sp %in% S$IV ~ "Insectivores")) %>% 
    filter(!is.na(taxa)) %>% 
    group_by(zone, biotop, line, year, method, taxa) %>% 
    summarise(abu = sum(abu), .groups = "drop") %>% 
    ggplot(aes(x = "", y = abu, color = method)) + 
    geom_boxplot() + 
    geom_jitter(alpha = 0.3) + 
    facet_wrap(vars(zone, taxa), scales = "free") + 
    labs(x = NULL, y = "Обилие, экз./100 лов.-сут.", 
         subtitle = "переделать перевод на 100 лов.сут")

df2 %>% 
    select(-year, -biotop) %>% # -`Arvicola terrestris L.`
    rbind(., .) %>% ##
    mutate(zone = c(df2$zone, rep("zones_all", nrow(df2)))) %>%  ##
    pivot_longer(names_to = "species", values_to = "abu", -c(1:2)) %>% 
    group_by(zone, method, species) %>% 
    summarise(abu = mean(abu), .groups = "drop_last") %>% 
    mutate(abu = abu/max(abu), .groups = "drop") %>% 
    ungroup %>% 
    filter(!(species %in% c("Tamias", "Neomus", "Mustela nivalis"))) %>% 
    ggplot(aes(x = species, y = abu, fill = method)) + 
    geom_col(position = "dodge") + 
    facet_wrap(~zone) + 
    coord_flip()

Разнообразие

div <- list()
div$zone <- df.w %>% 
    select(-line:-biotop_exact) %>% 
    group_by(zone) %>% 
    summarise_all(sum) %>% 
    transmute(zone, n.zone = apply(.[,-1], 1, function(a){length(a[a>0])}))
div$biotop <- df.w %>% 
    select(-line, -year, -biotop_exact) %>% 
    group_by(zone, biotop) %>% 
    summarise_all(sum) %>% 
    ungroup() %>% 
    transmute(zone, biotop, n.zn.biotop = apply(.[,-c(1:2)], 1, function(a){length(a[a>0])}))
div$line <- df.w %>% 
    select(-year, -biotop_exact) %>% 
    group_by(zone, biotop, line) %>% 
    summarise_all(sum) %>% 
    ungroup() %>% 
    transmute(zone, biotop, line, n.zn.btp.line = apply(.[,-c(1:3)], 1, function(a){length(a[a>0])}))
cat("Бета-разнообразие (Бета Уиттекера), макромасштаб:")
## Бета-разнообразие (Бета Уиттекера), макромасштаб:
left_join(div$biotop, div$zone, by = "zone") %>% 
    group_by(zone) %>% 
    summarise(b = mean(n.zone)/mean(n.zn.biotop))
## # A tibble: 2 × 2
##   zone      b
##   <chr> <dbl>
## 1 F      2.14
## 2 I      2.76
cat("Бета-разнообразие (Mean Bray dissimilarity), макромасштаб:")
## Бета-разнообразие (Mean Bray dissimilarity), макромасштаб:
df.w %>%
    select(-line, -year, -biotop_exact) %>% 
    group_by(zone, biotop) %>% 
    summarise_all(mean) %>% 
    ungroup() %>% 
    split(., .$zone) %>% 
    lapply(function(a){select(a, -zone, -biotop)}) %>% 
    lapply(function(a){vegan::vegdist(a)}) %>% 
    sapply(mean)
##         F         I 
## 0.8110515 0.8247600
cat("Бета-разнообразие (Бета Уиттекера), микромасштаб:")
## Бета-разнообразие (Бета Уиттекера), микромасштаб:
left_join(div$line, div$zone, by = "zone") %>% 
    group_by(zone) %>% 
    summarise(b = mean(n.zone)/mean(n.zn.btp.line))
## # A tibble: 2 × 2
##   zone      b
##   <chr> <dbl>
## 1 F      3.07
## 2 I      4.4
cat("Бета-разнообразие (Mean Bray dissimilarity), микромасштаб:")
## Бета-разнообразие (Mean Bray dissimilarity), микромасштаб:
df.w %>%
    select(-year, -biotop_exact) %>% 
    group_by(zone, biotop, line) %>% 
    summarise_all(mean) %>% 
    ungroup() %>%  split(., .$zone) %>% 
    lapply(function(a){select(a, -zone, -biotop, -line)}) %>% 
    lapply(function(a){
        d <- vegan::vegdist(a)
        d[is.nan(d)] <- 1
        mean(d) 
    })
## $F
## [1] 0.8238319
## 
## $I
## [1] 0.8581699
div$beta.whittaker <- div$line %>% 
    left_join(div$biotop, by = c("zone", "biotop")) %>% 
    left_join(div$zone, by = "zone") %>% 
    group_by(zone, biotop) %>% 
    summarise(y = mean(n.zn.biotop)/mean(n.zn.btp.line), .groups = "drop", type = "Whittaker's W")
div$beta.meanbray <- df.w %>%
    select(-year, -biotop_exact) %>% 
    unite("id", zone, biotop) %>% 
    group_by(id, line) %>% 
    summarise_all(mean) %>% 
    ungroup() %>%  
    split(., .$id) %>% 
    lapply(function(a){
        a <- select(a, -id, -line)
        d <- vegan::vegdist(a)
        d[is.na(d) | is.nan(d)] <- 1
        mean(d)
        }) %>% 
    as_vector() %>% 
    tibble(id = names(.), y = ., type = "Mean Bray dissimilarity") %>% 
    separate(id, into = c("zone", "biotop"))

p <- div$line %>% 
    group_by(zone, biotop) %>% 
    summarise(n.zn.biotop = mean(n.zn.btp.line), .groups = "drop") %>% 
    rbind(div$biotop) %>% 
    transmute(zone, biotop, y = n.zn.biotop, 
              type = c(rep("Всреднием видов в линии", 15), rep("Всего видов в биотопе", 15))) %>% 
    rbind(div$beta.whittaker, div$beta.meanbray) %>% 
    filter(biotop != "W") %>% 
    ggplot(aes(x = biotop,
               y = y, fill = zone)) + 
    geom_col(position = "dodge", width = 0.6) + 
    scale_fill_manual(values = hcl(h = seq(15, 375, length = 2 + 1), l = 65, c = 100)[2:1])  + 
    facet_grid(rows = vars(type), scales = "free") +
    labs(x = NULL, y = NULL, subtitle = "Вверху 1: бета-разнообразие в пределах биотопа (среднее растояние Брея-Кёртиса)
Вверху 2: бета-разнообразие в пределах биотопа по Уиттекеру (sample = line)
Внизу 1: альфа-разнообразие - всего видов в биотопе
Внизу 2: альфа-разнообразие - в среднем в линии биотопа")
plotly::ggplotly(p)
cat("Примечание.
Вверху 1: бета-разнообразие в пределах биотопа (среднее растояние Брея-Кёртиса)
Вверху 2: бета-разнообразие в пределах биотопа по Уиттекеру (sample = line)
Внизу 1: альфа-разнообразие - всего видов в биотопе
Внизу 2: альфа-разнообразие - в среднем в линии биотопа")
## Примечание.
## Вверху 1: бета-разнообразие в пределах биотопа (среднее растояние Брея-Кёртиса)
## Вверху 2: бета-разнообразие в пределах биотопа по Уиттекеру (sample = line)
## Внизу 1: альфа-разнообразие - всего видов в биотопе
## Внизу 2: альфа-разнообразие - в среднем в линии биотопа