[1] "/home/runner/work/gcrmn_model_dev/gcrmn_model_dev/docs"Comparing model types - missing years
1 Purpose
2 Preparations
Load the necessary R libraries
library(tidyverse) # for data manipulation and visualisation
library(sf) # for spatial data handling and visualisation
library(glmmTMB) # for frequentist GLMM
library(emmeans) # for post-hoc analysis
library(DHARMa) # for model diagnostics
library(patchwork)
library(brms) # for Bayesian GLMM
library(rstan) # for Bayesian modelling
library(bayesplot) # for Bayesian plotting
library(tidybayes) # for tidy bayesian analysis
library(posterior) # for posterior analysis
library(gbm) # for gradient boosted regression trees
library(dbarts) # for Bayesian additive regression trees
library(HDInterval) # for Bayesian credible intervalsEstablish the global environment and paths
assign(x = "data_path", value = "../data/", envir = .GlobalEnv)
assign(x = "output_path", value = "../output/", envir = .GlobalEnv)
paths <- list(
data_path = data_path,
synthetic_path = paste0(data_path, "synthetic/"),
output_path = output_path,
fig_path = paste0(output_path, "figures/")
)
lapply(paths, function(x) {
if (!dir.exists(x)) {
dir.create(x)
}
})Load any helper functions
3 Import data
These data represent the synthesized “truth” and are thereby useful to compare to the various modelled outcomes.
Simple feature collection with 267132 features and 8 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 0.1086254 ymin: -20.96124 xmax: 9.508625 ymax: -9.981244
Geodetic CRS: WGS 84
# A tibble: 267,132 × 9
Year Reef geometry HCC SC MA CYC DHW OTHER
<int> <chr> <POINT [°]> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 Reef1 (9.398625 -19.95124) 0.310 0.0750 0.565 0 0.399 0.512
2 2 Reef1 (9.398625 -19.95124) 0.270 0.0627 0.592 0.570 0.379 0.745
3 3 Reef1 (9.398625 -19.95124) 0.279 0.0652 0.586 0 0.417 0.486
4 4 Reef1 (9.398625 -19.95124) 0.291 0.0691 0.578 0.0180 0.340 0.615
5 5 Reef1 (9.398625 -19.95124) 0.301 0.0722 0.571 0 0.419 0.429
6 6 Reef1 (9.398625 -19.95124) 0.316 0.0769 0.561 0 0.318 0.728
7 7 Reef1 (9.398625 -19.95124) 0.328 0.0811 0.553 0 0.355 0.648
8 8 Reef1 (9.398625 -19.95124) 0.309 0.0748 0.565 0 0.637 0.691
9 9 Reef1 (9.398625 -19.95124) 0.261 0.0599 0.599 0 0.961 0.575
10 10 Reef1 (9.398625 -19.95124) 0.246 0.0557 0.609 0 0.650 0.519
# ℹ 267,122 more rowsShow the code
benthos_reefs_temporal_summary <- benthos_reefs_sf |>
st_drop_geometry() |>
group_by(Year) |>
summarise(Mean = mean(HCC),
Median = median(HCC),
SD = sd(HCC),
Lower = quantile(HCC, 0.025),
Upper = quantile(HCC, 0.975))
benthos_reefs_temporal_summary |>
saveRDS(
file = paste0(
data_path,
"synthetic/benthos_reefs_temporal_summary.rds"
)
)Show the code
benthos_reefs_temporal_summary <-
readRDS(
file = paste0(
data_path,
"synthetic/benthos_reefs_temporal_summary.rds"
)
)
g <- benthos_reefs_temporal_summary |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = Lower, ymax = Upper), alpha = 0.2) +
geom_line(aes(x = Year, y = Mean, colour = "mean")) +
geom_line(aes(x = Year, y = Median, colour = "median")) +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_all_temporal_summary_plot.png"
),
g,
width = 8, height = 6, dpi = 72
)
These data represent a simulated sampling design in which 25 fixed (though randomly selected from among all the defined reefs in the synthetic spatial domain) reefs are monitored over a continuous 12 year period.
benthos_fixed_locs_obs <- readRDS(file = paste0(
data_path,
"synthetic/benthos_fixed_locs_obs.rds"
))
benthos_fixed_locs_obs <- benthos_fixed_locs_obs |>
left_join(
benthos_reefs_sf |>
st_drop_geometry() |>
dplyr::select(Year, Reef, CYC, DHW, OTHER) |>
group_by(Year, Reef) |>
summarise(across(c(CYC, DHW, OTHER), mean)),
by = c("Year", "Reef")
)`summarise()` has grouped output by 'Year'. You can override using the
`.groups` argument.# A tibble: 3,000 × 13
Reef Longitude Latitude Site Transect Year HCC SC MA
<chr> <dbl> <dbl> <chr> <chr> <int> <dbl> <dbl> <dbl>
1 Reef103 3.66 -19.5 S1 T1 1 18.6 3.99 54.9
2 Reef103 3.66 -19.5 S1 T1 2 19.0 4.32 50.5
3 Reef103 3.66 -19.5 S1 T1 3 19.0 4.18 48.7
4 Reef103 3.66 -19.5 S1 T1 4 21.6 4.81 50.2
5 Reef103 3.66 -19.5 S1 T1 5 25.4 5.20 49.0
6 Reef103 3.66 -19.5 S1 T1 6 21.2 5.47 46.3
7 Reef103 3.66 -19.5 S1 T1 7 22.2 5.26 45.6
8 Reef103 3.66 -19.5 S1 T1 8 23.2 5.18 46.1
9 Reef103 3.66 -19.5 S1 T1 9 19.7 4.43 47.2
10 Reef103 3.66 -19.5 S1 T1 10 17.4 4.09 54.4
# ℹ 2,990 more rows
# ℹ 4 more variables: Date <dttm>, CYC <dbl>, DHW <dbl>, OTHER <dbl>Show the code
benthos_fixed_locs_obs <- benthos_fixed_locs_obs |>
mutate(
fYear = as.factor(Year),
Reef = as.factor(Reef),
Site = interaction(Reef, Site),
Transect = interaction(Site, Transect),
cover = HCC/100
)
print(benthos_fixed_locs_obs)
## Raw summary means
all_sampled_sum <-
benthos_fixed_locs_obs |>
group_by(fYear) |>
summarise(
mean_mean_response = mean(cover),
median_median_response = median(cover),
mean_sd = sd(cover),
mean_lower_lower = mean_mean_response - 1.96 * mean_sd,
mean_upper_upper = mean_mean_response + 1.96 * mean_sd,
median_lower_lower = quantile(cover, 0.025),
median_upper_upper = quantile(cover, 0.975)
) |>
dplyr::select(-mean_sd) |>
pivot_longer(
cols = c(
mean_mean_response, median_median_response, mean_lower_lower,
mean_upper_upper, median_lower_lower, median_upper_upper
),
names_to = c("type", "variable", "stat"),
names_sep = "_",
values_to = "values"
) |>
mutate(type = paste("all_sampled_", type)) |>
dplyr::select(-variable) |>
pivot_wider(
names_from = stat,
values_from = values
) |>
mutate(Year = as.numeric(as.character(fYear)))
print(all_sampled_sum)Show the code
g <- all_sampled_sum |>
ggplot(aes(x = Year, y = response, color = type, fill = type)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
scale_y_continuous(
name = "Hard Coral Cover (%)",
limits = c(0, 1),
labels = scales::percent_format(accuracy = 1)
) +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_full_simple_raw_means_plot.png"
),
g,
width = 8, height = 6, dpi = 72
)
These data represent the fixed 25 sites, however some of the reefs have a temporal gap between years 2 and 6.
benthos_fixed_locs_obs_3 <- readRDS(
file = paste0(
data_path,
"synthetic/benthos_fixed_locs_obs_3.rds"
)
)
benthos_fixed_locs_obs_3 <- benthos_fixed_locs_obs_3 |>
left_join(
benthos_reefs_sf |>
st_drop_geometry() |>
dplyr::select(Year, Reef, CYC, DHW, OTHER) |>
group_by(Year, Reef) |>
summarise(across(c(CYC, DHW, OTHER), mean)),
by = c("Year", "Reef")
)`summarise()` has grouped output by 'Year'. You can override using the
`.groups` argument.# A tibble: 2,850 × 13
Reef Longitude Latitude Site Transect Year HCC SC MA
<chr> <dbl> <dbl> <chr> <chr> <int> <dbl> <dbl> <dbl>
1 Reef103 3.66 -19.5 S1 T1 1 18.6 3.99 54.9
2 Reef103 3.66 -19.5 S1 T1 2 19.0 4.32 50.5
3 Reef103 3.66 -19.5 S1 T1 3 19.0 4.18 48.7
4 Reef103 3.66 -19.5 S1 T1 4 21.6 4.81 50.2
5 Reef103 3.66 -19.5 S1 T1 5 25.4 5.20 49.0
6 Reef103 3.66 -19.5 S1 T1 6 21.2 5.47 46.3
7 Reef103 3.66 -19.5 S1 T1 7 22.2 5.26 45.6
8 Reef103 3.66 -19.5 S1 T1 8 23.2 5.18 46.1
9 Reef103 3.66 -19.5 S1 T1 9 19.7 4.43 47.2
10 Reef103 3.66 -19.5 S1 T1 10 17.4 4.09 54.4
# ℹ 2,840 more rows
# ℹ 4 more variables: Date <dttm>, CYC <dbl>, DHW <dbl>, OTHER <dbl>Show the code
g <- benthos_fixed_locs_obs_3 |>
dplyr::select(-MA, -SC) |>
pivot_wider(
id_cols = c(Reef, Longitude, Latitude, Site, Transect),
names_from = Year,
values_from = HCC
) |>
pivot_longer(
cols = -c(Reef, Longitude, Latitude, Site, Transect),
names_to = "Year",
values_to = "HCC"
) |>
mutate(Year = as.numeric(Year)) |>
ggplot() +
geom_line(aes(
y = HCC, x = Year, colour = Site,
group = interaction(Site, Transect)
)) +
facet_wrap(~Reef) +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_sampled_reefs_3_plot.png"
),
g,
width = 8, height = 6, dpi = 72
)
These data represent the fixed 25 sites, however all of the reefs have a temporal gap between years 2 and 6.
benthos_fixed_locs_obs_4 <- readRDS(
file = paste0(
data_path,
"synthetic/benthos_fixed_locs_obs_4.rds"
)
)
benthos_fixed_locs_obs_4 <- benthos_fixed_locs_obs_4 |>
left_join(
benthos_reefs_sf |>
st_drop_geometry() |>
dplyr::select(Year, Reef, CYC, DHW, OTHER) |>
group_by(Year, Reef) |>
summarise(across(c(CYC, DHW, OTHER), mean)),
by = c("Year", "Reef")
)`summarise()` has grouped output by 'Year'. You can override using the
`.groups` argument.# A tibble: 2,250 × 13
Reef Longitude Latitude Site Transect Year HCC SC MA
<chr> <dbl> <dbl> <chr> <chr> <int> <dbl> <dbl> <dbl>
1 Reef103 3.66 -19.5 S1 T1 1 18.6 3.99 54.9
2 Reef103 3.66 -19.5 S1 T1 2 19.0 4.32 50.5
3 Reef103 3.66 -19.5 S1 T1 6 21.2 5.47 46.3
4 Reef103 3.66 -19.5 S1 T1 7 22.2 5.26 45.6
5 Reef103 3.66 -19.5 S1 T1 8 23.2 5.18 46.1
6 Reef103 3.66 -19.5 S1 T1 9 19.7 4.43 47.2
7 Reef103 3.66 -19.5 S1 T1 10 17.4 4.09 54.4
8 Reef103 3.66 -19.5 S1 T1 11 18.5 4.17 53.2
9 Reef103 3.66 -19.5 S1 T1 12 17.1 4.13 49.1
10 Reef103 3.66 -19.5 S1 T2 1 20.6 3.88 48.2
# ℹ 2,240 more rows
# ℹ 4 more variables: Date <dttm>, CYC <dbl>, DHW <dbl>, OTHER <dbl>Show the code
yrs <- benthos_fixed_locs_obs_4 |>
pull(Year) |>
full_seq(period = 1)
g <- benthos_fixed_locs_obs_4 |>
complete(Year = yrs, nesting(Reef, Longitude, Latitude, Site, Transect)) |>
ggplot() +
geom_line(aes(
y = HCC, x = Year, colour = Site,
group = interaction(Site, Transect)
)) +
facet_wrap(~Reef) +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_sampled_reefs_4_plot.png"
),
g,
width = 8, height = 6, dpi = 72
)
4 Data preparations
Perform the following data preparation steps:
- specifically declare any categorical variables (as factors in R, categories in python)
- ensure that all levels of hierarchical (varying) effects have unique IDs
- express cover as a proportion out of 1 as is necessary for modelling against a beta distribution
- create a prediction grid containing all years in the monitoring range
5 Modelled trends
The following table lists datasets that are used in the comparison figures along with their respective legend entries.
view table
| Dataset | Description | Legend label(s) |
|---|---|---|
benthos_reefs_sf |
All reef-level data | |
benthos_reefs_temporal_summary |
All global-level data | all mean / all median |
all_sampled_sum |
Simple stats of all sampled data (25 sites) | all_sampled_mean / all_sampled_median |
benthos_fixed_locs_obs_0 |
All sampled data (25 sites) | |
mod_glmmTMB_0 |
glmmTMB model of all sampled data (25 sites) | |
glmmTMB_0_sum |
Estimated marginal means from all sampled data (25 sites) | glmmTMB |
mod_brms_0 |
brms model of all sampled data (25 sites) | |
brms_0_sum |
Estimated marginal means from all sampled data (25 sites) | brms |
mod_stan_0 |
stan model of all sampled data (25 sites) | |
stan_0_sum |
Estimated marginal means from all sampled data (25 sites) | stan |
mod_gbm_0 |
gbm model of all sampled data (25 sites) | |
gbm_0_sum |
Estimated predictions from all sampled data (25 sites) | gbm |
mod_gbm_0b |
gbm model with covariates of all sampled data (25 sites) | |
gbm_0b_sum |
Estimated predictions with covariates from all sampled data (25 sites) | gbm |
gbm_0c_sum |
Estimated predictions with covariates from all sampled data (all sites) | gbm |
mod_dbarts_0 |
bart model with covariates of all sampled ata (25 sites) | |
dbarts_0_sum |
Estimated predictions from all sampled data (25 sites) | dbarts |
mod_dbarts_0b |
bart model with covariates of all sampled ata (25 sites) | |
dbarts_0b_sum |
Estimated predictions with covariates from all sampled data (25 sites) | dbarts |
dbarts_0c_sum |
Estimated predictions with covariates from all data (all sites) | dbarts |
benthos_fixed_locs_obs_3 |
Gapped sampled data (25 sites) | |
mod_glmmTMB_3 |
glmmTMB model of gapped sampled data (25 sites) | |
glmmTMB_3_sum |
Estimated marginal means from gapped sampled data (25 sites) | glmmTMB |
mod_brms_3 |
brms model of gapped sampled data (25 sites) | |
brms_3_sum |
Estimated marginal means from gapped sampled data (25 sites) | brms |
mod_stan_3 |
stan model of gapped sampled data (25 sites) | |
stan_3_sum |
Estimated marginal means from gapped sampled data (25 sites) | stan |
mod_gbm_3 |
gbm model of gapped sampled data (25 sites) | |
gbm_3_sum |
Estimated predictions from gapped sampled data (25 sites) | gbm |
mod_gbm_3b |
gbm model with covariates of gapped sampled data (25 sites) | |
gbm_3b_sum |
Estimated predictions with covariates from gapped sampled data (25 sites) | gbm |
gbm_3c_sum |
Estimated predictions with covariates from gapped sampled data (gapped sites) | gbm |
mod_dbarts_3 |
bart model with covariates of gapped sampled ata (25 sites) | |
dbarts_3_sum |
Estimated predictions from gapped sampled data (25 sites) | dbarts |
mod_dbarts_3b |
bart model with covariates of gapped sampled ata (25 sites) | |
dbarts_3b_sum |
Estimated predictions with covariates from gapped sampled data (25 sites) | dbarts |
dbarts_3c_sum |
Estimated predictions with covariates from gapped data (all sites) | dbarts |
benthos_fixed_locs_obs_4 |
Gapped sampled data (all sites) | |
mod_glmmTMB_4 |
glmmTMB model of gapped sampled data (all sites) | |
glmmTMB_4_sum |
Estimated marginal means from gapped sampled data (all sites) | glmmTMB |
mod_brms_4 |
brms model of gapped sampled data (all sites) | |
brms_4_sum |
Estimated marginal means from gapped sampled data (all sites) | brms |
mod_stan_4 |
stan model of gapped sampled data (all sites) | |
stan_4_sum |
Estimated marginal means from gapped sampled data (all sites) | stan |
mod_gbm_4 |
gbm model of gapped sampled data (all sites) | |
gbm_4_sum |
Estimated predictions from gapped sampled data (all sites) | gbm |
mod_gbm_4b |
gbm model with covariates of gapped sampled data (all sites) | |
gbm_4b_sum |
Estimated predictions with covariates from gapped sampled data (all sites) | gbm |
gbm_4c_sum |
Estimated predictions with covariates from gapped sampled data (gapped sites) | gbm |
mod_dbarts_4 |
bart model with covariates of gapped sampled ata (all sites) | |
dbarts_4_sum |
Estimated predictions from gapped sampled data (all sites) | dbarts |
mod_dbarts_4b |
bart model with covariates of gapped sampled ata (all sites) | |
dbarts_4b_sum |
Estimated predictions with covariates from gapped sampled data (all sites) | dbarts |
dbarts_4c_sum |
Estimated predictions with covariates from gapped data (all sites) | dbarts |
| Model | Description |
|---|---|
| Simple | Simple hierarchical means/medians |
Show the code
glmmTMB_0_sum <-
mod_glmmTMB_0 |> emmeans(~fYear, at = newdata_0, type = "response") |>
as.data.frame() |>
rename(median = response, lower = asymp.LCL, upper = asymp.UCL) |>
mutate(Year = as.numeric(as.character(fYear))) |>
mutate(type = "glmmTMB")
saveRDS(glmmTMB_0_sum,
file = paste0(data_path, "synthetic/glmmTMB_0_sum.rds")
) Show the code
g <- glmmTMB_0_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "glmmTMB")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_glmmTMB_0.png"
),
g,
width = 8, height = 6, dpi = 72
)
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- glmmTMB modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
Show the code
glmmTMB_0_dharma <- readRDS(
file = paste0(data_path, "synthetic/glmmTMB_0_dharma.rds")
)
g <- wrap_elements(~testUniformity(glmmTMB_0_dharma)) +
wrap_elements(~plotResiduals(glmmTMB_0_dharma)) +
wrap_elements(~testDispersion(glmmTMB_0_dharma))
ggsave(
filename = paste0(
fig_path, "R_dharma_mod_glmmTMB_0.png"
),
g,
width = 10, height = 4, dpi = 72
)
Show the code
Family: beta ( logit )
Formula: cover ~ fYear + (1 | Site) + (1 | Transect)
Data: benthos_fixed_locs_obs_0
AIC BIC logLik -2*log(L) df.resid
-13130.3 -13040.2 6580.2 -13160.3 2985
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 0.08491 0.2914
Transect (Intercept) 0.02920 0.1709
Number of obs: 3000, groups: Site, 50; Transect, 250
Dispersion parameter for beta family (): 320
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.26554 0.04346 -29.122 < 2e-16 ***
fYear2 0.01799 0.01205 1.493 0.1353
fYear3 0.04983 0.01199 4.154 3.26e-05 ***
fYear4 0.10552 0.01192 8.853 < 2e-16 ***
fYear5 0.14555 0.01186 12.275 < 2e-16 ***
fYear6 0.18332 0.01180 15.539 < 2e-16 ***
fYear7 0.21699 0.01176 18.457 < 2e-16 ***
fYear8 0.18269 0.01180 15.486 < 2e-16 ***
fYear9 0.03129 0.01201 2.605 0.0092 **
fYear10 -0.06817 0.01218 -5.596 2.19e-08 ***
fYear11 -0.02343 0.01211 -1.935 0.0530 .
fYear12 -0.02145 0.01212 -1.770 0.0767 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Show the code
benthos_fixed_locs_obs_0 |>
group_by(fYear) |>
summarise(mean_abundance = qlogis(mean(cover)),
median_abundance = qlogis(median(cover)),
sd_abundance = sd(qlogis(cover)) ,
mad_abundance = mad(qlogis(cover))
)
priors <- c(
prior(normal(0, 0.5), class = "b"),
prior(normal(-1.5, 0.5), class = "Intercept"),
prior(student_t(3, 0, 0.5), class = "sd")
)
mod_form <- bf(cover ~ fYear + (1 | Site) + (1 | Transect),
family = "Beta"
)Show the code
brms_0_sum <-
mod_brms_0 |> emmeans(~fYear, at = newdata_0, type = "response") |>
as.data.frame() |>
rename(median = response, lower = lower.HPD, upper = upper.HPD) |>
mutate(Year = as.numeric(as.character(fYear))) |>
mutate(type = "brms")
saveRDS(brms_0_sum,
file = paste0(data_path, "synthetic/brms_0_sum.rds")
) Show the code
brms_0_sum <- readRDS(
file = paste0(data_path, "synthetic/brms_0_sum.rds")
)
g <-
brms_0_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "brms")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_brms_0.png"
),
g,
width = 8, height = 6, dpi = 72
)
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- brms modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
Show the code
Family: beta
Links: mu = logit; phi = identity
Formula: cover ~ fYear + (1 | Site) + (1 | Transect)
Data: benthos_fixed_locs_obs_0 (Number of observations: 3000)
Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
total post-warmup draws = 2400
Multilevel Hyperparameters:
~Site (Number of levels: 50)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.30 0.03 0.24 0.37 1.00 2084 2289
~Transect (Number of levels: 250)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.17 0.01 0.15 0.19 1.00 2193 2210
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.27 0.04 -1.35 -1.18 1.00 1643 2218
fYear2 0.02 0.01 -0.01 0.04 1.00 2004 2227
fYear3 0.05 0.01 0.03 0.07 1.00 2036 2327
fYear4 0.10 0.01 0.08 0.13 1.00 2206 2315
fYear5 0.15 0.01 0.12 0.17 1.00 2127 2197
fYear6 0.18 0.01 0.16 0.21 1.00 2213 2212
fYear7 0.22 0.01 0.19 0.24 1.00 2159 2181
fYear8 0.18 0.01 0.16 0.21 1.00 2177 2452
fYear9 0.03 0.01 0.01 0.06 1.00 2189 2136
fYear10 -0.07 0.01 -0.09 -0.04 1.00 2052 2191
fYear11 -0.02 0.01 -0.05 0.00 1.00 2272 2281
fYear12 -0.02 0.01 -0.05 0.00 1.00 2336 2369
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 317.91 8.51 301.48 334.14 1.00 2363 2265
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Show the code
benthos_fixed_locs_obs_0 <-
benthos_fixed_locs_obs_0 |>
mutate(
cYear = fYear,
grid_id = factor(Reef),
cSite = factor(interaction(Reef, Site)),
cReplicate = ifelse(is.na(Transect),
interaction(Site, Year),
interaction(cSite, Transect)),
Cover = cover,
area = 1,
sum = 1
)
saveRDS(benthos_fixed_locs_obs_0,
file = paste0(data_path, "synthetic/saveRDS(benthos_fixed_locs_obs_0_forstan.rds")
)
stan_data <- prepare_data_for_stan(benthos_fixed_locs_obs_0, yrs = NULL)
model_stan <- cmdstanr::cmdstan_model(stan_file = "model1.stan")
## model_stan <- cmdstanr::cmdstan_model(stan_file = "mod1a.stan")
## model_stan <- cmdstanr::cmdstan_model(stan_file = "mod2a.stan")Show the code
Show the code
stan_0_sum <-
mod_stan_0$draws(variables = "cellmeans") |>
posterior::as_draws_df() |>
posterior::summarise_draws(
median,
HDInterval::hdi,
~ HDInterval::hdi(., credMass = c(0.9)),
rhat,
ess_bulk,
ess_tail
) |>
rename(lower_90 = V4, upper_90 = V5) |>
bind_cols(newdata_0) |>
mutate(Year = as.numeric(as.character(fYear)))
saveRDS(stan_0_sum,
file = paste0(data_path, "synthetic/stan_0_sum.rds")
) Show the code
stan_0_sum <- readRDS(
file = paste0(data_path, "synthetic/stan_0_sum.rds")
)
g <- stan_0_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "stan")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_stan_0.png"
),
g,
width = 8, height = 6, dpi = 72
)
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- stan modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
Show the code
mod_stan_0 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_0.rds")
)
color_scheme_set("viridis")
g <-
mod_stan_0$draws(variables = c("beta", "phi", "sd_1", "sd_2", "sd_3")) |>
mcmc_trace() +
theme_minimal()
ggsave(
filename = paste0(
fig_path, "R_stan_trace_0.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
mod_stan_0 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_0.rds")
)
color_scheme_set("viridis")
g <-
mod_stan_0$draws(variables = c("beta", "phi", "sd_1", "sd_2", "sd_3")) |>
mcmc_acf() +
theme_minimal()
ggsave(
filename = paste0(
fig_path, "R_stan_ac_0.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
mod_stan_0 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_0.rds")
)
g <-
bayesplot::pp_check(
benthos_fixed_locs_obs_0$cover,
mod_stan_0$draws("ypred", format = "matrix")[1:100, ],
ppc_dens_overlay
) +
theme_classic()
ggsave(
filename = paste0(
fig_path, "R_stan_ppc_0.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
phi 317.92 317.93 8.69 8.74 303.53 331.90 1.00 2854 2685
beta[1] -1.26 -1.26 0.05 0.04 -1.33 -1.18 1.00 2148 2386
beta[2] -1.24 -1.24 0.05 0.05 -1.32 -1.17 1.00 2210 2406
beta[3] -1.21 -1.21 0.05 0.05 -1.28 -1.13 1.00 2205 2393
beta[4] -1.15 -1.15 0.05 0.05 -1.23 -1.08 1.00 2164 2251
beta[5] -1.11 -1.11 0.05 0.05 -1.19 -1.04 1.00 2186 2374
beta[6] -1.08 -1.08 0.05 0.05 -1.15 -1.00 1.00 2149 2254
beta[7] -1.04 -1.04 0.05 0.05 -1.12 -0.97 1.00 2164 2571
beta[8] -1.08 -1.08 0.05 0.05 -1.15 -1.00 1.00 2168 2299
beta[9] -1.23 -1.23 0.05 0.05 -1.30 -1.15 1.00 2170 2365
beta[10] -1.33 -1.32 0.05 0.05 -1.40 -1.25 1.00 2174 2289
beta[11] -1.28 -1.28 0.05 0.05 -1.36 -1.21 1.00 2168 2391
beta[12] -1.28 -1.28 0.05 0.05 -1.35 -1.20 1.00 2206 2487
sd_1[1] 0.30 0.29 0.05 0.05 0.22 0.39 1.00 2892 2909
sd_2[1] 0.09 0.09 0.03 0.02 0.05 0.14 1.00 1525 1450
sd_3[1] 0.17 0.17 0.01 0.01 0.16 0.19 1.00 2469 2860 variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
cellmeans[1] 0.22 0.22 0.01 0.01 0.21 0.23 1.00 2148 2386
cellmeans[2] 0.22 0.22 0.01 0.01 0.21 0.24 1.00 2210 2393
cellmeans[3] 0.23 0.23 0.01 0.01 0.22 0.24 1.00 2205 2393
cellmeans[4] 0.24 0.24 0.01 0.01 0.23 0.25 1.00 2164 2251
cellmeans[5] 0.25 0.25 0.01 0.01 0.23 0.26 1.00 2186 2374
cellmeans[6] 0.25 0.25 0.01 0.01 0.24 0.27 1.00 2149 2254
cellmeans[7] 0.26 0.26 0.01 0.01 0.25 0.28 1.00 2164 2571
cellmeans[8] 0.25 0.25 0.01 0.01 0.24 0.27 1.00 2168 2299
cellmeans[9] 0.23 0.23 0.01 0.01 0.21 0.24 1.00 2170 2365
cellmeans[10] 0.21 0.21 0.01 0.01 0.20 0.22 1.00 2174 2289
cellmeans[11] 0.22 0.22 0.01 0.01 0.20 0.23 1.00 2168 2391
cellmeans[12] 0.22 0.22 0.01 0.01 0.21 0.23 1.00 2207 2487Show the code
gbm_0_sum <- readRDS(
file = paste0(data_path, "synthetic/gbm_0_sum.rds")
)
g <-
gbm_0_sum |>
ggplot() +
geom_line(aes(x = Year, y = median, color = "gbm")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_gbm_0.png"
),
g,
width = 8, height = 6, dpi = 72
)
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- gbm modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
- gbm derived predictions lack uncertainty estimates
Show the code
mod_gbm_0b <- gbm(cover ~ fYear + Latitude + Longitude + CYC + DHW + OTHER,
data = benthos_fixed_locs_obs_0,
distribution = "gaussian",
var.monotone = c(0, 0, 0, -1, -1, -1),
n.trees = 10000,
interaction.depth = 5,
shrinkage = 0.001,
bag.fraction = 0.5,
cv.folds = 5,
verbose = TRUE
)
saveRDS(mod_gbm_0b,
file = paste0(data_path, "synthetic/mod_gbm_0b.rds")
)Show the code
newdata_0b <- benthos_fixed_locs_obs_0
gbm_0b_sum <- newdata_0b |>
mutate(median = predict(mod_gbm_0b, newdata_0b, n.trees = n.trees, type = "response")) |>
mutate(Year = as.numeric(as.character(fYear))) |>
## group_by(Year, Site, Transect) |>
## summarise(median = median(median)) |>
## group_by(Year, Site, .add = FALSE) |>
## summarise(median = median(median)) |>
group_by(Year, .add = FALSE) |>
summarise(median = median(median)) |>
mutate(type = "gbm")
saveRDS(gbm_0b_sum,
file = paste0(data_path, "synthetic/gbm_0b_sum.rds")
) 
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- gbm modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
- gbm derived predictions lack uncertainty estimates
- relative influences mainly space/time
Show the code
dbarts_0c_sum <- preds_0c |>
t() |>
as_tibble(.name_repair = ~ paste0("V", seq_along(.))) |>
bind_cols(benthos_reefs_sf) |>
pivot_longer(
cols = matches("^V[0-9]*"),
names_to = ".draw", values_to = "fit"
) |>
group_by(Year, .draw) |>
summarise(fit = mean(fit)) |>
## dplyr::select(Year, fit, .draw) |>
group_by(Year, .add = FALSE) |>
summarise_draws(median, HDInterval::hdi)
saveRDS(dbarts_0c_sum,
file = paste0(data_path, "synthetic/dbarts_0c_sum.rds")
) 
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- gbm modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
- gbm derived predictions lack uncertainty estimates
- relative influences mainly space/time
Show the code
dbarts_0b_sum <- preds |>
t() |>
as_tibble(.name_repair = ~ paste0("V", seq_along(.))) |>
bind_cols(benthos_fixed_locs_obs_0) |>
pivot_longer(
cols = matches("^V[0-9]*"),
names_to = ".draw", values_to = "fit"
) |>
group_by(Year, .draw) |>
summarise(fit = mean(fit)) |>
## dplyr::select(Year, fit, .draw) |>
group_by(Year, .add = FALSE) |>
summarise_draws(median, HDInterval::hdi)
saveRDS(dbarts_0b_sum,
file = paste0(data_path, "synthetic/dbarts_0b_sum.rds")
) 
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- gbm modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
- gbm derived predictions lack uncertainty estimates
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- dbarts modelled estimates are consistent with the 25 site same data, yet biased lower than the “truth”
- sampled reefs (blue and green dashed lines) are lower than the true global average (red and khaki)
- dbarts modelled estimates are slightly elevated above the 25 site same data, yet still biased lower than the “truth”
Show the code
glmmTMB_3_sum <-
mod_glmmTMB_3 |> emmeans(~fYear, at = newdata_3, type = "response") |>
as.data.frame() |>
rename(median = response, lower = asymp.LCL, upper = asymp.UCL) |>
mutate(Year = as.numeric(as.character(fYear))) |>
mutate(type = "glmmTMB")
saveRDS(glmmTMB_3_sum,
file = paste0(data_path, "synthetic/glmmTMB_3_sum.rds")
) Show the code
g1 <- glmmTMB_3_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "glmmTMB")) +
geom_line(data = mod_simple_3,
aes(x = Year, y = Mean, colour = "simple data mean"), linetype = "dashed") +
geom_line(data = mod_simple_3,
aes(x = Year, y = Median, colour = "simple data median"), linetype = "dashed") +
## geom_line(data = all_sampled_sum,
## aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
g2 <- glmmTMB_3_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "glmmTMB")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "true mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "true median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_glmmTMB_3.png"
),
g1 + g2,
width = 12, height = 6, dpi = 72
)
- left figure compares glmmTMB modelled estimates to the sample data
- right figure compares glmmTMB modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and khaki dashed lines) are lower than the true global average (blue and purple) - right side figure
- glmmTMB modelled estimates are consistent with the 25 site, replacement data, yet biased lower than the “truth”
Show the code
glmmTMB_3_dharma <- readRDS(
file = paste0(data_path, "synthetic/glmmTMB_3_dharma.rds")
)
g <- wrap_elements(~testUniformity(glmmTMB_3_dharma)) +
wrap_elements(~plotResiduals(glmmTMB_3_dharma)) +
wrap_elements(~testDispersion(glmmTMB_3_dharma))
ggsave(
filename = paste0(
fig_path, "R_dharma_mod_glmmTMB_3.png"
),
g,
width = 10, height = 4, dpi = 72
)
Show the code
Family: beta ( logit )
Formula: cover ~ fYear + (1 | Site) + (1 | Transect)
Data: benthos_fixed_locs_obs_3
AIC BIC logLik -2*log(L) df.resid
-12463.9 -12374.6 6247.0 -12493.9 2835
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 0.08730 0.2955
Transect (Intercept) 0.02897 0.1702
Number of obs: 2850, groups: Site, 50; Transect, 250
Dispersion parameter for beta family (): 322
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.26650 0.04399 -28.791 < 2e-16 ***
fYear2 0.01813 0.01202 1.507 0.13170
fYear3 0.05521 0.01279 4.316 1.59e-05 ***
fYear4 0.11058 0.01270 8.706 < 2e-16 ***
fYear5 0.14189 0.01265 11.220 < 2e-16 ***
fYear6 0.18377 0.01177 15.608 < 2e-16 ***
fYear7 0.21746 0.01173 18.532 < 2e-16 ***
fYear8 0.18328 0.01177 15.565 < 2e-16 ***
fYear9 0.03173 0.01199 2.646 0.00815 **
fYear10 -0.06779 0.01216 -5.576 2.46e-08 ***
fYear11 -0.02289 0.01208 -1.894 0.05822 .
fYear12 -0.02085 0.01209 -1.724 0.08467 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Show the code
benthos_fixed_locs_obs_3 |>
group_by(fYear) |>
summarise(mean_abundance = qlogis(mean(cover)),
median_abundance = qlogis(median(cover)),
sd_abundance = sd(qlogis(cover)) ,
mad_abundance = mad(qlogis(cover))
)
priors <- c(
prior(normal(0, 0.5), class = "b"),
prior(normal(-1.5, 0.5), class = "Intercept"),
prior(student_t(3, 0, 0.5), class = "sd")
)
mod_form <- bf(cover ~ fYear + (1 | Site) + (1 | Transect),
family = "Beta"
)Show the code
brms_3_sum <-
mod_brms_3 |> emmeans(~fYear, at = newdata_3, type = "response") |>
as.data.frame() |>
rename(median = response, lower = lower.HPD, upper = upper.HPD) |>
mutate(Year = as.numeric(as.character(fYear))) |>
mutate(type = "brms")
saveRDS(brms_3_sum,
file = paste0(data_path, "synthetic/brms_3_sum.rds")
) Show the code
brms_3_sum <- readRDS(
file = paste0(data_path, "synthetic/brms_3_sum.rds")
)
g1 <-
brms_3_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "brms")) +
geom_line(data = mod_simple_3,
aes(x = Year, y = Mean, colour = "simple data mean"), linetype = "dashed") +
geom_line(data = mod_simple_3,
aes(x = Year, y = Median, colour = "simple data median"), linetype = "dashed") +
theme_bw()
g2 <-
brms_3_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "brms")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_brms_3.png"
),
g1 + g2,
width = 12, height = 6, dpi = 72
)
- left figure compares brms modelled estimates to the sample data
- right figure compares brms modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- brms modelled estimates are consistent with the 25 site, replacement data, yet biased lower than the “truth”
Show the code
Family: beta
Links: mu = logit; phi = identity
Formula: cover ~ fYear + (1 | Site) + (1 | Transect)
Data: benthos_fixed_locs_obs_3 (Number of observations: 2850)
Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
total post-warmup draws = 2400
Multilevel Hyperparameters:
~Site (Number of levels: 50)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.30 0.03 0.24 0.38 1.00 1469 1851
~Transect (Number of levels: 250)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.17 0.01 0.15 0.19 1.00 2025 2060
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.27 0.05 -1.36 -1.18 1.00 1950 2142
fYear2 0.02 0.01 -0.01 0.04 1.00 2408 2367
fYear3 0.05 0.01 0.03 0.08 1.00 2219 2061
fYear4 0.11 0.01 0.09 0.14 1.00 2286 2147
fYear5 0.14 0.01 0.12 0.17 1.00 2071 1986
fYear6 0.18 0.01 0.16 0.21 1.00 2198 2117
fYear7 0.22 0.01 0.19 0.24 1.00 2249 2179
fYear8 0.18 0.01 0.16 0.21 1.00 2348 2389
fYear9 0.03 0.01 0.01 0.06 1.00 2320 2329
fYear10 -0.07 0.01 -0.09 -0.04 1.00 2361 2329
fYear11 -0.02 0.01 -0.05 -0.00 1.00 2269 2365
fYear12 -0.02 0.01 -0.05 0.00 1.00 2310 2406
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 319.22 8.86 302.48 336.81 1.00 2333 1974
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Show the code
benthos_fixed_locs_obs_3 <-
benthos_fixed_locs_obs_3 |>
mutate(
cYear = fYear,
grid_id = factor(Reef),
cSite = factor(interaction(Reef, Site)),
cReplicate = ifelse(is.na(Transect),
interaction(Site, Year),
interaction(cSite, Transect)),
Cover = cover,
area = 1,
sum = 1
)
saveRDS(benthos_fixed_locs_obs_3,
file = paste0(data_path, "synthetic/saveRDS(benthos_fixed_locs_obs_3_forstan.rds")
)
stan_data <- prepare_data_for_stan(benthos_fixed_locs_obs_3, yrs = NULL)
model_stan <- cmdstanr::cmdstan_model(stan_file = "model1.stan")Show the code
Show the code
stan_3_sum <-
mod_stan_3$draws(variables = "cellmeans") |>
posterior::as_draws_df() |>
posterior::summarise_draws(
median,
HDInterval::hdi,
~ HDInterval::hdi(., credMass = c(0.9)),
rhat,
ess_bulk,
ess_tail
) |>
rename(lower_90 = V4, upper_90 = V5) |>
bind_cols(newdata_3) |>
mutate(Year = as.numeric(as.character(fYear)))
saveRDS(stan_3_sum,
file = paste0(data_path, "synthetic/stan_3_sum.rds")
) Show the code
stan_3_sum <- readRDS(
file = paste0(data_path, "synthetic/stan_3_sum.rds")
)
g1 <- stan_3_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "stan")) +
geom_line(data = mod_simple_3,
aes(x = Year, y = Mean, colour = "simple data mean"), linetype = "dashed") +
geom_line(data = mod_simple_3,
aes(x = Year, y = Median, colour = "simple data median"), linetype = "dashed") +
theme_bw()
g2 <- stan_3_sum |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "stan")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_stan_3.png"
),
g1 + g2,
width = 12, height = 6, dpi = 72
)
- left figure compares stan modelled estimates to the sample data
- right figure compares stan modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- stan modelled estimates are consistent with the 25 site, replacement data, yet biased lower than the “truth”
Show the code
mod_stan_3 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_3.rds")
)
color_scheme_set("viridis")
g <-
mod_stan_3$draws(variables = c("beta", "phi", "sd_1", "sd_2", "sd_3")) |>
mcmc_trace() +
theme_minimal()
ggsave(
filename = paste0(
fig_path, "R_stan_trace_3.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
mod_stan_3 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_3.rds")
)
color_scheme_set("viridis")
g <-
mod_stan_3$draws(variables = c("beta", "phi", "sd_1", "sd_2", "sd_3")) |>
mcmc_acf() +
theme_minimal()
ggsave(
filename = paste0(
fig_path, "R_stan_ac_3.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
mod_stan_3 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_3.rds")
)
g <-
bayesplot::pp_check(
benthos_fixed_locs_obs_3$cover,
mod_stan_3$draws("ypred", format = "matrix")[1:100, ],
ppc_dens_overlay
) +
theme_classic()
ggsave(
filename = paste0(
fig_path, "R_stan_ppc_3.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
phi 319.44 319.49 8.79 8.80 305.40 333.97 1.00 3042 2868
beta[1] -1.26 -1.26 0.05 0.05 -1.33 -1.18 1.00 2433 2866
beta[2] -1.24 -1.24 0.05 0.05 -1.32 -1.17 1.00 2408 2720
beta[3] -1.20 -1.20 0.05 0.05 -1.28 -1.13 1.00 2441 2755
beta[4] -1.15 -1.15 0.05 0.05 -1.23 -1.07 1.00 2463 2673
beta[5] -1.12 -1.12 0.05 0.05 -1.19 -1.04 1.00 2420 2732
beta[6] -1.08 -1.08 0.05 0.05 -1.15 -1.00 1.00 2457 2604
beta[7] -1.04 -1.04 0.05 0.05 -1.12 -0.97 1.00 2415 2793
beta[8] -1.08 -1.08 0.05 0.05 -1.15 -1.00 1.00 2458 2902
beta[9] -1.23 -1.23 0.05 0.05 -1.30 -1.15 1.00 2464 2687
beta[10] -1.33 -1.33 0.05 0.05 -1.40 -1.25 1.00 2417 2829
beta[11] -1.28 -1.28 0.05 0.05 -1.36 -1.21 1.00 2475 2738
beta[12] -1.28 -1.28 0.05 0.05 -1.36 -1.20 1.00 2499 2827
sd_1[1] 0.30 0.29 0.05 0.05 0.23 0.39 1.00 2850 2828
sd_2[1] 0.09 0.09 0.03 0.02 0.05 0.13 1.00 1866 2047
sd_3[1] 0.17 0.17 0.01 0.01 0.16 0.19 1.00 2420 2743 variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
cellmeans[1] 0.22 0.22 0.01 0.01 0.21 0.23 1.00 2433 2866
cellmeans[2] 0.22 0.22 0.01 0.01 0.21 0.24 1.00 2409 2720
cellmeans[3] 0.23 0.23 0.01 0.01 0.22 0.24 1.00 2441 2755
cellmeans[4] 0.24 0.24 0.01 0.01 0.23 0.25 1.00 2463 2673
cellmeans[5] 0.25 0.25 0.01 0.01 0.23 0.26 1.00 2420 2732
cellmeans[6] 0.25 0.25 0.01 0.01 0.24 0.27 1.00 2457 2604
cellmeans[7] 0.26 0.26 0.01 0.01 0.25 0.28 1.00 2415 2793
cellmeans[8] 0.25 0.25 0.01 0.01 0.24 0.27 1.00 2458 2902
cellmeans[9] 0.23 0.23 0.01 0.01 0.21 0.24 1.00 2464 2687
cellmeans[10] 0.21 0.21 0.01 0.01 0.20 0.22 1.00 2417 2829
cellmeans[11] 0.22 0.22 0.01 0.01 0.20 0.23 1.00 2475 2738
cellmeans[12] 0.22 0.22 0.01 0.01 0.20 0.23 1.00 2499 2827Show the code
gbm_3_sum <- readRDS(
file = paste0(data_path, "synthetic/gbm_3_sum.rds")
)
g1 <-
gbm_3_sum |>
ggplot() +
geom_line(aes(x = Year, y = median, color = "gbm")) +
geom_line(data = mod_simple_3,
aes(x = Year, y = Mean, colour = "simple data mean"), linetype = "dashed") +
geom_line(data = mod_simple_3,
aes(x = Year, y = Median, colour = "simple data median"), linetype = "dashed") +
theme_bw()
g2 <-
gbm_3_sum |>
ggplot() +
geom_line(aes(x = Year, y = median, color = "gbm")) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_gbm_3.png"
),
g1 + g2,
width = 12, height = 6, dpi = 72
)
- left figure compares gbm modelled estimates to the sample data
- right figure compares gbm modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- gbm modelled estimates are consistent with the 25 site, replacement data, yet biased lower than the “truth”
- gbm derived predictions lack uncertainty estimates
Show the code
mod_gbm_3b <- gbm(cover ~ fYear + Latitude + Longitude + CYC + DHW + OTHER,
data = benthos_fixed_locs_obs_3,
distribution = "gaussian",
var.monotone = c(0, 0, 0, -1, -1, -1),
n.trees = 10000,
interaction.depth = 5,
shrinkage = 0.001,
bag.fraction = 0.5,
cv.folds = 5,
verbose = TRUE
)
saveRDS(mod_gbm_3b,
file = paste0(data_path, "synthetic/mod_gbm_3b.rds")
)Show the code
newdata_3b <- benthos_fixed_locs_obs_3
gbm_3b_sum <- newdata_3b |>
mutate(median = predict(mod_gbm_3b, newdata_3b, n.trees = n.trees, type = "response")) |>
mutate(Year = as.numeric(as.character(fYear))) |>
## group_by(Year, Site, Transect) |>
## summarise(median = median(median)) |>
## group_by(Year, Site, .add = FALSE) |>
## summarise(median = median(median)) |>
group_by(Year, .add = FALSE) |>
summarise(median = median(median)) |>
mutate(type = "gbm")
saveRDS(gbm_3b_sum,
file = paste0(data_path, "synthetic/gbm_3b_sum.rds")
) 
- left figure compares gbm modelled estimates to the sample data
- right figure compares gbm modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- gbm modelled estimates are consistent with the 25 site, replacement data, yet biased lower than the “truth”
- gbm derived predictions lack uncertainty estimates
- relative influences mainly space/time
Show the code
dbarts_3c_sum <- preds_3c |>
t() |>
as_tibble(.name_repair = ~ paste0("V", seq_along(.))) |>
bind_cols(benthos_reefs_sf) |>
pivot_longer(
cols = matches("^V[0-9]*"),
names_to = ".draw", values_to = "fit"
) |>
group_by(Year, .draw) |>
summarise(fit = mean(fit)) |>
## dplyr::select(Year, fit, .draw) |>
group_by(Year, .add = FALSE) |>
summarise_draws(median, HDInterval::hdi)
saveRDS(dbarts_3c_sum,
file = paste0(data_path, "synthetic/dbarts_3c_sum.rds")
) 
- left figure compares gbm modelled estimates to the sample data
- right figure compares gbm modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- gbm modelled estimates are consistent with the 25 site, replacement data, yet biased lower than the “truth”
- gbm derived predictions lack uncertainty estimates
- relative influences mainly space/time (though note that this is because this partial is the result of the same model, just predicted onto the full data grid)
Show the code
dbarts_3b_sum <- preds |>
t() |>
as_tibble(.name_repair = ~ paste0("V", seq_along(.))) |>
bind_cols(benthos_fixed_locs_obs_3) |>
pivot_longer(
cols = matches("^V[0-9]*"),
names_to = ".draw", values_to = "fit"
) |>
group_by(Year, .draw) |>
summarise(fit = mean(fit)) |>
## dplyr::select(Year, fit, .draw) |>
group_by(Year, .add = FALSE) |>
summarise_draws(median, HDInterval::hdi)
saveRDS(dbarts_3b_sum,
file = paste0(data_path, "synthetic/dbarts_3b_sum.rds")
) 
- left figure compares dbarts modelled estimates to the sample data
- right figure compares dbarts modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- dbarts modelled estimates are slightly lower than the 25 site, replacement data, and thus biased lower than the “truth”
- left figure compares dbarts modelled estimates to the sample data
- right figure compares dbarts modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- dbarts modelled estimates are consistent with the 25 site, replacement data, yet biased lower than the “truth”
- relative influences mainly space/time
- left figure compares dbarts modelled estimates to the sample data
- right figure compares dbarts modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- dbarts modelled estimates are slightly higher than the 25 site, replacement data, yet biased lower than the “truth”
- relative influences mainly space/time (though note that this is because this partial is the result of the same model, just predicted onto the full data grid)
Show the code
glmmTMB_4_sum <-
mod_glmmTMB_4 |> emmeans(~fYear, at = newdata_4, type = "response") |>
as.data.frame() |>
rename(median = response, lower = asymp.LCL, upper = asymp.UCL) |>
mutate(Year = as.numeric(as.character(fYear))) |>
mutate(type = "glmmTMB")
saveRDS(glmmTMB_4_sum,
file = paste0(data_path, "synthetic/glmmTMB_4_sum.rds")
) Show the code
all_years <- glmmTMB_4_sum |>
reframe(Year = full_seq(Year, period = 1))
glmmTMB_4_sum_sam_yrs <- glmmTMB_4_sum |>
full_join(all_years, by = "Year")
gap_years_boundary <- glmmTMB_4_sum |>
reframe(Year = range(setdiff(full_seq(Year, period = 1), Year)) + c(-1, 1))
glmmTMB_4_sum_gap_yrs <- glmmTMB_4_sum |>
right_join(gap_years_boundary, by = "Year")
g1 <- glmmTMB_4_sum_sam_yrs |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "glmmTMB")) +
geom_ribbon(data = glmmTMB_4_sum_gap_yrs,
aes(x = Year, ymin = lower, ymax = upper), alpha = 0.1) +
geom_line(data = glmmTMB_4_sum_gap_yrs,
aes(x = Year, y = median, color = "glmmTMB"), alpha = 0.4) +
geom_line(data = mod_simple_4,
aes(x = Year, y = Mean, colour = "simple data mean"), linetype = "dashed") +
geom_line(data = mod_simple_4,
aes(x = Year, y = Median, colour = "simple data median"), linetype = "dashed") +
## geom_line(data = all_sampled_sum,
## aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
g2 <- glmmTMB_4_sum_sam_yrs |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "glmmTMB")) +
geom_ribbon(data = glmmTMB_4_sum_gap_yrs,
aes(x = Year, ymin = lower, ymax = upper), alpha = 0.1) +
geom_line(data = glmmTMB_4_sum_gap_yrs,
aes(x = Year, y = median, color = "glmmTMB"), alpha = 0.4) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "true mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "true median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_glmmTMB_4.png"
),
g1 + g2,
width = 12, height = 6, dpi = 72
)
- left figure compares glmmTMB modelled estimates to the sample data
- right figure compares glmmTMB modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- glmmTMB modelled estimates are slightly higher than the 5 site, replacement data, yet biased lower than the “truth”
- the uncertainty is relatively high
Show the code
glmmTMB_4_dharma <- readRDS(
file = paste0(data_path, "synthetic/glmmTMB_4_dharma.rds")
)
g <- wrap_elements(~testUniformity(glmmTMB_4_dharma)) +
wrap_elements(~plotResiduals(glmmTMB_4_dharma)) +
wrap_elements(~testDispersion(glmmTMB_4_dharma))
ggsave(
filename = paste0(
fig_path, "R_dharma_mod_glmmTMB_4.png"
),
g,
width = 10, height = 4, dpi = 72
)
Show the code
Family: beta ( logit )
Formula: cover ~ fYear + (1 | Site) + (1 | Transect)
Data: benthos_fixed_locs_obs_4
AIC BIC logLik -2*log(L) df.resid
-9675.2 -9606.5 4849.6 -9699.2 2238
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Site (Intercept) 0.08850 0.2975
Transect (Intercept) 0.02837 0.1684
Number of obs: 2250, groups: Site, 50; Transect, 250
Dispersion parameter for beta family (): 314
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.26755 0.04425 -28.643 < 2e-16 ***
fYear2 0.01827 0.01217 1.502 0.13316
fYear6 0.18491 0.01191 15.521 < 2e-16 ***
fYear7 0.21852 0.01187 18.405 < 2e-16 ***
fYear8 0.18444 0.01191 15.482 < 2e-16 ***
fYear9 0.03309 0.01213 2.728 0.00638 **
fYear10 -0.06658 0.01230 -5.413 6.19e-08 ***
fYear11 -0.02140 0.01223 -1.750 0.08008 .
fYear12 -0.01944 0.01223 -1.589 0.11214
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Show the code
benthos_fixed_locs_obs_4 |>
group_by(fYear) |>
summarise(mean_abundance = qlogis(mean(cover)),
median_abundance = qlogis(median(cover)),
sd_abundance = sd(qlogis(cover)) ,
mad_abundance = mad(qlogis(cover))
)
priors <- c(
prior(normal(0, 0.5), class = "b"),
prior(normal(-1.5, 0.5), class = "Intercept"),
prior(student_t(3, 0, 0.5), class = "sd")
)
mod_form <- bf(cover ~ fYear + (1 | Site) + (1 | Transect),
family = "Beta"
)Show the code
brms_4_sum <-
mod_brms_4 |> emmeans(~fYear, at = newdata_4, type = "response") |>
as.data.frame() |>
rename(median = response, lower = lower.HPD, upper = upper.HPD) |>
mutate(Year = as.numeric(as.character(fYear))) |>
mutate(type = "brms")
saveRDS(brms_4_sum,
file = paste0(data_path, "synthetic/brms_4_sum.rds")
) Show the code
brms_4_sum <- readRDS(
file = paste0(data_path, "synthetic/brms_4_sum.rds")
)
all_years <- brms_4_sum |>
reframe(Year = full_seq(Year, period = 1))
brms_4_sum_sam_yrs <- brms_4_sum |>
full_join(all_years, by = "Year")
gap_years_boundary <- brms_4_sum |>
reframe(Year = range(setdiff(full_seq(Year, period = 1), Year)) + c(-1, 1))
brms_4_sum_gap_yrs <- brms_4_sum |>
right_join(gap_years_boundary, by = "Year")
g1 <-
brms_4_sum_sam_yrs |>
full_join(all_years, by = "Year") |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "brms")) +
geom_ribbon(data = brms_4_sum_gap_yrs,
aes(x = Year, ymin = lower, ymax = upper), alpha = 0.1) +
geom_line(data = brms_4_sum_gap_yrs,
aes(x = Year, y = median, color = "brms"), alpha = 0.4) +
geom_line(data = mod_simple_4,
aes(x = Year, y = Mean, colour = "simple data mean"), linetype = "dashed") +
geom_line(data = mod_simple_4,
aes(x = Year, y = Median, colour = "simple data median"), linetype = "dashed") +
theme_bw()
g2 <-
brms_4_sum_sam_yrs |>
full_join(all_years, by = "Year") |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "brms")) +
geom_ribbon(data = brms_4_sum_gap_yrs,
aes(x = Year, ymin = lower, ymax = upper), alpha = 0.1) +
geom_line(data = brms_4_sum_gap_yrs,
aes(x = Year, y = median, color = "brms"), alpha = 0.4) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_brms_4.png"
),
g1 + g2,
width = 12, height = 6, dpi = 72
)
- left figure compares brms modelled estimates to the sample data
- right figure compares brms modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- brms modelled estimates are slightly higher than the 5 site, replacement data, yet biased lower than the “truth”
- the uncertainty is relatively high
Show the code
Family: beta
Links: mu = logit; phi = identity
Formula: cover ~ fYear + (1 | Site) + (1 | Transect)
Data: benthos_fixed_locs_obs_4 (Number of observations: 2250)
Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
total post-warmup draws = 2400
Multilevel Hyperparameters:
~Site (Number of levels: 50)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.31 0.03 0.25 0.38 1.00 1858 2204
~Transect (Number of levels: 250)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.17 0.01 0.15 0.19 1.00 1940 2275
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.27 0.05 -1.35 -1.18 1.00 1973 2216
fYear2 0.02 0.01 -0.01 0.04 1.00 2155 2279
fYear6 0.19 0.01 0.16 0.21 1.00 2303 2287
fYear7 0.22 0.01 0.19 0.24 1.00 2056 2291
fYear8 0.18 0.01 0.16 0.21 1.00 2105 2177
fYear9 0.03 0.01 0.01 0.06 1.00 2244 2237
fYear10 -0.07 0.01 -0.09 -0.04 1.00 2054 2276
fYear11 -0.02 0.01 -0.04 0.00 1.00 2269 2290
fYear12 -0.02 0.01 -0.04 0.00 1.00 2333 2456
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 312.27 10.09 293.15 332.33 1.00 2340 2269
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Show the code
benthos_fixed_locs_obs_4 <-
benthos_fixed_locs_obs_4 |>
mutate(fYear = factor(fYear, levels = 1:12)) |>
mutate(
cYear = fYear,
grid_id = factor(Reef),
cSite = factor(interaction(Reef, Site)),
cReplicate = ifelse(is.na(Transect),
interaction(Site, Year),
interaction(cSite, Transect)),
Cover = cover,
area = 1,
sum = 1
)
saveRDS(benthos_fixed_locs_obs_4,
file = paste0(data_path, "synthetic/saveRDS(benthos_fixed_locs_obs_4_forstan.rds")
)
stan_data <- prepare_data_for_stan(benthos_fixed_locs_obs_4, yrs = 1:12)
model_stan <- cmdstanr::cmdstan_model(stan_file = "model2.stan")Show the code
Show the code
stan_4_sum <-
## mod_stan_4$draws(variables = "cellmeans") |>
mod_stan_4$draws(variables = "Years") |>
posterior::as_draws_df() |>
mutate(across(starts_with("Years"), plogis)) |>
posterior::summarise_draws(
median,
HDInterval::hdi,
~ HDInterval::hdi(., credMass = c(0.9)),
rhat,
ess_bulk,
ess_tail
) |>
rename(lower_90 = V4, upper_90 = V5) |>
bind_cols(newdata_4a) |>
mutate(Year = as.numeric(as.character(fYear)))
saveRDS(stan_4_sum,
file = paste0(data_path, "synthetic/stan_4_sum.rds")
) Show the code
stan_4_sum <- readRDS(
file = paste0(data_path, "synthetic/stan_4_sum.rds")
)
all_years <- stan_4_sum |>
reframe(Year = full_seq(Year, period = 1))
stan_4_sum_sam_yrs <- newdata_4 |>
left_join(stan_4_sum, by = "fYear")
gap_years <- stan_4_sum_sam_yrs |>
reframe(Year = setdiff(full_seq(Year, period = 1), Year))
gap_years <- gap_years |>
reframe(Year = full_seq(range(Year) + c(-1, 1), period = 1))
stan_4_sum_gap_yrs <- stan_4_sum |>
right_join(gap_years, by = "Year")
stan_4_sum_sam_yrs <- stan_4_sum_sam_yrs |>
full_join(all_years, by = "Year")
g1 <- stan_4_sum_sam_yrs |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "stan")) +
geom_ribbon(data = stan_4_sum_gap_yrs,
aes(x = Year, ymin = lower, ymax = upper), alpha = 0.1) +
geom_line(data = stan_4_sum_gap_yrs,
aes(x = Year, y = median, color = "stan"), alpha = 0.4) +
geom_line(data = mod_simple_4,
aes(x = Year, y = Mean, colour = "simple data mean"), linetype = "dashed") +
geom_line(data = mod_simple_4,
aes(x = Year, y = Median, colour = "simple data median"), linetype = "dashed") +
theme_bw()
g2 <- stan_4_sum_sam_yrs |>
ggplot() +
geom_ribbon(aes(x = Year, ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(aes(x = Year, y = median, color = "stan")) +
geom_ribbon(data = stan_4_sum_gap_yrs,
aes(x = Year, ymin = lower, ymax = upper), alpha = 0.1) +
geom_line(data = stan_4_sum_gap_yrs,
aes(x = Year, y = median, color = "stan"), alpha = 0.4) +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Mean, colour = "all mean"), linetype = "dashed") +
geom_line(data = benthos_reefs_temporal_summary,
aes(x = Year, y = Median, colour = "all median"), linetype = "dashed") +
geom_line(data = all_sampled_sum,
aes(x = Year, y = response, colour = type), linetype = "dashed") +
theme_bw()
ggsave(
filename = paste0(
fig_path, "R_pdp_mod_stan_4.png"
),
g1 + g2,
width = 12, height = 6, dpi = 72
)
- left figure compares brms modelled estimates to the sample data
- right figure compares brms modelled estimates to the simple hierarchical sampled means/medians and the true means/medians
- sampled reefs (blue and green dashed lines) are lower than the true global average (orange and khaki) - right side figure
- brms modelled estimates are slightly higher than the 5 site, replacement data, yet biased lower than the “truth”
- the uncertainty is relatively high, but not as high as brms or glmmTMB
Show the code
mod_stan_4 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_4.rds")
)
color_scheme_set("viridis")
g <-
mod_stan_4$draws(variables = c("beta", "phi", "sd_1", "sd_2", "sd_3")) |>
mcmc_trace() +
theme_minimal()
ggsave(
filename = paste0(
fig_path, "R_stan_trace_4.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
mod_stan_4 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_4.rds")
)
color_scheme_set("viridis")
g <-
mod_stan_4$draws(variables = c("beta", "phi", "sd_1", "sd_2", "sd_3")) |>
mcmc_acf() +
theme_minimal()
ggsave(
filename = paste0(
fig_path, "R_stan_ac_4.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
mod_stan_4 <- readRDS(
file = paste0(data_path, "synthetic/mod_stan_4.rds")
)
g <-
bayesplot::pp_check(
benthos_fixed_locs_obs_4$cover,
mod_stan_4$draws("ypred", format = "matrix")[1:100, ],
ppc_dens_overlay
) +
theme_classic()
ggsave(
filename = paste0(
fig_path, "R_stan_ppc_4.png"
),
g,
width = 10, height = 8, dpi = 72
)
Show the code
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
phi 312.18 312.27 10.07 10.54 295.93 328.76 1.00 2674 2868
beta[1] -1.26 -1.26 0.05 0.05 -1.34 -1.18 1.00 1509 2227
beta[2] -1.24 -1.24 0.05 0.05 -1.32 -1.16 1.00 1506 2095
beta[3] -1.07 -1.07 0.05 0.05 -1.15 -0.99 1.00 1491 2235
beta[4] -1.04 -1.04 0.05 0.05 -1.12 -0.96 1.00 1489 2327
beta[5] -1.07 -1.07 0.05 0.05 -1.16 -0.99 1.00 1479 2274
beta[6] -1.22 -1.22 0.05 0.05 -1.31 -1.14 1.00 1485 2123
beta[7] -1.32 -1.32 0.05 0.05 -1.40 -1.24 1.00 1526 2225
beta[8] -1.28 -1.28 0.05 0.05 -1.36 -1.20 1.00 1522 2285
beta[9] -1.28 -1.28 0.05 0.05 -1.36 -1.20 1.00 1515 2431
sd_2[1] 0.09 0.09 0.03 0.02 0.05 0.14 1.00 1363 1104
sd_3[1] 0.17 0.17 0.01 0.01 0.16 0.19 1.00 2254 2608 variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
Years[1] -1.24 -1.24 0.05 0.05 -1.33 -1.16 1.00 2038 2309
Years[2] -1.24 -1.24 0.06 0.06 -1.34 -1.14 1.00 1927 2459
Years[3] -1.18 -1.17 0.15 0.15 -1.42 -0.93 1.00 2535 2773
Years[4] -1.12 -1.12 0.15 0.15 -1.36 -0.88 1.00 2837 2739
Years[5] -1.08 -1.08 0.09 0.09 -1.23 -0.94 1.00 2324 2461
Years[6] -1.08 -1.07 0.07 0.07 -1.19 -0.96 1.00 2045 2344
Years[7] -1.08 -1.08 0.07 0.06 -1.18 -0.97 1.00 1893 2496
Years[8] -1.12 -1.11 0.06 0.06 -1.22 -1.01 1.00 2050 2352
Years[9] -1.20 -1.20 0.06 0.06 -1.31 -1.10 1.00 1954 2664
Years[10] -1.27 -1.26 0.07 0.06 -1.37 -1.16 1.00 1912 2720
Years[11] -1.27 -1.27 0.07 0.06 -1.38 -1.16 1.00 1923 2571
Years[12] -1.27 -1.27 0.07 0.07 -1.39 -1.16 1.00 2007 2520