Measurement Error Model in Variables

Modeling
1

Hi everyone

I am working on a measurement error model in variables for the validation of satellite-derived evapotranspiration (ET) estimates against ground-truth observations collected in perennial crop fields. Validating satellite ET products is inherently challenging because ground-based measurements, typically obtained from eddy covariance stations, are not error-free.

The model I adopt is a Bayesian Simple Linear Regression (SLR) with measurement error in both response and explanatory variables:

  1. Outcome model, relates observed satellite ET to the true latent ET via a Student-t likelihood, providing robustness to outliers.
  2. To account for inter-annual variability in the satellite-ground ET relationship, the model adopts a hierarchical structure over years (K = 3 years). All variables are standardized prior to fitting to improve sampler geometry, and a non-centered parameterization is used for all year-level random effects.

The Stan model is reported below:


data {
  int<lower=0> N;
  int<lower=0> K;                          // Number of years
  array[N] int<lower=1, upper=K> year_id;  // Year index for each observation
  vector[N] x_meas;                 // Ground-station ET (observed, with error)
  vector[N] y_obs;                         // Satellite ET (observed)
  vector[N] sigma_u;                       // Known ground measurement error SD
  real mean_x;
  real sd_x;
  real mean_y;
  real sd_y;
}
transformed data {
  vector[N] x_meas_std  = (x_meas - mean_x) / sd_x;
  vector[N] y_obs_std   = (y_obs  - mean_y) / sd_y;
  vector[N] sigma_u_std = sigma_u / sd_x;
  real slope_prior_center = sd_x / sd_y;  // ~1 on standardised scale when x ~ y
}
parameters {
  // ── Population-level hyperparameters ─────────────────────────────────
  real          mu_beta0;               // population mean intercept
  real<lower=0> tau_beta0;              // population SD intercept
  real          mu_beta1;               // population mean slope
  real<lower=0> tau_beta1;              // population SD slope
  real          mu_log_sigma_e;         // log population mean of residual SD
  real<lower=0> tau_log_sigma_e;        // SD of log residual SD
  real          mu_log_sigma_x;         // log population mean of latent ET SD
  real<lower=0> tau_log_sigma_x;        // SD of log latent ET SD
  real<lower=2> nu;                     // Student-t degrees of freedom (shared)

  // ── Year-level random effects (non-centered) ──────────────────────────
  vector[K] z_beta0;
  vector[K] z_beta1;
  vector[K] z_log_sigma_e;
  vector[K] z_log_sigma_x;

  // ── Latent true ET ───────────────────────────────────────────────────
  vector[N] x_true_raw;                // std_normal auxiliary variable
}
transformed parameters {
  // Year-specific parameters (non-centered reparameterisation)
  vector[K] beta0_std   = mu_beta0 + tau_beta0 * z_beta0;
  vector[K] beta1_std   = mu_beta1 + tau_beta1 * z_beta1;
  vector[K] sigma_e_std = exp(mu_log_sigma_e + tau_log_sigma_e * z_log_sigma_e);
  vector[K] sigma_x_std = exp(mu_log_sigma_x + tau_log_sigma_x * z_log_sigma_x);

  // Latent true ET scaled by year-specific spread
  vector[N] x_true_std;
  for (i in 1:N)
    x_true_std[i] = sigma_x_std[year_id[i]] * x_true_raw[i];
}
model {
  // ── Hyperpriors (Gelman 2006; Vehtari et al. 2021) ───────────────────

  // Intercept: centred at 0 on standardised scale
  mu_beta0        ~ normal(0, 1.0);
  tau_beta0       ~ normal(0, 0.2);    // half-normal (lower=0 constraint)

  // Slope: centred at slope_prior_center (~1 on std scale, expecting sat ≈ ground)
  mu_beta1        ~ normal(slope_prior_center, 0.3);
  tau_beta1       ~ normal(0, 0.1);    // tight: slope not expected to vary strongly across years

  // log(σ_e): log-scale prior, centred around log(0.5) ≈ -0.7 on std scale
  mu_log_sigma_e  ~ normal(-0.7, 0.7);
  tau_log_sigma_e ~ normal(0, 0.15);

  // log(σ_x): log-scale prior for latent ET spread
  mu_log_sigma_x  ~ normal(0, 0.7);
  tau_log_sigma_x ~ normal(0, 0.15);

  // ν (degrees of freedom): gamma(10, 0.5) → mean = 20; allows heavy tails
  nu ~ gamma(10, 0.5);

  // ── Non-centered auxiliary priors ────────────────────────────────────
  z_beta0       ~ std_normal();
  z_beta1       ~ std_normal();
  z_log_sigma_e ~ std_normal();
  z_log_sigma_x ~ std_normal();
  x_true_raw    ~ std_normal();

  // ── Likelihood ────────────────────────────────────────────────────────
  for (i in 1:N) {
    // Measurement model: observed ground ET | true ET ~ Normal(x_true, σ_u²)
    x_meas_std[i] ~ normal(x_true_std[i], sigma_u_std[i]);

    // Outcome model: satellite ET | true ET ~ Student-t(ν, β₀ + β₁·x_true, σ_e)
    y_obs_std[i]  ~ student_t(nu, beta0_std[year_id[i]] + beta1_std[year_id[i]]   * x_true_std[i],sigma_e_std[year_id[i]]);
  }
}
generated quantities {
  // ── Back-transform year-specific parameters to original scale ─────────
  vector[K] beta_1;
  vector[K] beta_0;
  vector[K] sigma_e;
  vector[K] sigma_x;

  for (k in 1:K) {
    beta_1[k] = beta1_std[k] * sd_y / sd_x;
    beta_0[k] = mean_y - beta_1[k] * mean_x + beta0_std[k] * sd_y;
    sigma_e[k] = sigma_e_std[k] * sd_y;
    sigma_x[k] = sigma_x_std[k] * sd_x;
  }

  // ── Population-level parameters (original scale) ──────────────────────
  real pop_beta_1  = mu_beta1 * sd_y / sd_x;
  real pop_beta_0  = mean_y - pop_beta_1 * mean_x + mu_beta0 * sd_y;
  real pop_sigma_e = exp(mu_log_sigma_e) * sd_y;   // population median σ_e
  real pop_sigma_x = exp(mu_log_sigma_x) * sd_x;

  // ── Predictions ───────────────────────────────────────────
  vector[N] y_pred;
  vector[N] y_rep;
  vector[N] log_lik;

  for (i in 1:N) {
    real x_orig = x_true_std[i] * sd_x + mean_x;
    int  k      = year_id[i];
    y_pred[i]  = beta_0[k] + beta_1[k] * x_orig;
    y_rep[i]   = student_t_rng(nu, beta_0[k] + beta_1[k] * x_orig, sigma_e[k]);
    log_lik[i] = student_t_lpdf(y_obs_std[i] | nu,
                                beta0_std[k] + beta1_std[k] * x_true_std[i],
                                sigma_e_std[k]);
  }
}

Could anyone help me check whether the code is correct?Modeling hierarchical-model

2

Does the Stan model compile?

It looks like you have possibly leapt to the final version of your model, instead of starting with a simple model and progressively adding more components. A very important part of this iterative process is confirming that your models are able to estimate the parameters of synthetic datasets. An additional benefit of this approach is that you develop an understanding of the generative model that you are using to approximate the data generation process.

3

Yes, It compiles and gives good results.

4

Claude (?) is pretty good at writing Stan models, but not perfect (the tell here is the block comment style and aligned twiddles—all the Stan code that Claude writes look like that). There’s a lot here that I’d have done differently. Also, I’m not sure we’re going to be able to keep up verifying people writing complicated models with AI.

  1. I would define the means and standard deviation in the transformed data block to prevent user error sending them in.

  2. You can use offset/multiplier to simplify non-centering. Also, non-centering is only good if there’s relatively little data coupled with a weakly informative prior—when there’s a lot of information in data and prior, centered is better.

  3. You can vectorize the likelihood functions, e.g.

x_meas_std ~ normal(x_true_std, sigma_u_std);
  1. You can also vectorize the assignments, so
vector[N] x_true_std = sigma_x_std[year_id] .* x_true_raw;

I don’t like the name “true” in parameters. We never know the true values of anything latent unless we’re simulating.

  1. The variables like beta0_std are misnamed. The variables z_beta0 are the standardized variables, so these should just be called beta0, etc.

  2. You might want to convert the variables back to their natural scale rather than their centered scale.

What you really want to do is run posterior predictive checks to make sure the model’s capturing what you want it to capture.