Joint Genus 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_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)
trajectories <- read_tsv("chains/joint-diversity_q0_gen-trajectories.tsv")
parameters <-  read_tsv("chains/joint-diversity_q0_gen-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__        1.78e+4  1.78e+4 1.11e+2 1.17e+2  1.76e+4  1.80e+4  1.93     5.67
 2 gamma      -3.40e-2 -3.32e-2 1.73e-1 1.65e-1 -3.12e-1  2.46e-1  1.01   718.  
 3 gamma       1.03e-1  1.02e-1 8.88e-2 8.58e-2 -4.22e-2  2.51e-1  1.02   234.  
 4 gamma      -2.95e-1 -2.95e-1 1.16e-1 1.13e-1 -4.85e-1 -1.01e-1  1.01  1102.  
 5 gamma       3.82e-2  3.44e-2 2.53e-1 2.35e-1 -3.70e-1  4.45e-1  1.01  2257.  
 6 gamma      -2.32e-1  1.34e-1 9.30e-1 3.50e-1 -2.43e+0  5.24e-1  1.70     6.26
 7 gamma      -1.88e-1 -1.84e-1 3.36e-1 3.08e-1 -7.56e-1  3.49e-1  1.02   363.  
 8 gamma       2.15e-1  2.04e-1 1.46e-1 1.40e-1 -2.97e-4  4.67e-1  1.01  1086.  
 9 gamma      -1.00e-1 -1.02e-1 9.35e-2 9.25e-2 -2.51e-1  5.35e-2  1.00  2626.  
10 gamma      -2.18e-3  1.42e-3 1.29e-1 1.20e-1 -2.22e-1  2.05e-1  1.00  1735.  
# ℹ 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 genus 0") +
  scale_x_reverse() +
  scale_y_reverse()

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