Hello,
I’ve been trying to wrap my head around the measurement error options in brms for the response variable, and I’m still unsure about the practical differences between modelling a response with measurement error in brms using, e.g. y | mi(sd) vs. y | se(sd, sigma = TRUE).
For example, some sources seem to suggest that they are equivalent (e.g., here and here) in practice, even if the underlying likelihood statements in the stan code are different. In fact, the second link seems to suggest that using y | se(sd, sigma = TRUE) might often bypass convergence issues encountered with y | mi(sd) (and I have experienced this recently).
So, imagine a dataset that contains a continuous response y which was obtained from an instrument that also provides uncertainty for each observation, and a continuous predictor x — for simplicity of this example, let’s assume that it is extremely precise and carries no error:
y | mi(sd)
library(tidyverse)
library(brms)
set.seed(10)
dat <- data.frame(x = rnorm(100, 10, 2)) |>
dplyr::mutate(
y = rnorm(100, -10 + x * 5, 1), sd_y = rgamma(100, 2, 5)
)
brms::make_stancode(y | mi(sd_y) ~ x, data = dat)
In this case in which no missing observations exist in the dataset, the key elements in the stancode would be
data {
....
// data for measurement-error in the response
vector<lower=0>[N] noise;
....
}
....
parameters {
vector[N] Yl; // latent variable
....
}
....
model {
// likelihood including constants
if (!prior_only) {
target += normal_id_glm_lpdf(Yl | Xc, Intercept, b, sigma);
}
// priors including constants
target += lprior;
target += normal_lpdf(Y[Jme] | Yl[Jme], noise[Jme]);
}
....
y | se(sd, sigma = TRUE)
brms::make_stancode(y | se(sd_y, sigma = TRUE) ~ x, data = dat)
data {
....
vector<lower=0>[N] se; // known sampling error
....
}
transformed data {
vector<lower=0>[N] se2 = square(se);
....
}
....
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = rep_vector(0.0, N);
mu += Intercept + Xc * b;
target += normal_lpdf(Y | mu, sqrt(square(sigma) + se2));
}
// priors including constants
target += lprior;
}
....
The first likelihood declaration with a latent step states
Y_l \sim N(\mu, \sigma)
Y \sim N(Y_l, \sigma_{Y})
and the second states
Y \sim N\left(\mu, \sqrt{\sigma^2 + \sigma_{Y}^2}\right)
Noting that se == noise between the two stan codes (i.e. they are the same vector, just named differently depending on the input formula syntax), and borrowing rationale from here, the two statements could be made equivalent if:
U \sim N(0, \sigma_{Y})
Y_l \sim N(\mu, \sigma)
Y = U + Y_l \rightarrow Y \sim N\left(\mu, \sqrt{\sigma^2 + \sigma_{Y}^2}\right)
But I am not sure if there’s some flawed logic in there. I was hoping someone might be able to provide some definitive clarity on the issue, and whether my interpretation of the equivalence is correct.
Thanks!!