Joint Species 0

Fit

Code

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

# A tibble: 709 × 14
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__     16799.      1.68e+4 108.    103.     1.66e+4  1.70e+4  1.97     5.52
 2 gamma       -0.285  -2.80e-1   0.191   0.182 -5.94e-1  2.88e-2  1.00  1551.  
 3 gamma       -0.104  -9.94e-2   0.123   0.122 -3.11e-1  8.83e-2  1.01   927.  
 4 gamma       -0.314  -3.19e-1   0.168   0.164 -5.81e-1 -3.06e-2  1.02   255.  
 5 gamma       -0.266  -2.55e-1   0.291   0.272 -7.45e-1  1.97e-1  1.00  2126.  
 6 gamma       -0.0332 -2.45e-2   0.189   0.188 -3.53e-1  2.63e-1  1.01   963.  
 7 gamma       -0.640  -6.97e-1   0.959   1.18  -1.91e+0  9.66e-1  1.73     6.22
 8 gamma        0.103   9.63e-2   0.187   0.188 -1.98e-1  4.06e-1  1.12    24.3 
 9 gamma       -0.179  -1.88e-1   0.154   0.155 -4.15e-1  8.24e-2  1.23    12.6 
10 gamma       -0.220  -2.17e-1   0.152   0.150 -4.73e-1  2.73e-2  1.00  2078.  
# ℹ 699 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 0") +
  scale_x_reverse() +
  scale_y_reverse()

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