Joint Species

Fit

Code

data <- read_tsv("../modelling/data/derived_data/data_rec.tsv", col_types = cols()) %>%
  select(-y) %>% 
  filter(variable %in% c("ba", "diversity_q1_sp")) %>%
  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)
trajectories <- read_tsv("chains/joint-diversity_q1_sp-trajectories.tsv")
parameters <-  read_tsv("chains/joint-diversity_q1_sp-parameters.tsv")
parameters

# A tibble: 727 × 14
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__     16454.      1.65e+4 121.    135.     1.62e+4  1.66e+4  2.34     4.98
 2 gamma        0.0483  4.88e-2   0.218   0.204 -3.05e-1  4.03e-1  1.00  2068.  
 3 gamma       -0.146  -1.59e-1   0.166   0.163 -4.00e-1  1.41e-1  1.07    44.9 
 4 gamma       -0.454  -4.53e-1   0.252   0.243 -8.91e-1 -5.08e-2  1.02  1016.  
 5 gamma       -0.131  -1.11e-1   0.419   0.384 -8.29e-1  4.97e-1  1.01  1733.  
 6 gamma       -0.0785 -8.15e-2   0.302   0.289 -5.87e-1  4.32e-1  1.02   524.  
 7 gamma        0.145   1.55e-1   0.336   0.334 -4.22e-1  6.87e-1  1.02   218.  
 8 gamma        0.188   1.72e-1   0.211   0.195 -1.21e-1  5.66e-1  1.01  1215.  
 9 gamma       -0.0867 -5.45e-2   0.316   0.381 -6.07e-1  3.67e-1  1.65     6.53
10 gamma       -0.217  -2.13e-1   0.173   0.170 -5.03e-1  6.41e-2  1.01  1787.  
# ℹ 717 more rows
# ℹ 5 more variables: ess_tail <dbl>, sitenum <dbl>, site <chr>, plotnum <dbl>,
#   plot <chr>

IDH - H1

Code

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

Diveristy hump delta with basal area disturbance intensity phi for each plot in each site.

IDH - H2

Code

area <- data %>% 
  filter(variable == "area") %>% 
  select(site, plot, y) %>% 
  unique() %>% 
  group_by(site) %>% 
  summarise(area = mean(y))
filter(parameters, variable == "gamma") %>% 
  left_join(area) %>% 
  ggplot(aes(median, reorder(site, median))) +
  geom_vline(xintercept = 0, col = "darkgrey") +
  geom_segment(aes(x = (median - sd),
                   xend = (median + sd))) +
  geom_segment(aes(x = q5, xend = q95), size = .1) +
  geom_point() +
  geom_text(aes(label = paste(round(area, 1), "ha")), nudge_y = .5, size = 2) +
  theme_bw() +
  scale_x_reverse(expression(gamma~"["~"%"~"%"^{-1}~"]")) +
  ylab("") +
  scale_color_viridis_c()

Posteriors of the parameter gamma representing basal are disturbance intensity effect on genus diversity short term increase with plot areas in label.

Comparison

Code

read_tsv("chains/joint-diversity_q1_gen-parameters.tsv") %>% 
  filter(variable == "gamma") %>% 
  left_join(filter(parameters, variable == "gamma") %>% 
              select(site, median, sd, q5, q95) %>% 
              rename(median_new = median, sd_new = sd,
                     q5_new = q5, q95_new = q95)
  ) %>% 
  filter(!is.na(median_new)) %>% 
  ggplot(aes(median, median_new)) +
  geom_hline(yintercept = 0, col = "grey") +
  geom_vline(xintercept = 0, col = "grey") +
  geom_point() +
  geom_segment(aes(x = median - sd, xend = median + sd), alpha = .5) +
   geom_segment(aes(y = median_new - sd_new, yend = median_new + sd_new), alpha = .5) +
  ggrepel::geom_text_repel(aes(label = site)) +
  theme_bw() +
  geom_abline() +
  xlab("gamma from genus 1") +
  ylab("gamma from species 1") +
  scale_x_reverse() +
  scale_y_reverse()

Comparison with joint inference of genus order 1 diversity and basal area.