Joint Genus 0

Code

recovery <- cmdstan_model("models/recovery.stan")
data <- read_tsv("../modelling/data/derived_data/data_rec.tsv", col_types = cols()) %>%
    select(-y) %>% 
    filter(variable %in% c("ba", "diversity_q0_gen")) %>%
    select(-sitenum, -plotnum) %>% 
    pivot_wider(names_from = variable, values_from = value) %>% 
    na.omit() %>% 
    filter(!(site %in% c("Manare", "Montagne Tortue", "Nelliyampathy",
                         "Uppangala", "BAFOG", "Sao Nicolau",
                         "Kabo", "Mil", "Corinto",
                         "Peixoto", "Iwokrama", "Antimary", "Peteco"))) %>% 
    mutate(sitenum = as.numeric(as.factor(site))) %>% 
    mutate(plotnum = as.numeric(as.factor(paste(site, plot)))) %>% 
    arrange(sitenum, plotnum, rel_year)
ind <- data %>%
    select(site, plot, sitenum, plotnum) %>%
    unique()
ba_priors <- read_tsv(file.path("..", "modelling",
                                "chains", paste0("equ_", "ba"), "mu.tsv")) %>% 
    select(-sitenum) %>% 
    left_join(select(ind, site, sitenum) %>% unique()) %>% 
    na.omit()
div_priors <- read_tsv(file.path("..", "modelling",
                                 "chains", paste0("equ_", "diversity_q0_gen"),
                                 "mu.tsv")) %>% 
    select(-sitenum) %>% 
    left_join(select(ind, site, sitenum) %>% unique()) %>% 
    na.omit()
mdata <- list(
    n = nrow(data),
    s = max(data$sitenum),
    p = max(data$plotnum),
    ba = data$ba,
    y = data$diversity_q0_gen,
    t = data$rel_year - 3,
    site = data$sitenum,
    plot = data$plotnum,
    site_plot = ind$sitenum,
    mu_ba_theta_s = ba_priors$median,
    sigma_ba_theta_s = ba_priors$sd,
    mu_theta_s = div_priors$median,
    sigma_theta_s = div_priors$sd
)
fit <- recovery$sample(
    mdata,
    chains = 4,
    parallel_chains = 4,
    save_warmup = FALSE,
    refresh = 100,
    output_dir = "chains/jointgen0/"
)

Fit

Code

fit <- as_cmdstan_fit(list.files("chains/jointgen0/", full.names = TRUE))
data <- read_tsv("../modelling/data/derived_data/data_rec.tsv", col_types = cols()) %>%
  select(-y) %>% 
  filter(variable %in% c("ba", "diversity_q0_gen")) %>%
  select(-sitenum, -plotnum) %>% 
  pivot_wider(names_from = variable, values_from = value) %>% 
  na.omit() %>% 
filter(!(site %in% c("Manare", "Montagne Tortue", "Nelliyampathy",
                         "Uppangala", "BAFOG", "Sao Nicolau",
                         "Kabo", "Mil", "Corinto",
                         "Peixoto", "Iwokrama", "Antimary", "Peteco"))) %>% 
  mutate(sitenum = as.numeric(as.factor(site))) %>% 
  mutate(plotnum = as.numeric(as.factor(paste(site, plot)))) %>% 
  arrange(sitenum, plotnum, rel_year)  %>% 
  gather(variable, y, -plot, -year, -site, -treatment, -harvest_year,
         -longitude, -harvest_year_min, -rel_year, -sitenum, -plotnum)
ind <- data %>%
    select(site, plot, sitenum, plotnum) %>%
    unique()
trajs <- fit$summary(c("y_pred", "ba_pred"), "median") %>% 
  separate(variable,
      c("variable", "pred", "rel_year", "plotnum"),
      convert = TRUE
  ) %>%
  select(-pred) %>% 
  mutate(variable = recode(variable, "y" = "diversity_q0_gen")) %>% 
  rename(y = median) %>%
  mutate(rel_year = rel_year + 3) %>%
  left_join(ind)
fit

    variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
 lp__        9568.72 9564.72 92.86 94.55 9422.53 9727.89 1.20       16       64
 phi_ba_p[1]    0.74    0.74  0.04  0.04    0.67    0.80 1.01     1245     2390
 phi_ba_p[2]    0.84    0.84  0.04  0.04    0.76    0.91 1.00     1216     2565
 phi_ba_p[3]    0.96    0.96  0.03  0.03    0.89    1.00 1.00     1527     1929
 phi_ba_p[4]    0.76    0.76  0.04  0.04    0.69    0.83 1.00     1327     2072
 phi_ba_p[5]    0.73    0.73  0.04  0.04    0.66    0.80 1.00     1326     2408
 phi_ba_p[6]    0.87    0.87  0.05  0.05    0.79    0.94 1.01      841     1308
 phi_ba_p[7]    0.92    0.92  0.04  0.04    0.84    0.98 1.01     1021     1197
 phi_ba_p[8]    0.87    0.87  0.04  0.04    0.80    0.94 1.01     1010     1307
 phi_ba_p[9]    0.76    0.76  0.04  0.04    0.69    0.83 1.00     1104     1608

 # showing 10 of 49672 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

Trajectories

Code

filter(trajs, variable == "ba") %>% 
  ggplot(aes(rel_year, y)) +
  geom_line(aes(group = paste(site, plotnum), col = as.character(plotnum))) +
  geom_point(aes(group = paste(site, plotnum), col = as.character(plotnum)),
      data = data %>% filter(variable == "ba"), alpha = .5
  ) +
  facet_wrap(~site, scales = "free") +
  theme_bw() +
  xlab("") + ylab("") +
  theme(legend.position = "bottom") +
  scale_color_discrete(guide = "none")

Basal area jointly inferred trajectories for all sites.

Code

filter(trajs, variable == "diversity_q0_gen") %>% 
  ggplot(aes(rel_year, y)) +
  geom_line(aes(group = paste(site, plotnum), col = as.character(plotnum))) +
  geom_point(aes(group = paste(site, plotnum), col = as.character(plotnum)),
      data = data %>% filter(variable == "diversity_q0_gen"), alpha = .5
  ) +
  facet_wrap(~site, scales = "free") +
  theme_bw() +
  xlab("") + ylab("") +
  theme(legend.position = "bottom") +
  scale_color_discrete(guide = "none")

Genus diversity jointly inferred trajectories for all sites.

Disturbance

Code

dist_obs <- data %>% 
  filter(variable == "ba") %>% 
  group_by(site, plot) %>%
  summarise(ba_post = min(y, na.rm = TRUE)) %>% 
  left_join(read_tsv("../modelling/data/derived_data/data_equ.tsv",
                     col_types = cols()) %>%
              filter(variable == "ba") %>% 
              filter(rel_year <= 0) %>%
              group_by(site, plot) %>%
              summarise(ba_pre = median(value, na.rm = TRUE))) %>% 
    na.omit() %>% 
  mutate(dist_obs = (ba_pre - ba_post) / ba_pre * 100) %>%
  select(site, plot, dist_obs) %>%
  filter(dist_obs > 0)
fit$summary("phi_ba_p") %>% 
  separate(variable, c("X1", "X2", "X3", "plotnum"), convert = TRUE) %>% 
  select(-X1, -X2, -X3) %>% 
  left_join(ind) %>% 
  left_join(dist_obs) %>%
  na.omit() %>%
  ggplot(aes((1-median)*100, dist_obs)) +
  geom_abline(col = "darkgrey") +
  geom_segment(aes(x = median - sd, xend = median + sd)) +
  geom_segment(aes(x = q5, xend = q95), size = .1) +
  geom_point() +
  facet_wrap(~site) +
  theme_bw() +
  ggpubr::stat_cor(size = 3) +
  xlab(expression(phi[BA] ~ "inferred [ % lost ]")) +
  ylab(expression(dist[BA] ~ "observed [ % lost ]")) +
  theme(panel.spacing = unit(0, "lines")) +
  xlim(0, 100) + ylim(0, 100)

Evaluation of inferred disturbance intensity against observed with pre-logging inventories.

IDH - H1

Code

pars <- ind %>% 
  left_join(fit$summary("delta_p") %>% 
  separate(variable, c("parameter", "p", "plotnum"), convert = TRUE) %>% 
              select(-p)) %>% 
  left_join(fit$summary("phi_ba_p") %>% 
              mutate(variable = gsub("phi_ba_p", "phiba_p", variable)) %>% 
              separate(variable, c("parameter", "p", "plotnum"), convert = TRUE) %>% 
              select(-p), by = "plotnum")

pars %>% 
  ggplot(aes(median.x*100, (1-median.y)*100)) +
  geom_vline(xintercept = 0) +
  geom_segment(aes(x = (median.x - sd.x)*100,
                   xend = (median.x + sd.x)*100), col = "darkgrey") +
  geom_segment(aes(x = q5.x*100, xend = q95.x*100), size = .1, col = "darkgrey") +
  geom_segment(aes(y = (1-(median.y - sd.y))*100,
                   yend = (1-(median.y + sd.y))*100), col = "darkgrey") +
    geom_point() +
  facet_wrap(~ site, scales = "free") +
  theme_bw() +
  ylab(expression(phi[BA] ~ "[ % lost ]")) +
  xlab(expression(delta[Diversity] ~ "[ % ]"))

Code

pars %>% 
  ggplot(aes(median.x)) +
  geom_histogram(fill = "lightgrey", col = "black") +
  geom_vline(xintercept = 0) +
  theme_bw() +
  facet_wrap(~ site, scales = "free_y") +
  xlab(expression(delta[Diversity] ~ "[ % ]"))

Code

pars %>% 
  ggplot(aes(median.x)) +
  geom_histogram(aes(fill = site), col = "black") +
  geom_vline(xintercept = 0) +
  theme_bw() +
  xlab(expression(delta[Diversity] ~ "[ % ]"))

Code

pars %>% 
  filter(!grepl("Tapajos", site)) %>% 
  ggplot(aes(median.x)) +
  geom_histogram(aes(fill = site), col = "black") +
  geom_vline(xintercept = 0) +
  theme_bw() +
  xlab(expression(delta[Diversity] ~ "[ % ]"))

IDH - H2

Code

fit$draws(c("gamma_s")) %>% 
  mcmc_intervals() +
  scale_y_discrete(labels = unique(ind$site)) +
  scale_x_reverse(expression(gamma~"["~"%"~"%"^{-1}~"]"), lim = c(4,-4)) 

Posterior of the parameter gamma representing basal are disturbance intensity effect on genus diversity short term increase

Code

pars <- ind %>% 
  left_join(fit$summary(c("delta0_s"), "median") %>% 
              separate(variable, c("parameter", "s", "sitenum"), convert = TRUE) %>% 
              select(-s) %>% 
              pivot_wider(names_from = parameter, values_from = median)) %>% 
  left_join(fit$summary("gamma_s") %>% 
              separate(variable, c("parameter", "s", "sitenum"), convert = TRUE) %>% 
              select(-s) %>% 
              pivot_wider(names_from = parameter, values_from = median)) %>% 
  left_join(fit$summary(c("phi_ba_p", "delta_p"), "median") %>% 
              mutate(variable = gsub("phi_ba_p", "phiba_p", variable)) %>% 
              separate(variable, c("parameter", "p", "plotnum"), convert = TRUE) %>% 
              select(-p) %>% 
              pivot_wider(names_from = parameter, values_from = median))
pars %>% 
  ggplot(aes((1-phiba), delta)) +
  facet_wrap(~ site) +
  geom_abline(aes(slope = -mean, 
                  intercept = mean+delta0, group = site), 
              data = pars) +
    geom_abline(aes(slope = -mean-sd, 
                  intercept = mean-sd+delta0, group = site), 
              data = pars, linetype = "dashed") +
  geom_abline(aes(slope = -mean+sd, 
                  intercept = mean+sd+delta0, group = site), 
              data = pars, linetype = "dashed") +
  geom_point() +
  theme_bw() +
  xlab(expression(phi[BA] ~ "[ % lost ]")) +
  ylab(expression(delta[Diversity] ~ "[ % ]")) +
  theme(panel.spacing = unit(0, "lines")) +
  scale_x_continuous(labels = ~.*100) +
  scale_y_continuous(labels = ~.*100)

Inferred relationship between disturbance intensity and diversity short term increase.

Code

fit$summary(c("gamma_s")) %>% 
  ggplot(aes(median)) +
  geom_histogram(fill = "lightgrey", col = "black", binwidth = 0.2) +
  geom_vline(xintercept = 0) +
  theme_bw() +
  scale_x_reverse(expression(gamma~"["~"%"~"%"^{-1}~"]"))