Propagation of Measurement Error on the outcome within mi()

Modeling
3

Dear Angelos,
Thank you for the advise. Cleaver and elegant. I have some questions:

  1. Why do you suggest to include as random effect the “site”? Are you proposing to include “site” level because the model needs to know how to pair the DOC_{upstream} and DOC_{downstream} observations?
    In fact, we sampled one in summer and one in winter for the same site. I am attaching part of the data. Please see the plot too.

    Rplot
    Rplot967×566 11.5 KB

  2. Following your advice, I am modelling the DOC_{upstream} and DOC_{downstream} observations using the log-normal family. This is because DOC is continuous positive outcome. So I consider this model:

\begin{align*}\color{red}{\text{DOC}_{\text{OBS}, i}} & \color{red}\sim \color{red}{\operatorname{Log-Normal}(\text{DOC}_{\text{TRUE}, i}, \text{DOC}_{\text{SE}, i})} \\\color{red}{\text{DOC}_{\text{TRUE}, i}} & \sim \operatorname{Log-Normal}(\mu_i, \sigma) \\\mu_i & = e^{\alpha_{\text{location}[season]} } \\\alpha_{\text{upstream}, [summer]} & \sim \operatorname{Normal}(0.5, 1.3) \\\alpha_{\text{upstream}[winter]} & \sim \operatorname{Normal}(0, 1.2) \\\alpha_{\text{downstream}, [summer]} & \sim \operatorname{Normal}(0.5, 1.3) \\\alpha_{\text{downstream}, [winter]} & \sim \operatorname{Normal}(0, 1.2) \\\sigma & \sim \operatorname{Student-t}(3, 0, 0.5).\end{align*}

And the code is

bf_mi <-  bf(DOC | mi(0.2)  ~ 0 + season:location) + lognormal(link = "identity")

priors <- 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") ,
#~~~~~~~~~~~ the  residual standard deviation ~~~~~~~~~
set_prior("student_t(3, 0, 0.5)", class = "sigma" ))

# fit the model 

b1 <- brm(formula = bf_mi,
        data = d4, # observed data 
        prior = priors,
        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 
        file = here("results","Leonardo","models_fit","DOC","b1"),
        sample_prior = TRUE, # sample also the priors 
         )

The data:
d4.csv (21.5 KB)