Загрузка данных

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(tidyverse) # Проблемы: 2000 ничем не окружен, цикл с 2005 по 2008 четырехлетний
library(foreach)
library(parallel)
library(doParallel)
nspec <- function(d, go = 1, finish = ncol(d)){ 
    apply(d[,go:finish], 1, FUN = function(a){length(a[a>0])})
}
wide.df <- readxl::read_xlsx("Mukhacheva_data_18022021.xlsx") %>% 
    filter(field == 1) %>%  # there were no samples at 6&30 before 2004
    select(-field)
long.df <- wide.df %>% 
    pivot_longer(names_to = "spec", values_to = "val", 
               -c("id", "year", "period", "cycle", "zone", "km", "total"))

Рис. 1. Ординация

Обилие*цикл

w.cycle.n <- long.df %>%
    # filter(total > 0) %>% 
    group_by(zone, cycle, spec) %>%
    summarise(val = mean(val), .groups = "drop") %>%
    pivot_wider(names_from = spec, values_from = val) %>%
    as.data.frame() %>% 
    filter(cycle != "2000")
dist.cn <- vegan::vegdist(w.cycle.n[,3:16], method = "bray", binary = FALSE)
pcoa.cn <- cbind(w.cycle.n[,1:2], ape::pcoa(dist.cn)$vectors[,1:2]) 
eig.cn <- data.frame(eig = ape::pcoa(dist.cn)$values$Eigenvalues) %>% filter(eig > 0) %>% 
    mutate(eig = round(eig/sum(eig)*100)) %>% slice(1:5) %>% pull(eig)
cn <- pcoa.cn %>% 
    ggplot(aes(x = Axis.1, y = Axis.2, fill = zone, shape = zone)) + 
    geom_path(mapping = aes(x = Axis.1, y = Axis.2, data = ), size = 1, alpha = 0.8) +
    geom_point(size = 3)+
    geom_text(mapping = aes(label = cycle), color = "black", hjust=0, vjust=0) +
    scale_shape_manual(values = c(21, 24, 22)) +
    scale_fill_manual(values = c("white", "grey", "black")) +
    theme_bw() +
    theme(legend.position = "bottom", panel.grid.major = element_blank(), 
         panel.grid.minor = element_blank()) +
        labs(title = "А", 
        x = paste0("Axis 1, ", eig.cn[1], "%"), y = paste0("Axis 2, ", eig.cn[2], "%"))
cn

Доли * цикл

w.cycle.p <- long.df %>%
    filter(total > 0) %>% 
    mutate(val = val/total) %>%
    group_by(zone, cycle, spec) %>%
    summarise(val = mean(val), .groups = "drop") %>%
    #mutate(val = case_when(is.na(val) ~ 1, TRUE ~ val)) %>% #####
    pivot_wider(names_from = spec, values_from = val)   %>%
    as.data.frame() %>% 
    filter(cycle != "2000")
dist.cp <- vegan::vegdist(w.cycle.p[,3:16], method = "bray", binary = FALSE)
pcoa.cp <- cbind(w.cycle.p[,1:2], ape::pcoa(dist.cp)$vectors[,1:2])
eig.cp <- data.frame(eig = ape::pcoa(dist.cp)$values$Eigenvalues) %>% filter(eig > 0) %>% 
    mutate(eig = round(eig/sum(eig)*100)) %>% slice(1:5) %>% pull(eig)
cp <- pcoa.cp %>% 
    mutate(Axis.2 = Axis.2*-1) %>% 
    ggplot(aes(x = Axis.1, y = Axis.2, fill = zone, shape = zone)) + 
    geom_path(mapping = aes(x = Axis.1, y = Axis.2, data = ), size = 1, alpha = 0.8) +
    geom_point(size = 3)+
    geom_text(mapping = aes(label = cycle), color = "black", hjust=0, vjust=0) +
    scale_shape_manual(values = c(21, 24, 22)) +
    scale_fill_manual(values = c("white", "grey", "black")) +
    theme_bw() +
    theme(legend.position = "bottom", panel.grid.major = element_blank(), 
         panel.grid.minor = element_blank()) +
        labs(title = "Б", 
        x = paste0("Axis 1, ", eig.cp[1], "%"), y = paste0("Axis 2, ", eig.cp[2], "%"))
cp

Для публикации

plot1 <- gridExtra::grid.arrange(cn, cp, ncol = 1)

arj <- list(w.cycle.n = w.cycle.n, w.cycle.p = w.cycle.p, dist.cn = dist.cn, dist.cp = dist.cp, 
          pcoa.cn = pcoa.cn, pcoa.cp = pcoa.cp, eig.cn = eig.cn, eig.cp = eig.cp)
rm(w.cycle.n, w.cycle.p, cn, cp, dist.cn, dist.cp, pcoa.cn, pcoa.cp, eig.cn, eig.cp)
ggsave(plot = plot1, filename = "rpot1.svg", width = 6.5, height = 9, dpi = 900)

Табл. 1. PERMANOVA

Доли*год

w.year.p <- long.df %>%
    filter(total>0) %>% 
    mutate(val = val/total) %>%
    group_by(zone, year, spec, period) %>%
    summarise(val = mean(val), .groups = "drop") %>%
    pivot_wider(names_from = spec, values_from = val) %>%
    mutate(period = as.factor(period)) %>% 
    filter(year != 2000, year != 1996, year != 2012)
dist.yp <- vegan::vegdist(w.year.p[,4:17], method = "bray", binary = FALSE)
set.seed(25)
tab1 <- vegan::adonis(dist.yp ~ zone * period, data = w.year.p, permutations = 9999) ################ MUST BE 9999
# arj$w.year.p <- w.year.p
# arj$dist.yp <- dist.yp
#rm(w.year.p, dist.yp)
tab1

## 
## Call:
## vegan::adonis(formula = dist.yp ~ zone * period, data = w.year.p,      permutations = 9999) 
## 
## Permutation: free
## Number of permutations: 9999
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## zone         2    2.4227 1.21134 14.0705 0.23893 0.0001 ***
## period       2    0.8695 0.43474  5.0498 0.08575 0.0005 ***
## zone:period  4    0.9072 0.22681  2.6346 0.08947 0.0062 ** 
## Residuals   69    5.9403 0.08609         0.58584           
## Total       77   10.1397                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Рис. 3. Пространственная изменчивость во времени

Параллельные вычисления

cl <- makeCluster(16) # up to 16
doParallel::registerDoParallel(cl)
parallel::clusterEvalQ(cl, set.seed(10))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
## 
## [[7]]
## NULL
## 
## [[8]]
## NULL
## 
## [[9]]
## NULL
## 
## [[10]]
## NULL
## 
## [[11]]
## NULL
## 
## [[12]]
## NULL
## 
## [[13]]
## NULL
## 
## [[14]]
## NULL
## 
## [[15]]
## NULL
## 
## [[16]]
## NULL

r3 <- foreach(i=1:10000) %dopar% { ################ MUST BE 5000
# sample
    library(dplyr)
    if(sample(1:0, 1) == 1) {
        dfa <- wide.df %>% 
            group_by(year) %>% 
            sample_n(size = n()-1) %>% 
            ungroup
    } else { 
        dfa <- wide.df %>% 
            select(1:7, sort(sample(8:21, 13)))
    }
# all gradient year by year: BW
    bw.years <- dfa %>% 
        select(-id, -period, -cycle, -zone, -total) %>% 
        group_by(year, km) %>% 
        summarise_all(sum) %>% 
        ungroup %>% 
        mutate(ss = nspec(., 3, ncol(.))) %>% 
        group_by(year) %>% 
        summarise_all(mean) %>% 
        transmute(year, by.km = nspec(., 3, ncol(.) - 1) / ss)
# all gradient year by year: Ibc
    dis <- as.matrix(vegan::vegdist(dfa[,8:ncol(dfa)], method = "bray", binary = FALSE))
    rownames(dis) <- colnames(dis) <- dfa$id
    dis <- reshape2::melt(dis, varnames = c("id1", "id2")) %>% 
        as_tibble() %>% 
        mutate(id1 = as.character(id1), id2 = as.character(id2), 
            value = case_when(is.nan(value) ~ 1, TRUE ~ value)) %>% 
        filter(as.numeric(substr(id1, 4, nchar(id1))) > as.numeric(substr(id2, 4, nchar(id2)))) %>% 
        left_join(select(dfa, id1 = id, km1 = km, year1 = year, zone1 = zone), by = "id1") %>% 
        left_join(select(dfa, id2 = id, km2 = km, year2 = year, zone2 = zone), by = "id2") 
    ibc.years <- dis %>% 
        filter(year1 == year2) %>% 
        group_by(year1) %>% 
        summarise(val = mean(value), .groups = "drop")
# Results
    list(bw.years = bw.years, 
        ibc.years = ibc.years)
}
r3 <- transpose(r3)

Computing is hidden

Расчёт показателей

fit <- list()
bw.years <- map_dfr(r3$bw.years, rbind) %>% 
    group_by(year) %>% 
    summarise(low = as.numeric(quantile(by.km, 0.025)), up = as.numeric(quantile(by.km, 0.975)), .groups = "drop")
bw.years <- wide.df %>% 
    select(-id, -period, -cycle, -zone, -total) %>% 
    group_by(year, km) %>% 
    summarise_all(mean) %>% 
    ungroup %>% 
    mutate(ss = nspec(., 3, ncol(.))) %>% 
    group_by(year) %>%
    summarise_all(mean) %>% 
    transmute(year, B = nspec(., 3, ncol(.) - 1) / ss) %>% 
    left_join(bw.years, by = "year")
fit$bw1 <- lm(B ~ year, data = bw.years)
fit$bw2 <- bw.years %>% mutate(y2 = year^2) %>% lm(B ~ year + y2 , data = .)
bw.years <- bw.years %>% mutate(y2 = year^2) %>% 
    transmute(year, B, low, up, B2 = as.numeric(predict(fit$bw2, newdata = .)), type = "bw")
ibc.years <- as.matrix(vegan::vegdist(wide.df[,8:ncol(wide.df)], method = "bray", binary = FALSE))
rownames(ibc.years) <- colnames(ibc.years) <- wide.df$id
ibc.years <- reshape2::melt(ibc.years, varnames = c("id1", "id2")) %>% 
    as_tibble() %>% 
    mutate(id1 = as.character(id1), id2 = as.character(id2)) %>% 
    filter(as.numeric(substr(id1, 4, nchar(id1))) > as.numeric(substr(id2, 4, nchar(id2)))) %>% 
    left_join(select(wide.df, id1 = id, km1 = km, year1 = year, zone1 = zone), by = "id1") %>% 
    left_join(select(wide.df, id2 = id, km2 = km, year2 = year, zone2 = zone), by = "id2") %>% 
    filter(year1 == year2) %>% 
    mutate(value = case_when(is.nan(value) ~ 1, TRUE ~ value)) %>% 
    group_by(year1) %>% 
    summarise(B = mean(value), .groups = "drop") 
ibc.years <- map_dfr(r3$ibc.years, rbind) %>% 
    group_by(year1) %>% 
    summarise(low = as.numeric(quantile(val, 0.025)), up = as.numeric(quantile(val, 0.975)), .groups = "drop") %>% 
    left_join(ibc.years, by = "year1") %>% 
    transmute(year = year1, B, low, up)
fit$ibc1 <- lm(B ~ year, data = ibc.years)
fit$ibc2 <- ibc.years %>% mutate(y2 = year^2) %>% lm(B ~ year + y2 , data = .)
ibc.years <- ibc.years %>% mutate(y2 = year^2) %>% 
    transmute(year, B, low, up, B2 = as.numeric(predict(fit$ibc2, newdata = .)), type = "ibc")
years <- rbind(bw.years, ibc.years) %>% mutate(zone = "0.all")
bw.z <- wide.df %>% 
    select(-id, -period, -cycle, -total) %>% 
    group_by(year, zone, km) %>% 
    summarise_all(mean) %>% 
    ungroup %>% 
    mutate(ss = nspec(., 4, ncol(.))) %>% 
    group_by(year, zone) %>%
    summarise_all(mean) %>% 
    ungroup() %>% 
    transmute(year, B = nspec(., 3, ncol(.) - 1) / ss, low = NA, up = NA, B2 = NA, type = "bw", zone) %>% 
    filter(zone == "3.imp" | ((year >= 2004) & (zone == "1.fon")) | ((year >= 2009) & (zone == "2.buf")) , is.finite(B)) 
ibc.z <- as.matrix(vegan::vegdist(wide.df[,8:ncol(wide.df)], method = "bray", binary = FALSE))
rownames(ibc.z) <- colnames(ibc.z) <- wide.df$id
ibc.z <- reshape2::melt(ibc.z, varnames = c("id1", "id2")) %>% 
    as_tibble() %>% 
    mutate(id1 = as.character(id1), id2 = as.character(id2), 
          value = case_when(is.nan(value) ~ 1, TRUE ~ value)) %>% 
    filter(as.numeric(substr(id1, 4, nchar(id1))) > as.numeric(substr(id2, 4, nchar(id2)))) %>% 
    left_join(select(wide.df, id1 = id, km1 = km, year1 = year, zone1 = zone), by = "id1") %>% 
    left_join(select(wide.df, id2 = id, km2 = km, year2 = year, zone2 = zone), by = "id2") %>% 
    filter(zone1 == zone2, year1 == year2) %>% 
    group_by(zone1, year1) %>% 
    summarise(ibc = mean(value), .groups = "drop") %>% 
    filter(zone1 == "3.imp" | year1 >= 2004) %>% 
    transmute(year = year1, B = ibc, low = NA, up = NA, B2 = NA, type = "ibc", zone = zone1)
years <- rbind(years, bw.z, ibc.z)

Computing is hidden

Бета Уиттекера

plot3a

p3a <- years %>% filter(type == "bw") %>% 
    ggplot(aes(x = year, y = B, ymin = low, ymax = up)) + 
    geom_ribbon(alpha = 0.4, fill = "lightgrey", color = "darkgrey") + 
    geom_line(color = "black") + 
    geom_point(size = 2.5, color = "black", fill = "darkgrey", shape = 21) + 
    facet_grid(rows = vars(zone)) +
    theme_bw() + 
    theme(legend.position = "none", 
         axis.text.x = element_text(angle=90), 
         panel.grid.minor.x = element_blank()) + 
    scale_x_continuous(breaks = 1990:2020) + 
    coord_fixed(ratio = 1.5) +
    geom_vline(aes(xintercept = 1997.5), size = 1, linetype = "dashed") +
    geom_vline(aes(xintercept = 2009.5), size = 1, linetype = "dashed") +
    labs(x = "", y = "Бета Уиттекера", title = "А")
p3a

plot3a + trends

p3a + 
    geom_line(mapping = aes(x = year, y = B2), linetype = "dashed", data = filter(years, type == "bw") ) + 
    geom_abline(mapping = aes(slope = ss, intercept = ii), linetype = "dashed", data = data.frame(
        ss = fit$bw1$coefficients[[2]], ii = fit$bw1$coefficients[[1]], zone = "0.all"))

Bw: линейная модель

Информационный критерий Акаике = 44.25

Параметры модели:

summary(fit$bw1)

## 
## Call:
## lm(formula = B ~ year, data = bw.years)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7849 -0.2874 -0.1061  0.1435  1.4813 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.534755  19.750138  -0.280    0.781
## year         0.003754   0.009849   0.381    0.706
## 
## Residual standard error: 0.4849 on 27 degrees of freedom
## Multiple R-squared:  0.005353,   Adjusted R-squared:  -0.03149 
## F-statistic: 0.1453 on 1 and 27 DF,  p-value: 0.706

Bw: полиномиальная модель

Информационный критерий Акаике = 46.07

Параметры модели:

summary(fit$bw2)

## 
## Call:
## lm(formula = B ~ year + y2, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7444 -0.2952 -0.1100  0.1623  1.4981 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.019e+03  5.085e+03   0.397    0.695
## year        -2.016e+00  5.073e+00  -0.397    0.694
## y2           5.037e-04  1.265e-03   0.398    0.694
## 
## Residual standard error: 0.4927 on 26 degrees of freedom
## Multiple R-squared:  0.01138,    Adjusted R-squared:  -0.06467 
## F-statistic: 0.1496 on 2 and 26 DF,  p-value: 0.8618

Расстояние Брея-Кёртиса

plot3b

p3b <- years %>% filter(type == "ibc") %>% 
    ggplot(aes(x = year, y = B, ymin = low, ymax = up)) + 
    geom_ribbon(alpha = 0.4, fill = "lightgrey", color = "darkgrey") + 
    geom_line(color = "black") + 
    geom_point(size = 2.5, color = "black", fill = "darkgrey", shape = 21) + 
    facet_grid(rows = vars(zone)) +
    theme_bw() + 
    theme(legend.position = "none", 
         axis.text.x = element_text(angle=90), 
         panel.grid.minor.x = element_blank()) + 
    scale_x_continuous(breaks = 1990:2020) + 
    coord_fixed(ratio = 8.1) +
    geom_vline(aes(xintercept = 1997.5), size = 1, linetype = "dashed") +
    geom_vline(aes(xintercept = 2009.5), size = 1, linetype = "dashed") +
    labs(x = "", y = "Расстояние Брея-Кёртиса", title = "Б")
p3b

plot3b + trends

p3b + 
    geom_line(mapping = aes(x = year, y = B2), linetype = "dashed", data = filter(years, type == "ibc") ) + 
    geom_abline(mapping = aes(slope = ss, intercept = ii), linetype = "dashed", data = data.frame(
        ss = fit$ibc1$coefficients[[2]], ii = fit$ibc1$coefficients[[1]], zone = "0.all"))

Ibc: линейная модель

Информационный критерий Акаике = -29.73

Параметры модели:

summary(fit$ibc1)

## 
## Call:
## lm(formula = B ~ year, data = ibc.years)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22870 -0.08384 -0.02982  0.09614  0.33846 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.963305   5.515452  -1.081    0.289
## year         0.003278   0.002750   1.192    0.244
## 
## Residual standard error: 0.1354 on 27 degrees of freedom
## Multiple R-squared:  0.04999,    Adjusted R-squared:  0.01481 
## F-statistic: 1.421 on 1 and 27 DF,  p-value: 0.2436

Ibc: полиномиальная модель

Информационный критерий Акаике = -27.85

Параметры модели:

summary(fit$ibc2)

## 
## Call:
## lm(formula = B ~ year + y2, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22481 -0.08816 -0.03015  0.09412  0.33791 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.579e+02  1.422e+03  -0.322    0.750
## year         4.541e-01  1.418e+00   0.320    0.751
## y2          -1.124e-04  3.537e-04  -0.318    0.753
## 
## Residual standard error: 0.1377 on 26 degrees of freedom
## Multiple R-squared:  0.05367,    Adjusted R-squared:  -0.01912 
## F-statistic: 0.7373 on 2 and 26 DF,  p-value: 0.4881

Всё вместе

gridExtra::grid.arrange(p3a, p3b, ncol = 2, top = "Без линий тренда")

gridExtra::grid.arrange(p3a +
    geom_line(mapping = aes(x = year, y = B2), linetype = "dashed", data = filter(years, type == "bw") ) + 
    geom_abline(mapping = aes(slope = ss, intercept = ii), linetype = "dashed", data = data.frame(
        ss = fit$bw1$coefficients[[2]], ii = fit$bw1$coefficients[[1]], zone = "0.all")),
                    p3b +
    geom_line(mapping = aes(x = year, y = B2), linetype = "dashed", data = filter(years, type == "ibc") ) + 
    geom_abline(mapping = aes(slope = ss, intercept = ii), linetype = "dashed", data = data.frame(
        ss = fit$ibc1$coefficients[[2]], ii = fit$ibc1$coefficients[[1]], zone = "0.all")),
    ncol = 2, top = "С линиями тренда")

Таблицы с данными

Бета Уиттекера, всесь градиент

years %>% filter(type == "bw", zone == "0.all") %>% 
    transmute(year, betadiv = round(B, 2), low = round(low, 2), up = round(up, 2)) %>% 
    formattable::formattable()

year	betadiv	low	up
1990	1.90	1.54	2.08
1991	2.08	1.50	2.25
1992	1.88	1.45	2.27
1993	2.06	1.82	2.50
1994	1.71	1.20	2.25
1995	1.67	1.42	1.79
1996	2.67	2.00	4.00
1997	2.11	1.67	2.40
1998	2.22	1.80	2.57
2000	2.50	1.60	3.00
2002	1.88	1.43	2.40
2003	1.20	1.00	1.50
2004	1.44	1.20	1.52
2005	1.75	1.50	2.00
2006	1.64	1.29	1.93
2007	1.48	1.35	1.57
2008	1.43	1.33	1.50
2009	2.88	2.47	3.50
2010	1.80	1.60	1.93
2011	2.15	1.83	2.58
2012	3.50	3.00	4.20
2013	2.21	1.89	2.50
2014	1.65	1.50	1.81
2015	1.87	1.57	2.13
2016	2.14	1.71	2.41
2017	1.75	1.09	2.00
2018	2.50	2.00	3.50
2019	1.93	1.55	2.12
2020	1.84	1.64	2.03

Бета Уиттекера, по зонам

years %>% filter(type == "bw", zone != "0.all") %>% 
    transmute(year, B = round(B, 2), zone) %>% 
    reshape2::dcast(., year ~ zone, value.var="B") %>% 
    formattable::formattable()

year	1.fon	2.buf	3.imp
1990	NA	NA	1.69
1991	NA	NA	2.10
1992	NA	NA	1.71
1993	NA	NA	2.00
1994	NA	NA	1.33
1995	NA	NA	1.88
1997	NA	NA	2.00
1998	NA	NA	1.75
2000	NA	NA	2.00
2002	NA	NA	2.00
2003	NA	NA	2.00
2004	1.54	NA	1.14
2005	1.45	NA	1.50
2006	1.18	NA	2.00
2007	1.33	NA	1.20
2008	1.71	NA	1.33
2009	1.78	2.00	3.00
2010	1.50	1.50	1.91
2011	1.75	1.40	1.88
2012	2.00	2.00	NA
2013	1.56	2.00	3.00
2014	1.33	1.33	1.85
2015	1.50	1.33	1.67
2016	1.43	1.75	1.50
2017	1.78	1.25	1.36
2018	2.00	1.67	2.40
2019	1.38	1.23	1.64
2020	1.25	1.33	1.80

Расстояние Брея-Кёртиса, всесь градиент

years %>% filter(type == "ibc", zone == "0.all") %>% 
    transmute(year, betadiv = round(B, 2), low = round(low, 2), up = round(up, 2)) %>% 
    formattable::formattable()

year	betadiv	low	up
1990	0.44	0.38	0.51
1991	0.71	0.59	0.81
1992	0.59	0.38	0.69
1993	0.54	0.34	0.85
1994	0.35	0.28	0.79
1995	0.42	0.29	0.66
1996	0.92	0.84	1.00
1997	0.69	0.48	0.81
1998	0.56	0.41	0.83
2000	0.73	0.20	1.00
2002	0.53	0.41	0.94
2003	0.52	0.33	0.71
2004	0.41	0.26	0.47
2005	0.70	0.60	0.74
2006	0.68	0.39	0.82
2007	0.53	0.45	0.58
2008	0.55	0.38	0.65
2009	0.74	0.63	0.85
2010	0.61	0.55	0.66
2011	0.55	0.41	0.66
2012	0.92	0.89	0.96
2013	0.73	0.63	0.80
2014	0.55	0.47	0.66
2015	0.56	0.44	0.66
2016	0.63	0.53	0.70
2017	0.57	0.51	0.63
2018	0.78	0.70	0.91
2019	0.61	0.52	0.66
2020	0.60	0.52	0.68

Расстояние Брея-Кёртиса, по зонам

years %>% filter(type == "ibc", zone != "0.all") %>% 
    transmute(year, B = round(B, 2), zone) %>% 
    reshape2::dcast(., year ~ zone, value.var="B") %>% 
    formattable::formattable()

year	1.fon	2.buf	3.imp
1990	NA	NA	0.33
1991	NA	NA	0.86
1992	NA	NA	0.67
1993	NA	NA	0.73
1995	NA	NA	0.34
1996	NA	NA	1.00
1997	NA	NA	0.88
1998	NA	NA	0.28
2004	0.16	NA	NA
2005	0.68	NA	NA
2006	0.19	NA	NA
2007	0.31	NA	NA
2008	0.09	NA	NA
2009	0.52	0.75	0.82
2010	0.13	0.57	0.43
2011	0.27	0.32	0.37
2012	0.71	0.60	1.00
2013	0.36	0.75	0.85
2014	0.20	0.35	0.49
2015	0.29	0.40	0.58
2016	0.23	0.65	0.50
2017	0.12	0.15	0.51
2018	0.76	0.69	0.78
2019	0.35	0.37	0.36
2020	0.42	0.24	0.31

Рис. 2. Временная изменчивость в пространстве

Параллельные вычисления

cl <- makeCluster(16) # up to 16
doParallel::registerDoParallel(cl)
parallel::clusterEvalQ(cl, set.seed(10))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
## 
## [[7]]
## NULL
## 
## [[8]]
## NULL
## 
## [[9]]
## NULL
## 
## [[10]]
## NULL
## 
## [[11]]
## NULL
## 
## [[12]]
## NULL
## 
## [[13]]
## NULL
## 
## [[14]]
## NULL
## 
## [[15]]
## NULL
## 
## [[16]]
## NULL

r4 <- foreach(i=1:10000) %dopar% { ################ MUST BE 5000
    library(dplyr)
    if(sample(1:0, 1) == 1) {
        dfc <- wide.df %>% 
            group_by(km, period) %>% 
            sample_n(size = n()-1) %>% 
            ungroup
    } else { 
        dfc <- wide.df %>% 
            group_by(km, period) %>% 
            select(1:7, sort(sample(8:21, 13))) %>% 
            ungroup
    }
    bw.perkm <- dfc %>% 
        select(-id, -cycle, -zone, -total) %>% 
        mutate(nsp = nspec(., 3, ncol(.))) %>% 
        group_by(period, km) %>% 
        summarise_all(mean) %>% 
        ungroup() %>% 
        transmute(period, km, by.km.per = nspec(., 3, ncol(.) - 1) / nsp)
    
    dis <- as.matrix(vegan::vegdist(dfc[,8:ncol(dfc)], method = "bray", binary = FALSE))
    rownames(dis) <- colnames(dis) <- dfc$id
    dis <- reshape2::melt(dis, varnames = c("id1", "id2")) %>% 
        as_tibble() %>% 
        mutate(id1 = as.character(id1), id2 = as.character(id2), 
            value = case_when(is.nan(value) ~ 1, TRUE ~ value)) %>% 
        filter(as.numeric(substr(id1, 4, nchar(id1))) > as.numeric(substr(id2, 4, nchar(id2)))) %>% 
        left_join(select(dfc, id1 = id, km1 = km, year1 = year, per1 = period), by = "id1") %>% 
        left_join(select(dfc, id2 = id, km2 = km, year2 = year, per2 = period), by = "id2") %>% 
        filter(km1 == km2, per1 == per2) %>% 
        group_by(per1, km1) %>% 
        summarise(ibc = mean(value), .groups = "drop")
    list(bw.perkm = bw.perkm, ibc.perkm = dis)
}
stopCluster(cl) 
r4 <- transpose(r4)

Parallel computing is hidden

Бета Уиттекера

bwkm <- map_dfr(r4$bw.perkm, rbind) %>% 
    group_by(period, km) %>% 
    summarise(low = quantile(by.km.per, 0.025), up = quantile(by.km.per, 0.975), .groups = "drop")
bwkm <- wide.df  %>% 
    select(-id, -cycle, -zone, -total) %>% 
    mutate(nsp = nspec(., 3, ncol(.))) %>% 
    group_by(period, km) %>% 
    summarise_all(mean) %>% 
    ungroup() %>% 
    transmute(period, km, bw.perkm = nspec(., 3, ncol(.) - 1) / nsp, 
            zone = case_when(km == 20 | km == 30 ~ "1.fon", 
                                 km == 4 | km == 6 ~ "2.buf", 
                                 TRUE ~ "3.imp")) %>% 
    left_join(bwkm, by = c("period", "km")) %>% 
    mutate(km = as.factor(km))
plot2a <- ggplot() +
    geom_errorbar(aes(x = km, ymin = low, ymax = up), data = bwkm) + 
    facet_wrap(~period) + 
    geom_point(aes(x = km, y = bw.perkm, shape = zone, fill = zone), 
             data = bwkm, size = 4, color = "black") +
    theme_bw() +
    theme(legend.position = "bottom")+
    scale_shape_manual(values = c(21, 24, 22)) + 
    scale_fill_manual(values = c("white", "grey", "black")) + 
    labs(y = "Бета Уиттекера", x = "", title = "А")
plot2a

Расстояние Брея-Кёртиса

ibc.km <- as.matrix(vegan::vegdist(wide.df[,8:ncol(wide.df)], method = "bray", binary = FALSE))
rownames(ibc.km) <- colnames(ibc.km) <- wide.df$id
ibc.km <- reshape2::melt(ibc.km, varnames = c("id1", "id2")) %>% 
    as_tibble() %>% 
    mutate(id1 = as.character(id1), id2 = as.character(id2)) %>% 
    filter(as.numeric(substr(id1, 4, nchar(id1))) > as.numeric(substr(id2, 4, nchar(id2)))) %>% 
    left_join(select(wide.df, id1 = id, km1 = km, per1 = period), by = "id1") %>% 
    left_join(select(wide.df, id2 = id, km2 = km, per2 = period), by = "id2") %>% 
    filter(km1 == km2, per1 == per2) %>% 
    mutate(value = case_when(is.nan(value) ~ 1, TRUE ~ value)) %>% 
    group_by(per1, km1) %>% 
    summarise(ibc = mean(value), .groups = "drop") 
ibc.km <- map_dfr(r4$ibc.perkm, rbind) %>% 
    group_by(per1, km1) %>% 
    summarise(low = quantile(ibc, 0.025), up = quantile(ibc, 0.975), .groups = "drop") %>% 
    left_join(ibc.km, by = c("per1", "km1")) %>% 
    rename(period = per1, km = km1) %>% 
    mutate(zone = case_when(km == 20 | km == 30 ~ "1.fon", 
        km == 4 | km == 6 ~ "2.buf", TRUE ~ "3.imp"), km = as.factor(km))
plot2b <- ggplot(ibc.km, aes(x = km, ymin = low, ymax = up, y = ibc, shape = zone, fill = zone)) +
    geom_errorbar() + 
    facet_wrap(~period) + 
    geom_point(size = 4) +
    theme_bw() +
    theme(legend.position = "bottom")+
    scale_shape_manual(values = c(21, 24, 22)) + 
    scale_fill_manual(values = c("white", "grey", "black")) + 
    labs(y = "Расстояние Брея-Кёртиса", x = "", title = "Б")
plot2b

Вместе

plot2all <- gridExtra::grid.arrange(plot2a, plot2b, ncol = 1)

# ggsave(plot = plot2all, filename = "rpot2.svg", width = 7.5, height = 9, dpi = 500)

Таблицы с данными

Бета Уиттекера, все зоны (по км)

formattable::formattable(select(ibc.km, period, km, ibc, low, up))

period	km	ibc	low	up
1	1	0.8632193	0.8085071	0.9503054
1	1.5	0.7592627	0.6774281	0.8250385
1	2	0.6387517	0.5183356	0.7678596
1	4	0.5068144	0.4448033	0.7568949
1	20	0.6676890	0.5970366	0.8697692
2	1	1.0000000	1.0000000	1.0000000
2	2	0.7482978	0.7228334	0.8134757
2	4	0.6822947	0.6307425	0.9081259
2	6	0.6585416	0.5539124	0.7151281
2	20	0.6696325	0.6371717	0.7375415
2	30	0.5199612	0.4359269	0.6428391
3	1	0.7647557	0.7124792	0.8731461
3	1.5	0.6212653	0.5371020	0.8089971
3	2	0.6152604	0.5297627	0.6511490
3	4	0.5379929	0.4907400	0.5891965
3	6	0.6034536	0.5797945	0.6957623
3	20	0.5571585	0.4720257	0.7127891
3	30	0.5349196	0.4730093	0.7536196

Расстояние Брея-Кёртиса, все зоны

formattable::formattable(select(bwkm, period, km, bw.perkm, low, up))

period	km	bw.perkm	low	up
1	1	2.520000	2.250000	2.842105
1	1.5	1.923077	1.555556	2.142857
1	2	2.057143	1.814815	2.285714
1	4	1.920000	1.787234	2.095238
1	20	1.833333	1.673913	2.000000
2	1	1.600000	1.000000	1.600000
2	1.5	1.000000	1.000000	1.000000
2	2	2.051282	1.741935	2.187500
2	4	2.105263	1.750000	2.500000
2	6	1.935484	1.521739	2.076923
2	20	1.967213	1.800000	2.156863
2	30	1.463415	1.350000	1.562500
3	1	2.357143	2.000000	2.588235
3	1.5	1.650000	1.410256	1.833333
3	2	1.604167	1.395349	1.736842
3	4	1.718750	1.571429	1.867925
3	6	2.037037	1.867925	2.250000
3	20	2.062500	1.920635	2.283019
3	30	1.958904	1.791045	2.129032

Рис. 4. Попарные сравнения

Для соблюдения трехлетнести циклов, 2007 и 2008 усреднены. Другие ключевые решения:

Рис. 4a. Данные - обилие. Все расстояния до нулевых годов реконструированы пакетом zoo. Рис. 4b. Данные - доли. Все расстояния до нулевых годов реконструированы пакетом zoo. Рис. 4c. Данные - обилие. Все расстояния до нулевых годов приняты за 1 Рис. 4d. Данные - доли. Все расстояния до нулевых годов приняты за 1.

Функции для расчётов

my.dis <- function(df1, start = 5) { 
    df1[,start:ncol(df1)] %>% vegan::vegdist() %>% 
    as.matrix() %>% 
    as_tibble() %>% 
    rownames_to_column("r1") %>% 
    pivot_longer(names_to = "r2", values_to = "dis", -r1) %>% 
    mutate(r1 = as.numeric(r1), r2 = as.numeric(r2)) %>% 
    filter(r1 < r2) %>%  
    mutate(dis = case_when(is.na(dis) ~ 1, TRUE ~ dis)) %>% 
    
    left_join(transmute(rownames_to_column(df1, "r1"), 
        r1 = as.numeric(r1), zone1 = zone, year1 = year, km1 = km), by = "r1") %>%  
    left_join(transmute(rownames_to_column(df1, "r2"), 
        r2 = as.numeric(r2), zone2 = zone, year2 = year, km2 = km), by = "r2") %>% 
    filter(year1 == year2, zone1 != zone2) %>% 
    transmute(pair = paste0(zone1, "_", zone2), year = year1, dis) %>% 
    group_by(year, pair) %>% 
    summarise(dis = mean(dis), .groups = "drop") %>% 
    pivot_wider(names_from = pair, values_from = dis) %>% 
    rbind(c(1999, NA, NA, NA), c(2000, NA, NA, NA)) %>% 
    arrange(year)
}
my.ts <- function(df2) { 
    df1 <- df2 %>% select(-year) %>% 
        lapply(function(a){
            ts(a, frequency = 3) %>% 
                zoo::na.StructTS(.) %>% 
                decompose(type = "additive")
            }) %>% 
    map(.f = function(a){keep(a, is.ts)}) %>% 
    map(.f = function(a){map_dfr(a, cbind)})
    df1 %>% map_dfr(rbind) %>% 
        cbind(pair = rep(names(df1), each = nrow(df1[[1]])), 
             year = c(1990:2006, 2007.5, 2009:2020)) %>% 
        as_tibble() %>% # ^ year = 1990:2020) 
        pivot_longer(names_to = "type", values_to = "val", -c("pair", "year")) %>% 
        mutate(type = case_when(type == "x" ~ "1.x", type == "trend" ~ "2.trend",
                            type == "seasonal" ~ "3.cycle", type == "random" ~ "4.random"))
}
my.plot <- function(tmp) { 
    tmp %>% 
        ggplot(aes(x = year, y = val, color = type)) + 
        geom_line() + 
        facet_grid(cols = vars(pair), rows = vars(type), scales = "free") + 
        geom_point() + 
        labs(x = "", y = "")+
        theme_bw() + 
        theme(legend.position = "none")
}
wide78 <- wide.df %>% select(-id, -cycle, -period)
wide78 <- wide78 %>% filter(year == 2007 | year == 2008) %>% 
    group_by(zone, km) %>% 
    summarise_all(mean) %>% 
    select(year, zone, km, total, 5:18) %>% 
    rbind(filter(wide78, year != 2007 & year != 2008)) %>%
    arrange(year, km)

функции скрыты

4a. Отн.обилие, 0&NA <- zoo

df1 <- wide78 %>% 
    my.dis() %>% 
    my.ts() 
my.plot(df1)

Параметры моделей (линейных):

4a.1 - buf_fon

df1 %>% filter(pair == "2.buf_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14198 -0.06736 -0.01497  0.07323  0.13786 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -11.591947   3.689769  -3.142  0.00416 **
## year          0.006027   0.001840   3.275  0.00299 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08278 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.292,  Adjusted R-squared:  0.2648 
## F-statistic: 10.72 on 1 and 26 DF,  p-value: 0.002992

4a.2 - imp_fon

df1 %>% filter(pair == "3.imp_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.178605 -0.039420  0.005527  0.047803  0.124151 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -14.688791   3.240633  -4.533 0.000115 ***
## year          0.007706   0.001616   4.768 6.21e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0727 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.4664, Adjusted R-squared:  0.4459 
## F-statistic: 22.73 on 1 and 26 DF,  p-value: 6.208e-05

4a.3 - imp_buf

df1 %>% filter(pair == "3.imp_2.buf") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.172507 -0.076046  0.000943  0.064194  0.152857 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -11.399553   4.269339  -2.670   0.0129 * 
## year          0.005992   0.002129   2.814   0.0092 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09578 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2334, Adjusted R-squared:  0.2039 
## F-statistic: 7.917 on 1 and 26 DF,  p-value: 0.009205

4b. Доли, 0&NA <- zoo

df2 <- wide78[,5:18] %>% 
    mutate_all(.funs = function(a) {a/wide78$total}) %>% 
    as.matrix()
df2[is.nan(df2)] <- 0
df2 <- cbind(wide78[,1:4], df2) %>% 
    my.dis() %>% 
    my.ts() 
my.plot(df2)

Параметры моделей (линейных):

4b.1 - buf_fon

df2 %>% filter(pair == "2.buf_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.10442 -0.06682 -0.01334  0.06836  0.13681 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.941378   3.515782  -5.103 2.56e-05 ***
## year          0.009135   0.001754   5.209 1.94e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07887 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.5107, Adjusted R-squared:  0.4919 
## F-statistic: 27.14 on 1 and 26 DF,  p-value: 1.937e-05

4b.2 - imp_fon

df2 %>% filter(pair == "3.imp_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.238411 -0.053487 -0.004448  0.087130  0.178259 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -25.983430   4.963780  -5.235 1.81e-05 ***
## year          0.013271   0.002476   5.360 1.30e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1114 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.525,  Adjusted R-squared:  0.5067 
## F-statistic: 28.73 on 1 and 26 DF,  p-value: 1.301e-05

4b.3 - imp_buf

df2 %>% filter(pair == "3.imp_2.buf") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21596 -0.08742 -0.02471  0.09886  0.18089 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.763393   5.375646  -1.444    0.161
## year         0.004139   0.002681   1.544    0.135
## 
## Residual standard error: 0.1206 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.08394,    Adjusted R-squared:  0.04871 
## F-statistic: 2.382 on 1 and 26 DF,  p-value: 0.1348

4c. Отн.обилие, 99/01 как есть

df3 <- wide78 %>% 
    my.dis() 
df3[is.na(df3)] <- 1
df3 <- my.ts(df3)
my.plot(df3)

Параметры моделей (линейных):

4c.1 buf_fon

df3 %>% filter(pair == "2.buf_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21266 -0.08065 -0.01203  0.07802  0.23660 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.226008   4.958913  -1.256    0.220
## year         0.003370   0.002473   1.362    0.185
## 
## Residual standard error: 0.1112 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.06663,    Adjusted R-squared:  0.03073 
## F-statistic: 1.856 on 1 and 26 DF,  p-value: 0.1848

4c.2 imp_fon

df3 %>% filter(pair == "3.imp_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.214155 -0.057909 -0.008334  0.048255  0.249788 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -11.863043   4.379193  -2.709  0.01178 * 
## year          0.006307   0.002184   2.887  0.00772 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09824 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2428, Adjusted R-squared:  0.2137 
## F-statistic: 8.337 on 1 and 26 DF,  p-value: 0.007722

4c.3 imp_buf

df3 %>% filter(pair == "3.imp_2.buf") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.23071 -0.08958  0.00390  0.07321  0.37163 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.868454   5.885333  -1.167    0.254
## year         0.003748   0.002935   1.277    0.213
## 
## Residual standard error: 0.132 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.05901,    Adjusted R-squared:  0.02282 
## F-statistic: 1.631 on 1 and 26 DF,  p-value: 0.2129

4d. Доли, 99/01 как есть

df4 <- wide78[,5:18] %>% 
    mutate_all(.funs = function(a) {a/wide78$total}) %>% 
    as.matrix()
df4[is.nan(df4)] <- 0
df4 <- cbind(wide78[,1:4], df4) %>% 
    my.dis() 
df4[is.na(df4)] <- 1
df4 <- df4 %>% my.ts() 
my.plot(df4)

Параметры моделей (линейных):

4d.1 buf_fon

df4 %>% filter(pair == "2.buf_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.18417 -0.10433 -0.05255  0.06192  0.36323 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -11.928932   6.688871  -1.783   0.0862 .
## year          0.006158   0.003336   1.846   0.0763 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1501 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.1159, Adjusted R-squared:  0.08185 
## F-statistic: 3.407 on 1 and 26 DF,  p-value: 0.07634

4d.2 imp_fon

df4 %>% filter(pair == "3.imp_1.fon") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27285 -0.06442 -0.03112  0.06627  0.39741 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -21.569061   6.454023  -3.342  0.00253 **
## year          0.011086   0.003219   3.444  0.00196 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1448 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.3133, Adjusted R-squared:  0.2868 
## F-statistic: 11.86 on 1 and 26 DF,  p-value: 0.001956

4d.3 imp_buf

df4 %>% filter(pair == "3.imp_2.buf") %>% 
    pivot_wider(names_from = type, values_from = "val") %>% 
    lm(`2.trend` ~ year, data = .) %>% 
    summary

## 
## Call:
## lm(formula = `2.trend` ~ year, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.27850 -0.08891 -0.05250  0.07726  0.43825 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.895289   7.244263  -0.400    0.693
## year         0.001729   0.003613   0.478    0.636
## 
## Residual standard error: 0.1625 on 26 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.008725,   Adjusted R-squared:  -0.0294 
## F-statistic: 0.2289 on 1 and 26 DF,  p-value: 0.6364

все вместе

gridExtra::grid.arrange(my.plot(df1) + labs(title = "1"), my.plot(df2) + labs(title = "2"), my.plot(df3) + labs(title = "3"), my.plot(df4) + labs(title = "4"))

#ggsave(plot = p4b, filename = "rpot4b.svg", width = 9, height = 7, dpi = 500)

Beta-diversity, v 7

А. Созонтов

21 02 2021

Загрузка данных

Рис. 1. Ординация

Обилие*цикл

Доли * цикл

Для публикации

Табл. 1. PERMANOVA

Доли*год

Рис. 3. Пространственная изменчивость во времени

Параллельные вычисления

Расчёт показателей

Бета Уиттекера

plot3a

plot3a + trends

Bw: линейная модель

Bw: полиномиальная модель

Расстояние Брея-Кёртиса

plot3b

plot3b + trends

Ibc: линейная модель

Ibc: полиномиальная модель

Всё вместе

Таблицы с данными

Бета Уиттекера, всесь градиент

Бета Уиттекера, по зонам

Расстояние Брея-Кёртиса, всесь градиент

Расстояние Брея-Кёртиса, по зонам

Рис. 2. Временная изменчивость в пространстве

Параллельные вычисления

Бета Уиттекера

Расстояние Брея-Кёртиса

Вместе

Таблицы с данными

Бета Уиттекера, все зоны (по км)

Расстояние Брея-Кёртиса, все зоны

Рис. 4. Попарные сравнения

Функции для расчётов

4a. Отн.обилие, 0&NA <- zoo

4a.1 - buf_fon

4a.2 - imp_fon

4a.3 - imp_buf

4b. Доли, 0&NA <- zoo

4b.1 - buf_fon

4b.2 - imp_fon

4b.3 - imp_buf

4c. Отн.обилие, 99/01 как есть

4c.1 buf_fon

4c.2 imp_fon

4c.3 imp_buf

4d. Доли, 99/01 как есть

4d.1 buf_fon

4d.2 imp_fon

4d.3 imp_buf

все вместе