Bayesian Measurement Error Model in Stan

Modeling
,
1

Hi everyone,

I’m working on a Bayesian Measurement Error Model in Stan and would appreciate your feedback on whether my code is correctly specified. The objective is to account for measurement errors in both x (Method 1) and y (Method 2) while modeling their true relationship using a Bayesian approach.

  • Method 1 estimates evapotranspiration,
  • Method 2 measures evapotranspiration and is treated as the reference.

Both methods are assumed to have a 10% error. I’d greatly appreciate any insights or suggestions on improving the model specification.
Thank you in advance

Stan code:

data {
int<lower=0> N;                      // Number of data points (days)
vector[N] x_obs;                     // Observed values for x (Method 1)
vector[N] y_obs;                     // Observed values for y (Method 2)
real<lower=0> sigma_x;         // Known measurement error in Method 1 (x) (10%)
real<lower=0> sigma_y;         // Known measurement error in Method 2 (y) (10%)
}

parameters {
real alpha;                              // Intercept
real beta;                                // Slope
real<lower=0> sigma;             // Residual standard deviation
}

model {
// Priors
alpha ~ normal (0, 10);            // Weakly informative prior for intercept
beta ~ normal (0, 10);              // Weakly informative prior for slope
sigma ~ cauchy (0, 5);             // Weakly informative prior for residuals

// Structural model with error propagation for both x and y
y_obs ~ normal(alpha + beta * x_obs, sqrt(sigma^2 + sigma_y^2 + beta^2 * sigma_x^2));
}
2

Hi, @Alessandro, and sorry for missing this earlier.

The usual thing to do here is to treat the true x as an unkown. The measurement error in y is rolled into sigma_y. Then you’d have something like this, where you fit sigma_x, sigma_y and x,

parameters {
  vector[N] x;
  real<lower=0> sigma_x;
  real<lower=0> sigma_y;

model {
  sigma_x ~ normal(0.1, 1);
  sigma_y ~ normal(0.1, 1);
  x ~ normal(x_obs, sigma_x);
  ...
  y_obs ~ normal(alpha + beta * x, sigma_y);
}

The sigma_y variable soaks up modeling error and observation error in y_obs.

You’d want to adjust location and scale of sigma_x and sigma_y.