Отчёт 2022

library(tidyverse)
theme_set(theme_bw() + theme(
      text = element_text(family = "serif", size = 15),
      legend.position = "bottom")
      )
load(paste0(Sys.Date()-2, "_inext.RData"))
rar.yz <- rar.yz %>% map(~mutate(.x, 
                       zone = factor(zone, levels = c("фоновая", "буферная", "импактная")), 
                       zone = fct_relabel(zone, ~paste0(.x, " зона")))
               )

long <- readxl::read_excel("data/Carabidae_12.12.2022.xlsx", 
    sheet = "main_data", range = "A1:H887") %>% 
    mutate(tur = factor(tur), 
           zone = case_when(zone == "fon" ~ "фоновая", 
                                   zone == "bufer" ~ "буферная", 
                                   TRUE ~ "импактная"),
           zone = factor(zone, levels = c("фоновая", "буферная", "импактная")))
wide <- long %>% 
      select(-num) %>% 
      pivot_wider(names_from = taxa, values_from = abu, 
                  values_fn = sum, values_fill = 0) %>% 
      mutate(id = paste0(year, tur, site, plot), .before = year) %>% 
      select(-no_insects)
tr <- readxl::read_excel("data/Carabidae_12.12.2022.xlsx", 
    sheet = "traits") %>% 
    select(taxa, size_1, stratum.ext:wings)

div <- tibble(wide[,2:6], 
    abu = apply(wide[,7:ncol(wide)], 1, sum),
    nsp = apply(wide[,7:ncol(wide)], 1, function(a){length(a[a>0])}),
    shan= vegan::diversity(wide[,6:ncol(wide)], 
                           MARGIN = 1, index = "shannon", base = exp(1)), 
    dom = apply(wide[,7:ncol(wide)], 1, function(a){a <- a/sum(a); return(max(a))})
    ) %>% 
    mutate(km = as.numeric(substr(site, 2,3)), .before = plot) %>% 
      left_join(mutate(rar.yzstp, 
                       year = as.numeric(year), 
                       plot = as.numeric(plot)), 
                by = c("year", "tur", "zone", "site", "plot")) %>% 
      mutate(d50 = case_when(is.na(d50) ~ 0, TRUE ~ d50))

Базовые подсчёты

Виды

Количество видов

## total number of species: 58

Количество видов по зонам

## # A tibble: 3 × 2
##   zone      nsp_total
##   <fct>         <int>
## 1 фоновая          40
## 2 буферная         36
## 3 импактная        19

Количество видов по зонам и годам

## # A tibble: 3 × 3
##   zone      `2009` `2014`
##   <fct>      <int>  <int>
## 1 фоновая       35     30
## 2 буферная      30     27
## 3 импактная     15     13

Особи

Особи по зонам

## # A tibble: 3 × 2
##   zone      n_individuals
##   <fct>             <dbl>
## 1 фоновая            3024
## 2 буферная           2639
## 3 импактная           660

Особи по зонам и годам

## # A tibble: 3 × 3
##   zone      `2009` `2014`
##   <fct>      <dbl>  <dbl>
## 1 фоновая     1240   1784
## 2 буферная    1258   1381
## 3 импактная    357    303

Обилие

Обилие по зонам, годам и турам учётов

## # A tibble: 3 × 5
##   zone      `2009_tur_1` `2009_tur_2` `2014_tur_1` `2014_tur_2`
##   <fct>            <dbl>        <dbl>        <dbl>        <dbl>
## 1 фоновая           526.         42.4        561.         232  
## 2 буферная          401.         21.7        276.         185. 
## 3 импактная         153.          5.9         59.6         75.1

Разрежение и экстраполяция

Точка - наблюдаемое количество особей и видов.
Линия - разрежение и экстраполяция видового богатства.
Полоса - бутстрепный доверительный интервал (95%) для прогноза.

Туры объединены

rar.yz$line %>% 
      mutate(year = as.factor(year)) %>% 
             
             # zone = factor(zone, levels = c("фоновая", "буферная", "импактная")),
             # zone = fct_relabel(zone, ~paste0(.x, " зона"))) %>% 
      ggplot(aes(x = m, y = qD, ymin = qD.LCL, ymax = qD.UCL, color = year, fill = year)) + 
    geom_ribbon(alpha = 0.1, linetype = "blank") +
    geom_line(linewidth = 1, mapping = aes(linetype = Method)) +
    geom_point(data = mutate(rar.yz$point, year = as.factor(year)), 
               shape = 22, size = 2, color = "black") +
    # geom_label(data = mutate(filter(rar.yz$line, m == 2500), year = as.factor(year)),
    #            mapping = aes(x = 2800, y = qD, label = year), color = "black") +
    facet_grid(cols = vars(zone)) + 
    scale_linetype_manual(values = c("dotted", "solid")) +
      # scale_x_continuous(limits = c(0, 2500)) +
    guides(linetype = "none") +
      theme(panel.grid.minor = element_blank()) +
    labs(x = "Особей", y = "Видов", subtitle = "Разрежение/экстраполяция", 
         fill = "Год", color = "Год")

Туры разделены

Первый тур

rar.ytz$line %>% 
      filter(tur == "1") %>% 
      mutate(year = as.factor(year) 
             # zone = paste0(zone, " зона")
             ) %>% 
      ggplot(aes(x = m, y = qD, ymin = qD.LCL, ymax = qD.UCL, color = year, fill = year)) + 
    geom_ribbon(alpha = 0.1, linetype = "blank") +
    geom_line(linewidth = 1, mapping = aes(linetype = Method)) +
    geom_point(data = mutate(filter(rar.ytz$point, tur == "1"), 
                             # zone = paste0(zone, " зона"),
                             year = as.factor(year)),
               shape = 22, size = 2, color = "black") +
    geom_label(data = mutate(filter(rar.ytz$line, m == 1600, tur == "1"), 
                             # zone = paste0(zone, " зона"),
                             year = as.factor(year)),
               mapping = aes(x = 1200, y = qD, label = year), color = "black") +
    facet_grid(cols = vars(zone)) + 
    scale_linetype_manual(values = c("dotted", "solid")) +
      scale_x_continuous(limits = c(0, 1600)) +
    guides(linetype = "none") +
    labs(x = "Особей", y = "Видов", subtitle = "Разрежение/экстраполяция
Тур первый", 
         fill = "Год", color = "Год")

Второй тур

rar.ytz$line %>% 
      filter(tur == "2") %>% 
      mutate(year = as.factor(year) 
             # zone = paste0(zone, " зона")
             ) %>% 
      ggplot(aes(x = m, y = qD, ymin = qD.LCL, ymax = qD.UCL, color = year, fill = year)) + 
    geom_ribbon(alpha = 0.1, linetype = "blank") +
    geom_line(linewidth = 1, mapping = aes(linetype = Method)) +
    geom_point(data = mutate(filter(rar.ytz$point, tur == "2"), 
                             # zone = paste0(zone, " зона"),
                             year = as.factor(year)),
               shape = 22, size = 2, color = "black") +
    geom_label(data = mutate(filter(rar.ytz$line, m == 1600, tur == "2"), 
                             # zone = paste0(zone, " зона"),
                             year = as.factor(year)),
               mapping = aes(x = 1200, y = qD, label = year), color = "black") +
    facet_grid(cols = vars(zone)) + 
    scale_linetype_manual(values = c("dotted", "solid")) +
      scale_x_continuous(limits = c(0, 1600)) +
    guides(linetype = "none") +
    labs(x = "Особей", y = "Видов", subtitle = "Разрежение/экстраполяция
Тур второй", 
         fill = "Год", color = "Год")

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

Обилие

Участки

div %>% 
       mutate(tur = paste0("Тур ", tur), 
              km = as.factor(km)) %>% 
ggplot(aes(x = km, y = abu, fill = zone)) + 
    geom_boxplot() + 
    facet_grid(cols = vars(year)) +# , rows = vars(tur)) + 
    labs(x = NULL, y = "Обилие (особей на 100 лов.-сут.)", 
         title = "Обилие (особей на 100 лов.-сут.)",
         subtitle = "Столбцы - туры учёта, строки - годы") + 
    # theme(axis.text.x = element_text(angle = 90)) + 
    scale_fill_manual(values = c("green", "yellow", "red")) +
    guides(fill="none")

Километры

div %>% 
      mutate(year = factor(year), 
             tur = paste0("Тур ", as.character(tur))) %>% 
    ggplot(aes(x = km, y = abu, color = year, fill = year)) + 
    geom_smooth(formula = "y ~ x", #se = FALSE,
                method = "loess", linewidth = 0.5, linetype = "dashed") +
    geom_point(shape = 21, size = 2) +
    scale_x_continuous(breaks = c(1, 4, 5, 9, 11, 12, 18, 26, 27, 32)) +
      scale_y_log10() +
    facet_grid(rows = vars(tur), scales = "free") + #, rows = vars(zone), ) + 
    theme(axis.text.x = element_text(angle = 45), 
          panel.grid.minor.x = element_blank()) + 
    labs(x = "Расстояние до завода", y = "Обилие (особей на 100 лов.-сут.)", 
         title = "Обилие (особей на 100 лов.-сут.)",
         subtitle = "Столбцы - туры учёта, строки - годы
Пунктирные линии - регрессия линейная и скользящая", 
color = "Год", fill = "Год")

Количество видов

Участки

div %>% 
      mutate(tur = paste0("Тур ", tur), 
             km = as.factor(km)) %>% 
ggplot(aes(x = km, y = nsp, fill = zone)) + 
    geom_boxplot() + 
    facet_grid(cols = vars(year)) + #, rows = vars(tur)) + 
    labs(x = NULL, y = "Количество видов", 
         title = "Количество видов",
         subtitle = "Столбцы - туры учёта, строки - годы") + 
    # theme(axis.text.x = element_text(angle = 90)) + 
    scale_fill_manual(values = c("green", "yellow", "red")) +
    guides(fill="none")

Километры

div %>% 
      mutate(year = factor(year), 
             tur = paste0("Тур ", as.character(tur))) %>% 
    ggplot(aes(x = km, y = nsp, color = year, fill = year)) + 
    geom_smooth(formula = "y ~ x", #se = FALSE,
                method = "loess", linewidth = 0.5, linetype = "dashed") +
    geom_point(shape = 21, size = 2) +
    scale_x_continuous(breaks = c(1, 4, 5, 9, 11, 12, 18, 26, 27, 32)) +
      scale_y_log10() +
    facet_grid(rows = vars(tur), scales = "free") + #, rows = vars(zone), ) + 
    theme(axis.text.x = element_text(angle = 45), 
          panel.grid.minor.x = element_blank())+ 
    labs(x = "Расстояние до завода", y = "Количество видов", 
         title = "Количество видов",
         subtitle = "Столбцы - туры учёта, строки - годы
Пунктирные линии - регрессия линейная и скользящая", 
      fill = "Год", color = "Год") 

Индекс разнообразия Шеннона

Участки

div %>% 
      mutate(tur = paste0("Тур ", tur), 
             km = as.factor(km)) %>% 
      ggplot(aes(x = km, y = shan, fill = zone)) + 
    geom_boxplot() + 
    facet_grid(cols = vars(year)) + # , rows = vars(tur)) + 
    labs(x = NULL, y = "Индекс разнообразия Шеннона", 
         title = "Индекс разнообразия Шеннона",
         subtitle = "Столбцы - туры учёта, строки - годы") + 
    # theme(axis.text.x = element_text(angle = 90)) + 
    scale_fill_manual(values = c("green", "yellow", "red")) +
    guides(fill="none")

Километры

div %>% 
      mutate(year = factor(year), 
             tur = paste0("Тур ", as.character(tur))) %>% 
    ggplot(aes(x = km, y = shan, color = year, fill = year)) + 
    geom_smooth(formula = "y ~ x", #se = FALSE,
                method = "loess", linewidth = 0.5, linetype = "dashed") +
    geom_point(shape = 21, size = 2) +
    scale_x_continuous(breaks = c(1, 4, 5, 9, 11, 12, 18, 26, 27, 32)) +
      scale_y_log10() +
    facet_grid(rows = vars(tur), scales = "free") + #, rows = vars(zone), ) + 
    theme(axis.text.x = element_text(angle = 45), 
          panel.grid.minor.x = element_blank())+
    labs(x = "Расстояние до завода", y = "Индекс разнообразия Шеннона", 
         title = "Индекс разнообразия Шеннона",
         subtitle = "Столбцы - туры учёта, строки - годы
Пунктирные линии - регрессия линейная и скользящая", 
      color = "Год", fill = "Год")

Индекс доминирования Б.-Паркера

Участки

div %>% 
      mutate(tur = paste0("Тур ", tur), 
             km = as.factor(km)) %>% 
ggplot(aes(x = km, y = dom, fill = zone)) + 
    geom_boxplot() + 
    facet_grid(cols = vars(year)) + # , rows = vars(tur)) + 
    labs(x = NULL, y = "Индекс доминирования", 
         title = "Индекс доминирования Бергера-Паркера",
         subtitle = "Столбцы - туры учёта, строки - годы") + 
    # theme(axis.text.x = element_text(angle = 90)) + 
    scale_fill_manual(values = c("green", "yellow", "red")) +
    guides(fill="none")

Километры

div %>% 
      mutate(year = factor(year), 
             tur = paste0("Тур ", as.character(tur))) %>% 
    ggplot(aes(x = km, y = dom, color = year, fill = year)) + 
    geom_smooth(formula = "y ~ x", #se = FALSE,
                method = "loess", linewidth = 0.5, linetype = "dashed") +
    geom_point(shape = 21, size = 2) +
    scale_x_continuous(breaks = c(1, 4, 5, 9, 11, 12, 18, 26, 27, 32)) +
      # scale_y_log10() +
    facet_grid(rows = vars(tur), scales = "free") + #, rows = vars(zone), ) + 
    theme(axis.text.x = element_text(angle = 45), 
          panel.grid.minor.x = element_blank())+
    labs(x = "Расстояние до завода", y = "Индекс доминирования", 
         title = "Индекс доминирования Бергера-Паркера",
         subtitle = "Столбцы - туры учёта, строки - годы
Пунктирные линии - регрессия линейная и скользящая", 
color = "Год", fill = "Год")

Разнообразие: effect size

D <- c("abu", "Динамическая плотность", 
"nsp", "Количество видов",
"shan", "Индекс Шеннона", 
"dom", "Индекс доминирования",
"d50", "Количество видов на 50 особей",
"size1_small", "мелкие виды",
"size1_medium", "средние виды", 
"size1_big", "крупные виды", 
"class_zoophaga", "зоофаги",
"class_myxophaga", "миксофаги", 
"habitats_meadow", "луговые", 
"habitats_forest", "лесные", 
"habitats_forest.meadow", "эвритопные", 
"stratum_strato", "стратобионты",  #"LS",
"stratum_geohorto", "хортобионты",  # "HSD",
"stratum_epigeo", "эпигеобионты",   #"SS",
"stratum_stratohorto", "хортобионты",    #"HSD",
"fenol_all.seasons", "эврихронные", 
"humid_xero", "ксерофилы",
"humid_meso", "мезофилы", 
"humid_hygro", "гигрофилы", 
"wings_macropt", "высокомобильные", 
"wings_brachypter", "немобильные",
"wings_dimorph", "слабомобильные") %>% 
    matrix(ncol = 2, byrow = TRUE) %>% 
    as.data.frame() %>% 
    rename(eng2 = 1, rus = 2) %>% 
    separate(eng2, into = c("eng1", "eng2"), sep = "_")

dens <- function(x, y = 0.05){ 
      sh <- data.frame(
            x = 0, 
            y = y, 
            label = paste0("p.val = ", round(shapiro.test(x)$p.value, 3)))
      ggplot(data.frame(x), aes(x)) + 
            geom_density() +
            geom_label(mapping = aes(x, y, label = label), data = sh) +
            labs(x = "residuals", y = NULL)
      }
join <- function(a){ 
   tr %>% 
      select(taxa, type = {{a}}) %>% 
      left_join(long, ., by = "taxa") %>%
      filter(taxa != "no_insects") %>%       
      left_join(rename(D, type = eng2), by = "type") %>% ##
      select(-type, -eng1) %>% ##
      group_by(year, tur, site, plot, rus) %>% ##
      summarise(abu = sum(abu), .groups = "drop") %>% 
      # mutate(type = paste0(a, "_", rus)) %>% 
      pivot_wider(names_from = rus, values_from = abu, values_fill = 0)
}

div <- div %>% 
      left_join(join("size_1"), by = c("year", "tur", "site", "plot")) %>% 
      left_join(join("class"), by = c("year", "tur", "site", "plot")) %>% 
      left_join(join("habitats"), by = c("year", "tur", "site", "plot")) %>% #
      left_join(join("stratum"), by = c("year", "tur", "site", "plot")) %>% 
      left_join(join("humid"), by = c("year", "tur", "site", "plot")) %>% 
      left_join(join("wings"), by = c("year", "tur", "site", "plot"))

ES <- list()
ES$by_tur <- div %>% 
    # select(abu:wings_dimorph) %>% 
    select(abu:`слабомобильные`) %>% 
    colnames() %>% 
    map(function(V = "abu", df = div){
    df %>% 
        rename(VV = {{V}}) %>% 
        split(paste0(df$year, df$tur)) %>% 
        lapply(FUN = function(a){
            rbind(
            SingleCaseES::LRRi(a$VV[a$zone == "фоновая"], a$VV[a$zone == "буферная"]), 
            SingleCaseES::LRRi(a$VV[a$zone == "фоновая"], a$VV[a$zone == "импактная"])) %>% 
                mutate(zone = c("буферная", "импактная"), type = V)
            }) %>% 
        map_df(rbind, .id = "tur") %>% 
        mutate(year = substr(tur, 1, 4), tur = substr(tur, 5, 5), .before = 1)
}) %>% 
    map_df(rbind) %>% 
    as_tibble()

ES$year_only <- div %>% 
    select(abu:`слабомобильные`) %>% 
    #select(abu:wings_dimorph) %>% 
    colnames() %>% 
    map(function(V, df = div){
      # map(function(V = "abu", df = div){
    df %>% 
        rename(VV = {{V}}) %>% 
        split(paste0(df$year)) %>% 
        lapply(FUN = function(a){
            rbind(
            SingleCaseES::LRRi(a$VV[a$zone == "фоновая"], a$VV[a$zone == "буферная"]), 
            SingleCaseES::LRRi(a$VV[a$zone == "фоновая"], a$VV[a$zone == "импактная"])) %>% 
                mutate(zone = c("буферная", "импактная"), type = V)
            }) %>% 
        map_df(rbind, .id = "year")  
        # mutate(year = substr(tur, 1, 4), tur = substr(tur, 5, 5), .before = 1)
}) %>% 
    map_df(rbind) %>% 
    as_tibble()

ES <- ES %>% map(~mutate(.x, type2 = case_when(
      type == "abu" ~ "Динамическая плотность",
      type == "nsp" ~ "Количество видов",
      type == "shan" ~ "Индекс Шеннона",
      type == "dom" ~ "Индекс доминирования",
      type == "d50" ~ "Количество видов на 50 особей", 
      TRUE ~ type)))
      # type == "size_1_small" ~ "мелкие виды",
      # type == "size_1_medium" ~ "средние виды",
      # type == "size_1_big" ~ "крупные виды",
#       type == "class_zoophaga" ~ "зоофаги",
#       type == "class_myxophaga" ~ "миксофаги", 
#       type == "habitats_meadow" ~ "луговые", 
#       type == "habitats_forest" ~ "лесные", 
#       type == "habitats_forest.meadow" ~ "эвритопные", 
#       # type == "fenol_spring" ~ "весенние",
#       # type == "fenol_autumn" ~ "осенние", 
#       type == "stratum_strato" ~ "LS",
#       type == "stratum_geohorto" ~ "HSD",
#       type == "stratum_epigeo" ~ "SS",
#       type == "stratum_stratohorto" ~ "HSD",
#       type == "fenol_all.seasons" ~ "эврихронные", 
#       type == "humid_xero" ~ "ксерофилы",
#       type == "humid_meso" ~ "мезофилы", 
#       type == "humid_hygro" ~ "гигрофилы", 
#       type == "wings_macropt" ~ "высокомобильные", 
#       type == "wings_brachypter" ~ "немобильные",
#       type == "wings_dimorph" ~ "слабомобильные")
# ))

ES %>% 
      map(function(a){a %>% 
                  rename(variable = type) %>% 
                  mutate_if(., is.numeric, round, digits = 1.5) %>% 
                  mutate_at(.vars = 1:2, as.numeric) %>% 
                  select(-ES)}) %>% 
      writexl::write_xlsx(paste0("Effect size ", Sys.Date(), ".xlsx"))

Туры разделены

Тур 1

ES$by_tur %>% 
      filter(tur == "1", type %in% c("abu", "nsp", "shan", "dom","d50")) %>% 
    mutate(
        xx = as.vector(sapply(4:0*3, function(a)(rep(1:2, each = 2) + a))),
        zone = paste0(zone, " зона")) %>%  
    ggplot(aes(x = xx, y = Est, ymin = CI_lower, ymax = CI_upper, 
               fill = year, shape= year)) +
    geom_vline(xintercept = 1:4*3, color = "grey", alpha = 0.5) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_pointrange() + 
    scale_shape_manual(values = c(25, 22)) +
    scale_x_continuous(breaks = (1:5*3)-1.5,
                         labels = unique(ES$by_tur$type2)[5:1]) +
    facet_wrap(~zone, scales = "free_x") + 
    coord_flip() + 
    theme(panel.grid = element_blank()) + 
    labs(x = NULL, y = NULL, 
         fill = "Год", shape = "Год", 
         subtitle = "Тур 1")

Тур 2

ES$by_tur %>% 
      filter(tur == "2", type %in% c("abu", "nsp", "shan", "dom","d50")) %>% 
    mutate(
        xx = as.vector(sapply(4:0*3, function(a)(rep(1:2, each = 2) + a))),
        zone = paste0(zone, " зона")) %>%  
    ggplot(aes(x = xx, y = Est, ymin = CI_lower, ymax = CI_upper, 
               fill = year, shape= year)) +
    geom_vline(xintercept = 1:4*3, color = "grey", alpha = 0.5) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_pointrange() + 
    scale_shape_manual(values = c(25, 22)) +
    scale_x_continuous(breaks = (1:5*3)-1.5,
                         labels = unique(ES$by_tur$type2)[5:1]) +
    facet_wrap(~zone, scales = "free_x") + 
    coord_flip() + 
    theme(panel.grid = element_blank()) + 
    labs(x = NULL, y = NULL, 
         fill = "Год", shape = "Год", 
         subtitle = "Тур 2")

Туры усреднены

ES$year_only %>% 
      filter(type %in% colnames(select(div, abu:d50))) %>% 
    mutate(xx = as.vector(sapply(4:0*3, function(a)(rep(1:2, each = 2) + a))), 
           zone = paste0(zone, " зона")) %>%  
    ggplot(aes(x = xx, y = Est, ymin = CI_lower, ymax = CI_upper, 
        fill = year, shape = year)) +
    geom_vline(xintercept = 1:4*3, color = "grey", alpha = 0.5) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_pointrange() + 
    scale_shape_manual(values = c(25, 22)) +
    scale_x_continuous(breaks = (1:5*3)-1.5, 
                       labels = unique(ES$by_tur$type2)[5:1]) +
        # labels = colnames(select(div, abu:d50))[5:1] ) +
    facet_wrap(~zone, scales = "free_x") + 
    coord_flip() + 
    theme(panel.grid = element_blank()) + 
    labs(x = NULL, y = NULL, 
         fill = "Год", shape = "Год", 
         subtitle = "Туры усреднены")

## Модели обилия экол.групп: ZONE 

res_zone <- list(fit = 
      colnames(div)[8:ncol(div)] %>%
      paste0(., "~tur+zone") %>%
      data.frame(x = .) %>%
      split(.$x) %>% 
      lapply(function(a){pull(a, x)}) %>% 
      lapply(FUN = function(x){lm(formula = x, data = div)})
)

res_zone$Summary <- lapply(res_zone$fit, summary) # Summary
res_zone$Plots <- res_zone %>% # Plot
      pluck("fit") %>% 
      lapply(FUN = pluck("residuals")) %>% 
      lapply(FUN = function(a){dens(a)})
res_zone$Coeff <- res_zone$Summary %>% # Coefficients
      lapply(function(a){a %>% 
                  pluck("coefficients") %>% 
                  as.data.frame() %>% 
                  rownames_to_column("var") %>% 
                  mutate(rr = 3- ceiling(log10(abs(Estimate))), 
                         est = round(Estimate, rr), 
                         pval = round(`Pr(>|t|)`, 4), 
                         pval = case_when(pval > 0 ~ as.character(pval), TRUE ~ "<0.0001"),
                         pp = paste0(est, " (p.val ", pval, ")")) %>% 
                  select(var, pp) %>% 
                  pivot_wider(names_from = var, values_from = pp) 
      })
res_zone$Shapiro <- res_zone$Summary %>% # Shapiro test
      lapply(function(a){
            s <- pluck(a, "residuals")
            data.frame(shapiro = shapiro.test(s)$p.value)
            }
      )
res_zone$Pval <- res_zone$Summary %>% # general p.val
      lapply(function(a){
            s <- pluck(a, "fstatistic")
            data.frame(pval = pf(s[1],s[2],s[3],lower.tail=F))
            }
      )
res_zone$Rsquared <- res_zone$Summary %>% # R.squared
      lapply(function(a)data.frame(rsquared = pluck(a, "r.squared")))

res_zone[4:7] %>% 
      transpose() %>% 
      lapply(function(a)map_dfc(a, cbind)) %>% 
      map_dfr(rbind, .id = "model") %>% 
      as_tibble() %>% 
      rename("zone.fon_tur1" = 2) %>% 
      mutate(shapiro = round(shapiro, 3), 
            # p.adj= round(p.adjust(pval, "BH"), 3),
            pval = round(pval, 4), 
            rsquared = round(rsquared, 2), 
            modarrange = str_extract(
                  model, 
                  paste(colnames(div)[8:ncol(div)], collapse = "|"))) %>% 
      arrange(order(colnames(div)[8:ncol(div)], modarrange)) %>% 
      select(-modarrange) %>% 
      mutate_if(is.character, function(a){
            str_replace(a, "p.val ", "") %>% 
            str_replace(" ", "<br>") %>% 
            str_replace("~", "<br>~")
            }) %>%
      formattable::formattable()



## Модели обилия экол.групп: KM
res_km <- list(fit = 
                       colnames(div)[8:ncol(div)] %>%
                       paste0(., "~tur+km") %>%
                       data.frame(x = .) %>%
                       split(.$x) %>% 
                       lapply(function(a){pull(a, x)}) %>% 
                       lapply(FUN = function(x){lm(formula = x, data = div)})
)

res_km$Summary <- lapply(res_km$fit, summary) # Summary
res_km$Plots <- res_km %>% # Plot
      pluck("fit") %>% 
      lapply(FUN = pluck("residuals")) %>% 
      lapply(FUN = function(a){dens(a)})
res_km$Coeff <- res_km$Summary %>% # Coefficients
      lapply(function(a){a %>% 
                  pluck("coefficients") %>% 
                  as.data.frame() %>% 
                  rownames_to_column("var") %>% 
                  mutate(rr = 3- ceiling(log10(abs(Estimate))), 
                         est = round(Estimate, rr), 
                         pval = round(`Pr(>|t|)`, 4), 
                         pval = case_when(pval > 0 ~ as.character(pval), TRUE ~ "<0.0001"),
                         pp = paste0(est, " (p.val ", pval, ")")) %>% 
                  select(var, pp) %>% 
                  pivot_wider(names_from = var, values_from = pp) 
      })
res_km$Shapiro <- res_km$Summary %>% # Shapiro test
      lapply(function(a){
            s <- pluck(a, "residuals")
            data.frame(shapiro = shapiro.test(s)$p.value)
      }
      )
res_km$Pval <- res_km$Summary %>% # general p.val
      lapply(function(a){
            s <- pluck(a, "fstatistic")
            data.frame(pval = pf(s[1],s[2],s[3],lower.tail=F))
      }
      )
res_km$Rsquared <- res_km$Summary %>% # R.squared
      lapply(function(a)data.frame(rsquared = pluck(a, "r.squared")))

res_km[4:7] %>% 
      transpose() %>% 
      lapply(function(a)map_dfc(a, cbind)) %>% 
      map_dfr(rbind, .id = "model") %>% 
      as_tibble() %>% 
      rename("zone.fon_tur1" = 2) %>% 
      mutate(shapiro = round(shapiro, 3), 
             # p.adj= round(p.adjust(pval, "BH"), 3),
             pval = round(pval, 4), 
             rsquared = round(rsquared, 2), 
             modarrange = str_extract(
                   model, 
                   paste(colnames(div)[8:ncol(div)], collapse = "|"))) %>% 
      arrange(order(colnames(div)[8:ncol(div)], modarrange)) %>% 
      select(-modarrange) %>% 
      mutate_if(is.character, function(a){
            str_replace(a, "p.val ", "") %>% 
                  str_replace(" ", "<br>") %>% 
                  str_replace("~", "<br>~")
      }) %>%
      formattable::formattable()

Экол. группы: effect size

Туры разделены

Первый тур

ES$by_tur %>% 
      filter(tur == "1",
            type %in% colnames(select(div, `мелкие виды`:`слабомобильные`))) %>% 
    mutate(xx = as.vector(sapply(16:0*3, function(a)(rep(1:2, each = 2) + a))), 
           zone = paste0(zone, " зона")) %>%  
    ggplot(aes(x = xx, y = Est, ymin = CI_lower, ymax = CI_upper, 
        fill = year, shape = year)) +
    geom_vline(xintercept = 1:17*3, color = "grey", linetype = "dotted") +
      geom_vline(xintercept = c(9, 18, 27, 36, 42), color = "grey")+
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_pointrange() + 
    scale_shape_manual(values = c(25, 22)) +
    scale_x_continuous(breaks = (1:17*3)-1.5,
            labels = unique(ES$by_tur$type2)[22:6]) +
    facet_wrap(~zone, scales = "free_x") + 
    coord_flip() + 
    theme(panel.grid = element_blank()) + 
    labs(x = NULL, y = NULL, 
         fill = "Год", shape = "Год", 
         subtitle = "Первый тур")

Второй тур

ES$by_tur %>% 
      filter(tur == "2",
            type %in% colnames(select(div, `мелкие виды`:`слабомобильные`))) %>% 
    mutate(xx = as.vector(sapply(16:0*3, function(a)(rep(1:2, each = 2) + a))), 
           zone = paste0(zone, " зона")) %>%  
    ggplot(aes(x = xx, y = Est, ymin = CI_lower, ymax = CI_upper, 
        fill = year, shape = year)) +
    geom_vline(xintercept = 1:16*3, color = "grey", linetype = "dotted") +
      geom_vline(xintercept = c(9, 18, 27, 36, 42), color = "grey")+
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_pointrange() + 
    scale_shape_manual(values = c(25, 22)) +
    scale_x_continuous(breaks = (1:17*3)-1.5,
        labels = unique(ES$by_tur$type2)[22:6]) +
    facet_wrap(~zone, scales = "free_x") + 
    coord_flip() + 
    theme(panel.grid = element_blank()) + 
    labs(x = NULL, y = NULL, 
         fill = "Год", shape = "Год", 
         subtitle = "Второй тур")

Туры усреднены

ES$year_only %>% 
      filter(type %in% colnames(select(div, `мелкие виды`:`слабомобильные`))) %>% 
    mutate(xx = as.vector(sapply(16:0*3, function(a)(rep(1:2, each = 2) + a))), 
           zone = paste0(zone, " зона")) %>%  
    ggplot(aes(x = xx, y = Est, ymin = CI_lower, ymax = CI_upper, 
        fill = year, shape = year)) +
    geom_vline(xintercept = 1:16*3, color = "grey", linetype = "dotted") +
      geom_vline(xintercept = c(9, 18, 27, 36, 42), color = "grey")+
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_pointrange() + 
    scale_shape_manual(values = c(25, 22)) +
    scale_x_continuous(breaks = (1:17*3)-1.5,
                       labels = unique(ES$by_tur$type2)[22:6]) +
    facet_wrap(~zone, scales = "free_x") + 
    coord_flip() + 
    theme(panel.grid = element_blank()) + 
    labs(x = NULL, y = NULL, 
         fill = "Год", shape = "Год", 
         subtitle = "Туры усреднены")

Ординация

Динамическая плотность

Столбцы - туры, строки - годы учёта.

Туры разделены

# save.image(paste0("for_tables_", Sys.Date(), ".RData")) !!!!

wide <- wide %>% 
      filter(apply(wide[,6:ncol(wide)], 1, sum) > 0)
pcoa2 <- wide %>% 
    column_to_rownames("id") %>% 
    select(-1:-4) %>% 
    vegan::vegdist() %>% 
    ape::pcoa()
eig2 <- pcoa2$values$Eigenvalues
if(min(eig2) < 0) {eig2 <- eig2 + abs(min(eig2))}
eig2 <- round(eig2/sum(eig2)*100)

По-новому

pcoa2$vectors[,1:2] %>% 
    as.data.frame() %>% 
    rownames_to_column("id") %>% 
    mutate(year = as.numeric(substr(id, 1, 4)), 
           tur = factor(substr(id, 5, 5)),
           site = substr(id, 6, 9), 
           plot = as.numeric(substr(id, 10, 10))) %>% 
    as_tibble() %>% 
    left_join(select(div, year, tur, site, plot, zone), 
              by = c("year", "tur", "site", "plot")) %>% 
      mutate(Axis.1 = Axis.1 * -1, 
             tur = paste0("Тур ", tur), 
             zone = fct_relabel(zone, ~paste0(.x, " зона")),
             year = factor(year) 
             ) %>% 
      ggplot(aes(x = Axis.1, y = Axis.2, fill = year, shape = year, size = tur)) +
      geom_point(color = "black") + 
      stat_ellipse(mapping = aes(color = year), linewidth = 0.5, linetype = "dashed") +
      facet_grid(cols = vars(zone), rows = vars(tur)) +
      scale_size_manual(values = c(2,2)) +
      scale_shape_manual(values = c(21, 22, 23, 25)) + 
      # scale_color_manual(values = c("green", "orange", "red", "brown")) + 
      # scale_fill_manual(values = c("green", "orange", "red", "brown")) + 
      labs(subtitle = "По динамической плотности", 
           x = paste0("Ось 1 (", eig2[1], " %)"), 
           y = paste0("Ось 2 (", eig2[2], " %)"), 
           fill = "Год", color = "Год", shape = "Год") +
      ggplot2::guides(size = "none", color = "none")

Как раньше

pcoa2$vectors[,1:2] %>%     
    as.data.frame() %>% 
    rownames_to_column("id") %>% 
    mutate(year = as.numeric(substr(id, 1, 4)), 
           tur = factor(substr(id, 5, 5)),
           site = substr(id, 6, 9), 
           plot = as.numeric(substr(id, 10, 10))) %>% 
    as_tibble() %>% 
    left_join(select(div, year, tur, site, plot, zone), 
              by = c("year", "tur", "site", "plot")) %>% 
      mutate(Axis.1 = Axis.1 * -1, 
             tur = fct_relabel(tur, ~paste0("Тур ", .x))) %>% 
      ggplot(aes(x = Axis.1, y = Axis.2, fill = zone, shape = zone, size = tur)) +
      geom_point(color = "black") + 
      stat_ellipse(mapping = aes(color = zone), size = 0.5, linetype = "dashed") +
      facet_grid(cols = vars(year), rows = vars(tur)) +
      scale_size_manual(values = c(2,2)) +
      scale_shape_manual(values = c(21, 22, 23, 25)) + 
      scale_color_manual(values = c("green", "orange", "red", "brown")) + 
      scale_fill_manual(values = c("green", "orange", "red", "brown")) + 
      labs(subtitle = "По динамической плотности", 
           x = paste0("Ось 1 (", eig2[1], " %)"), 
           y = paste0("Ось 2 (", eig2[2], " %)"), 
           fill = "Зона", color = "Зона", shape = "Зона") +
      ggplot2::guides(size = "none", color = "none")

Туры объединены

pcoa3 <- wide %>% 
      select(-id, -tur) %>% 
      group_by(year, zone, site, plot) %>% 
      summarise_all(sum) %>% 
      unite(id, 1:4,  sep = "_") %>% 
    column_to_rownames("id") %>% 
    vegan::vegdist() %>% 
    ape::pcoa()
eig3 <- pcoa3$values$Eigenvalues
if(min(eig3) < 0) {eig3 <- eig3 + abs(min(eig3))}
eig3 <- round(eig3/sum(eig3)*100)

pcoa3 <- pcoa3$vectors[,1:2] %>% 
    as.data.frame() %>% 
    rownames_to_column("id") %>% 
      as_tibble() %>% 
      separate(id, into = c("year", "zone", "site", "plot"), sep = "_") %>% 
      mutate(zone = factor(zone, levels = c("фоновая", "буферная", "импактная")))

По-новому

pcoa3 %>% 
      mutate(zone = fct_relabel(zone, ~paste0(.x, " зона"))) %>% 
      ggplot(aes(x = Axis.1, y = Axis.2, fill = year, color = year, shape = year)) +
      geom_point(color = "black") +  
      stat_ellipse(mapping = aes(color = year), linewidth = 0.5, linetype = "dashed") +
      facet_grid(cols = vars(zone)) +
      scale_shape_manual(values = c(21, 22))+ 
      labs(subtitle = "По динамической плотности", 
           x = paste0("Ось 1 (", eig3[1], " %)"), 
           y = paste0("Ось 2 (", eig3[2], " %)"), 
           fill = "Год", color = "Год", shape = "Год") +
      ggplot2::guides(size = "none", color = "none")

Как раньше

pcoa3 %>% 
      ggplot(aes(x = Axis.1, y = Axis.2, fill = zone, shape = zone)) +
      geom_point(color = "black") + 
      stat_ellipse(mapping = aes(color = zone), size = 0.5, linetype = "dashed") +
      facet_grid(cols = vars(year)) +
      scale_size_manual(values = c(2,2)) +
      scale_shape_manual(values = c(21, 22, 23, 25)) + 
      scale_color_manual(values = c("green", "orange", "red", "brown")) + 
      scale_fill_manual(values = c("green", "orange", "red", "brown")) + 
      labs(subtitle = "По динамической плотности", 
           x = paste0("Ось 1 (", eig2[1], " %)"), 
           y = paste0("Ось 2 (", eig2[2], " %)"), 
           fill = "Зона", color = "Зона", shape = "Зона") +
      ggplot2::guides(size = "none", color = "none")

Лепестковые диаграммы

dummy <- c(
    "myxophaga", "Mx",
    "zoophaga", "Z",
    # "all.seasons", 
    # "autumn", 
    # "spring", 
    "forest", "F",
    "forest.meadow", "E",
    "meadow", "Oh",
    "hygro", "hyg",
    "meso", "mes",
    "xero", "xer", ##
    "big", "B",
    "medium", "M",
    "small", "S",
    "epigeo", "SS",
    "geohorto", "HSD",
    "stratohorto", "HSD",
    "strato", "LS",
    "brachypter", "lm",
    "dimorph", "mm",
    "macropt", "hm"
) %>% 
    matrix(ncol = 2, byrow = TRUE, dimnames = list(NULL, c("grp", "grp2"))) %>% 
    as.data.frame() 

spider_year <- long %>% 
    left_join(tr, by = "taxa") %>% 
    select(-site, -plot, -taxa, -stratum.ext, -abu, -habitats.ext, -fenol) %>% 
    pivot_longer(names_to = "cls", values_to = "grp", -c(1:4)) %>% 
    # filter(!is.na(grp), grp != "stratohorto") %>% 
    right_join(dummy, by = "grp") %>% 
    group_by(year, zone, cls, grp2) %>% 
    summarise(abtotal = sum(num), .groups = "drop_last") %>% 
    mutate(abtotal = round(abtotal/sum(abtotal)*100, 1)) %>% 
    ungroup() %>% 
    select(-cls) %>% 
    pivot_wider(names_from = grp2, values_from = abtotal, values_fill = 0)
    
spider_turyear <- long %>% 
    left_join(tr, by = "taxa") %>% 
    select(-site, -plot, -taxa, -stratum.ext, -abu, -habitats.ext) %>% 
    pivot_longer(names_to = "cls", values_to = "grp", -c(1:4)) %>% 
    # filter(!is.na(grp), grp != "stratohorto") %>% 
    right_join(dummy, by = "grp") %>% 
    group_by(year, tur, zone, cls, grp2) %>% 
    summarise(abtotal = sum(num), .groups = "drop_last") %>% 
    mutate(abtotal = round(abtotal/sum(abtotal)*100, 1)) %>% 
    ungroup() %>% 
    select(-cls) %>% 
    pivot_wider(names_from = grp2, values_from = abtotal, values_fill = 0)

dummy %>% 
      select(`Сокращение__` = grp2, `Обозначение в БД__` = grp) %>% 
      mutate(`Описание` = c("Питание: Миксофаги", "Питание: Зоофаги", 
                            "Биотопический преферендум: Лесные", 
                            "Биотопический преферендум: Эвритопные", 
                            "Биотопический преферендум: Луговые", 
                            "Отношение к влажности: Гигрофилы", 
                            "Отношение к влажности: Мезофилы", 
                            "Отношение к влажности: Ксерофилы", 
                            "Размерный класс: Крупные (>11 мм)", 
                            "Размерный класс: Средние (9-11 мм)", 
                            "Размерный класс: Мелкие (<9 мм)",
                            "Предпочитаемый ярус: поверхность почвы (эпигеобионты)", 
                            "Предпочитаемый ярус: растительность (хорто-, тамно- и дендробионты)",
                            "Предпочитаемый ярус: растительность (хорто-, тамно- и дендробионты)",
                            "Предпочитаемый ярус: подстилка, почва (стратобионт)", 
                            "Способность к расселению: низкая", 
                            "Способность к расселению: умеренная", 
                            "Способность к расселению: высковая")) %>% 
      
knitr::kable("html", caption = "Расшифровка сокращений", align = "ccl")

Расшифровка сокращений

Сокращение__	Обозначение в БД__	Описание
Mx	myxophaga	Питание: Миксофаги
Z	zoophaga	Питание: Зоофаги
F	forest	Биотопический преферендум: Лесные
E	forest.meadow	Биотопический преферендум: Эвритопные
Oh	meadow	Биотопический преферендум: Луговые
hyg	hygro	Отношение к влажности: Гигрофилы
mes	meso	Отношение к влажности: Мезофилы
xer	xero	Отношение к влажности: Ксерофилы
B	big	Размерный класс: Крупные (>11 мм)
M	medium	Размерный класс: Средние (9-11 мм)
S	small	Размерный класс: Мелкие (<9 мм)
SS	epigeo	Предпочитаемый ярус: поверхность почвы (эпигеобионты)
HSD	geohorto	Предпочитаемый ярус: растительность (хорто-, тамно- и дендробионты)
HSD	stratohorto	Предпочитаемый ярус: растительность (хорто-, тамно- и дендробионты)
LS	strato	Предпочитаемый ярус: подстилка, почва (стратобионт)
lm	brachypter	Способность к расселению: низкая
mm	dimorph	Способность к расселению: умеренная
hm	macropt	Способность к расселению: высковая

Туры объединены

Фоновая зона

op <- par(family = "serif", font = 1) 
spider_year %>% 
    filter(zone == "фоновая") %>% 
    select(-zone) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Фоновая зона", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     # family = "serif",
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Буферная зона

op <- par(family = "serif", font = 1) 
spider_year %>% 
    filter(zone == "буферная") %>% 
    select(-zone) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Буферная зона", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Импактная зона

op <- par(family = "serif", font = 1) 
spider_year %>% 
    filter(zone == "импактная") %>% 
    select(-zone) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Импактная зона", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Первый тур

Фоновая зона

op <- par(family = "serif", font = 1) 
spider_turyear %>% 
    filter(zone == "фоновая", tur == "1") %>% 
    select(-zone, -tur) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Фоновая зона. Тур 1", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Буферная зона

op <- par(family = "serif", font = 1) 
spider_turyear %>% 
    filter(zone == "буферная", tur == "1") %>% 
    select(-zone, -tur) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Буферная зона. Тур 1", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Импактная зона

op <- par(family = "serif", font = 1) 
spider_turyear %>% 
    filter(zone == "импактная", tur == "1") %>% 
    select(-zone, -tur) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Импактная зона. Тур 1", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Второй тур

Фоновая зона

op <- par(family = "serif", font = 1) 
spider_turyear %>% 
    filter(zone == "фоновая", tur == "2") %>% 
    select(-zone, -tur) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Фоновая зона. Тур 2", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Буферная зона

op <- par(family = "serif", font = 1) 
spider_turyear %>% 
    filter(zone == "буферная", tur == "2") %>% 
    select(-zone, -tur) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Буферная зона. Тур 2", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )

Импактная зона

op <- par(family = "serif", font = 1) 
spider_turyear %>% 
    filter(zone == "импактная", tur == "2") %>% 
    select(-zone, -tur) %>% 
    column_to_rownames("year") %>% 
    rbind(rep(round(max(.), -2), ncol(.)), 
          rep(0, ncol(.)), 
          .) %>% 
    select( LS, SS, mes, hyg, xer, `F`, Oh, E, Mx, Z, 
       hm, mm, lm, S, M, B, HSD
) %>% 
    fmsb::radarchart(title = "Импактная зона. Тур 2", 
                     axistype = 1, 
                     caxislabels=seq(0, 100, by = 25), #
                     axislabcol = "darkgrey", #make darker
                     cglcol = "darkgrey", 
                     pcol = c(rgb(1, 0.5, 0.15, 1), rgb(0, 0.8, 0.4, alpha = 1)),
                     pfcol = c(rgb(1, 0.5, 0.15, 0.2), rgb(0, 0.8, 0.4, alpha = 0.4)), 
                     plty = c(1,1)
    )#; legend(

    #     x = 0.9, y=0.9,  pch=19, bty  = "n",
    #     legend = c("2014","2009"),
    #     col = c(rgb(0, 0.8, 0.4), rgb(1, 0.5, 0.15))
    # )