I have fitted the same model with and without error measurement. What I got is a not different posterior predictins of the outcome by season and location. Isnâ€™t it quite surprising, Angelos? I am trying to figure out , why.

This is the code:

```
# brms formula with error measurement
bf_mi <- bf(DOC | mi(0.2) ~ 0 + season:location + (0 + season:location| site)) + lognormal(link = "identity")
# brms formula without error measurement
bf <- bf(DOC ~ 0 + season:location + (0 + season:location| site)) + lognormal(link = "identity")
# Set priors
priors_informative <- c(
#~~~~~~~~~~~ the average summer DOC concentration downstream ~~~~~~~~~
set_prior("normal(0.5, 1.3)", class = "b", coef = "seasonsummer:locationout" ),
#~~~~~~~~~~~ the average summer DOC concentration upstream ~~~~~~~~~
set_prior("normal(0.5, 1.3)", class = "b", coef = "seasonsummer:locationin"),
#~~~~~~~~~~~ the average winter DOC concentration upstream ~~~~~~~~~
set_prior("normal(0, 1.2)", class = "b", coef = "seasonwinter:locationin"),
#~~~~~~~~~~~ the average winter DOC concentration downstream ~~~~~~~~~
set_prior("normal(0, 1.2)", class = "b", coef = "seasonwinter:locationout") ,
# between-group (here, between-site) variability.
set_prior("normal(0, 0.2)", class = "sd", group= "site", coef = "seasonsummer:locationin") ,
set_prior("normal(0, 0.2)", class = "sd", group= "site", coef = "seasonsummer:locationout") ,
set_prior("normal(0, 0.1)", class = "sd", group= "site", coef = "seasonwinter:locationin") ,
set_prior("normal(0, 0.1)", class = "sd", group= "site", coef = "seasonwinter:locationout") ,
#~~~~~~~~~~~ the residual standard deviation ~~~~~~~~~
set_prior("student_t(3, 0, 0.5)", class = "sigma" )
)
# then fit the models
b1 <- brm(formula = bf_mi,
data = d4, # observed data
prior = priors_informative,
iter = 4000,
warmup = 2000, # number of warmup iterations
chains = 4,
cores = 4,
control = list(adapt_delta = 0.99, max_treedepth = 15), # facilitate convergence
backend = "cmdstanr",
save_pars = save_pars(latent = TRUE),
seed = 7, # allow exact replication of results
sample_prior = TRUE, # sample also the priors
)
b1_no_error <- brm(formula = bf,
data = d4, # observed data
prior = priors_informative,
iter = 4000,
warmup = 2000, # number of warmup iterations
chains = 4,
cores = 4,
control = list(adapt_delta = 0.99, max_treedepth = 15), # facilitate convergence
backend = "cmdstanr",
save_pars = save_pars(latent = TRUE),
seed = 7, # allow exact replication of results
sample_prior = TRUE, # sample also the priors
)
# generate new data frame
conditions <- expand.grid(season = c("summer","winter"), location = c("out","in"), site = unique(d4$site))
# predicted draws with error measurament
f_error <-
add_predicted_draws(b1, ndraws = 1500 , newdata = conditions, scale = "response") %>%
as_tibble() |>
mutate(model = "with error mi(0.2)")
# with no error measurament
f_no_error <-
add_predicted_draws(b1_no_error, ndraws = 1500, newdata = conditions, scale = "response") %>%
as_tibble() |>
mutate(model = "no error")
# bind df
dmerg <- rbind(f_no_error, f_error)
# Plot
ggplot(dmerg, aes(x = .prediction)) +
stat_slab(aes(fill = model), alpha = .53) +
stat_pointinterval(aes(color = model), position = position_dodge(width = .4, preserve = "single"), .width = c(0.66, 0.95)) +
facet_wrap(season ~ location, scales = "free_x", axes = "all", axis.labels = "all_x") +
labs(x = "DOC (mg/L)", y = NULL)+
theme_bw()+
facetted_pos_scales(x = list(
location == "in" & season == "summer" ~ scale_x_continuous(breaks = seq(from = 0, to = 10, by = 1), limits= c(0, 10)),
location == "in" & season == "winter" ~ scale_x_continuous(breaks = seq(from = 0, to = 7, by = 1), limits = c(0, 7)),
location == "out" & season == "summer" ~ scale_x_continuous(breaks = seq(from = 0, to = 12, by = 1), limits = c(0, 12)),
location == "out" & season == "winter" ~ scale_x_continuous(breaks = seq(from = 0, to = 10, by = 1), limits = c(0 , 10))
))
```