\[
\delta_{plot} \sim \mathcal N ( \delta0_{site}+\phi_{BA,plot} \times \gamma_{site}+\phi_{BA,plot}^2 \times \iota_{site}, \sigma_{\delta})
\]
Fit
Code
data <- read_tsv("../modelling/data/derived_data/data_rec.tsv", col_types = cols()) %>%
select(-y) %>%
filter(variable %in% c("ba", "diversity_q1_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_quad-diversity_q1_gen-trajectories.tsv")
parameters <- read_tsv("chains/joint_quad-diversity_q1_gen-parameters.tsv")
parameters
# A tibble: 745 × 14
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ 16844. 16844. 60.1 62.5 16747. 16945. 1.29 11.7
2 gamma -1.20 -1.43 1.82 1.95 -3.72 2.22 1.02 292.
3 gamma 0.502 0.594 1.76 1.90 -2.58 3.27 1.05 197.
4 gamma 1.21 1.50 1.86 1.98 -2.39 3.72 1.07 53.4
5 gamma 0.116 0.174 2.11 2.54 -3.37 3.45 1.01 386.
6 gamma -0.654 -0.716 0.981 0.831 -2.11 1.11 1.18 19.1
7 gamma 0.234 0.271 2.33 3.00 -3.52 3.73 1.23 14.0
8 gamma 0.702 1.08 2.08 2.17 -3.25 3.54 1.23 13.6
9 gamma -1.71 -2.10 1.78 1.81 -3.88 1.81 1.17 18.7
10 gamma 0.829 0.811 1.31 1.32 -1.34 3.03 1.04 232.
# ℹ 735 more rows
# ℹ 5 more variables: ess_tail <dbl>, sitenum <dbl>, site <chr>, plotnum <dbl>,
# plot <chr>
Trajectories
Code
trajectories %>%
filter(variable == "diversity_q1_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_q1_gen"), alpha = .5
) +
facet_wrap(~site, scales = "free") +
theme_bw() +
xlab("") + ylab("") +
theme(legend.position = "bottom") +
scale_color_discrete(guide = "none")
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()
Code
quad <- filter(parameters, variable %in% c("delta0", "gamma", "iota")) %>%
select(site, variable, median) %>%
pivot_wider(names_from = variable, values_from = median) %>%
mutate(phiba = list(0:100/100)) %>%
unnest(phiba) %>%
mutate(delta = delta0 + gamma*phiba + iota*phiba^2)
filter(parameters, variable %in% c("delta", "phiba")) %>%
select(site, plot, variable, median) %>%
pivot_wider(names_from = variable, values_from = median) %>%
ggplot(aes((1-phiba), delta)) +
facet_wrap(~ site) +
geom_point() +
geom_line(data = quad) +
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)
