Correct log likelihood formulation for non-centered model with response uncertainty

Jason_Mercer · February 8, 2021, 11:13pm

Hi,
I’m trying to figure out the appropriate log likelihood formulation for a model in which the response is uncertain. Further, I’m using a non-centered parameterization to estimate the “true” value of the response as my parameter estimates using the centered formulation are expressing a form of Neal’s funnel.

The reason I ask is because my impression is that if I have a normal log likelihood, I should calculate that likelihood using the observed response (y_meas), the expected value (mu) and a scaling parameter (sigma), whatever the expected “true” value (y_true) is determined from the model, accounting for measurement uncertainty in the response (y_noise). However, when I look at the distribution of log likelihoods for the results (or try to run them in the {loo} package), I do not get values that make any sense. The only way I get likelihoods that seem reasonable is to compare y_meas to y_true and y_noise.

Below is a toy linear model highlighting what I’m talking about. I have three scenarios. A) how I thought I was supposed to calculate the log likelihood in Stan, B) an alternative that includes the scaling factor on sigma (helps prevent the “funnel”) – may be identical to A, but more explicit, and C) the parameterization that seems to produce the most reasonable results, but I don’t know why.

This is just a toy example (e.g., these data don’t necessitate a non-centered parameterization), but here is some code to produce some data, focusing a linear model for the data generation process:

library(rstan)
library(tidybayes)
library(tidyverse)

set.seed(20821)
#Number of samples
n <- 30
#Some explanatory
x <- rnorm(n, mean = 0, sd = 5)
#Intercept
alpha <- 0.5
#Slope
beta <- 2.2
#Relationship uncertainty
epsilon <- rnorm(n = n, mean = 0, sd = 1)
#Perfect fit between response and explanatory
y_perfect <- alpha + beta * x
#Fit after including relationship uncertainty
y_true <- alpha + beta * x + epsilon
#Measurement noise (equivalent to standard error)
y_noise <- rgamma(n, shape = 10)
#Measurement error (and fit uncertainty)
y_meas <- rnorm(n, mean = y_true, sd = y_noise)

#Model data
stan_data_uncertain <- list(N = n,
  y_meas = y_meas,
  y_noise = y_noise,
  x = x)

And the Stan code (with the generated quantities block highlighting the three scenarios):

data {
  int<lower=0> N; //Number of observations
  vector[N] y_meas; //Observed response
  vector<lower=0>[N] y_noise; //Standard error of response
  vector[N] x; //Explanatory matrix
}

parameters {
  real alpha; //Intercept
  real beta; //Slope
  real<lower=0> z_sigma; //Scales non-centered standard deviation
  real<lower=0> sigma; //Standard deviation
}

transformed parameters {
  vector[N] mu; //Expected response from model
  vector[N] y_true; //The "true" response value
  mu = alpha + beta * x; //Model
  y_true = mu + z_sigma * sigma; //Non-centered estimate of y_true
}

model {
  target += std_normal_lpdf(alpha); //Standard normal intercept prior
  target += std_normal_lpdf(beta); //Standard normal slope prior
  target += std_normal_lpdf(z_sigma); //Non-centered scaling factor for sigma
  target += exponential_lpdf(sigma | 1); //Exponential prior for sigma
  target += normal_lpdf(y_meas | y_true, y_noise); //Likelihood using latent
//variable and measurement uncertainty
}

generated quantities {
  vector[N] pred_A; //Predictions from scenario A
  vector[N] pred_B; //Predictions from scenario B
  vector[N] pred_C; //Predictions from scenario C
  vector[N] log_lik_A; //Log likelihood for scenario A
  vector[N] log_lik_B; //Log likelihood for scenario B
  vector[N] log_lik_C; //Log likelihood for scenario C
  
  for(n in 1:N) {
    pred_A[n] = normal_rng(mu[n], sigma);
    pred_B[n] = normal_rng(mu[n], z_sigma * sigma);
    pred_C[n] = normal_rng(y_true[n], y_noise[n]);
    log_lik_A[n] = normal_lpdf(y_meas[n] | mu[n], sigma);
    log_lik_B[n] = normal_lpdf(y_meas[n] | mu[n], z_sigma * sigma);
    log_lik_C[n] = normal_lpdf(y_meas[n] | y_true[n], y_noise[n]);
  }
}

Here is a visualization (and code; Stan model results are stored in model_samples_uncertain) of the log likelihoods (note the different axis scales):

##Log likelihoods
gather_draws(model_samples_uncertain, log_lik_A[n], log_lik_B[n], log_lik_C[n]) %>%
  ggplot(aes(.value)) +
  geom_density() +
  labs(x = NULL) +
  facet_wrap(~ .variable, scales = "free", strip.position = "bottom") +
  theme_bw() +
  theme(strip.background = element_blank(),
    strip.placement = "outside",
    panel.grid.minor = element_blank()) +
  scale_fill_brewer(palette = "Set1")

The first two scenarios seem pretty extreme, while the scenario that estimates the “true” response value and predetermined uncertainty look reasonable. Just to help confirm these whacky values, below are the predictions across scenarios, compared to the measured response values.

##Response predictions
gather_draws(model_samples_uncertain, pred_A[n], pred_B[n], pred_C[n]) %>%
  group_by(n, .variable) %>%
  summarise(y_pred = quantile(.value, probs = 0.5),
    y_lo = quantile(.value, probs = 0.05),
    y_hi = quantile(.value, probs = 0.95)) %>%
  full_join(tibble(y_meas) %>%
      mutate(n = 1:n())) %>%
  ggplot(aes(y_pred, y_meas)) +
  geom_abline(intercept = 0, slope = 1) +
  geom_errorbar(aes(ymin = y_lo, ymax = y_hi)) +
  geom_point(size = 3, shape = 21, fill = "gray", color = "black") +
  labs(y = "Data generating response",
    x = "Predicted response") +
  facet_wrap(~ .variable) +
  theme_bw() +
  theme(strip.background = element_blank())

To be clear, the x-axis represents the predicted values from the three different scenarios, the y-axis is the data used to fit the model, the vertical lines are the 90th percent credible intervals, and the lines bisecting the plots are the 1:1 line. I’m not actually sure why the model is so confident about the results in the first two scenarios, but if we compare them to the “prefect” response values (no relationship or measurement uncertainties), the results actually make a lot of sense:

gather_draws(model_samples_uncertain, pred_A[n], pred_B[n], pred_C[n]) %>%
  group_by(n, .variable) %>%
  summarise(y_pred = quantile(.value, probs = 0.5),
    y_lo = quantile(.value, probs = 0.05),
    y_hi = quantile(.value, probs = 0.95)) %>%
  full_join(tibble(y_perfect) %>%
      mutate(n = 1:n())) %>%
  ggplot(aes(y_pred, y_perfect)) +
  geom_abline(intercept = 0, slope = 1) +
  geom_errorbar(aes(ymin = y_lo, ymax = y_hi)) +
  geom_point(size = 3, shape = 21, fill = "gray", color = "black") +
  labs(y = "Data generating response",
    x = "Predicted response") +
  facet_wrap(~ .variable) +
  theme_bw() +
  theme(strip.background = element_blank())

So, I’m basically confused as to what I’m supposed to be comparing to what now.

Does anyone know how one is actually supposed to calculate the log likelihoods for a model like this?

bbbales2 · February 9, 2021, 1:04pm

I think a linear regression will work here:

data {
  int<lower=0> N; //Number of observations
  vector[N] y_meas; //Observed response
  vector<lower=0>[N] y_noise; //Standard error of response
  vector[N] x; //Explanatory matrix
}

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

transformed parameters {
  vector[N] mu; //Expected response from model
  mu = alpha + beta * x; //Model
}

model {
  target += std_normal_lpdf(alpha); //Standard normal intercept prior
  target += std_normal_lpdf(beta); //Standard normal slope prior
  target += exponential_lpdf(sigma | 1); //Exponential prior for sigma
  target += normal_lpdf(y_meas | mu, y_noise + sigma); //Likelihood using latent
}

The way I read what you wrote was:

y_{meas, i} \sim N(\mu + \phi, y_{noise, i}) \\ \phi \sim N(0, \sigma) \\ y_{noise, i} \sim Gamma(\alpha, \beta)

You’re passing in y_{noise, i} as data, so that’s known. Because \phi is a normal, this is the same as saying:

y_{meas, i} \sim N(\mu, \sigma + y_{noise, i})

This should be what you get from integrating \phi out of the likelihood.

Jason_Mercer · February 9, 2021, 7:28pm

Oh wow. That’s a really interesting way of approaching the problem I’ve not seen before, but I really like it. And both the approach and results make sense, which is great! But I also like that inclusion of a gamma distribution when the measurement uncertainty isn’t known.

Building on the code I’ve already posted, here’s my take on the code:

data {
  int<lower=0> N; //Number of observations
  vector[N] y_meas; //Observed response
  vector<lower=0>[N] y_noise; //Standard error of response
  vector[N] x; //Explanatory matrix
}

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

transformed parameters {
  vector[N] mu;
  mu = alpha + beta * x;
}

model {
  //Priors
  target += std_normal_lpdf(alpha);
  target += std_normal_lpdf(beta);
  target += exponential_lpdf(sigma | 1);
  //Likelihood
  target += normal_lpdf(y_meas | mu, y_noise + sigma);
}

generated quantities {
  vector[N] pred; //Predictions
  vector[N] log_lik; //Log-likelihood

  for(n in 1:N) {
    pred[n] = normal_rng(mu[n], y_noise[n] + sigma);
    log_lik[n] = normal_lpdf(y_meas[n] | mu[n], y_noise[n] + sigma);
  }
}

The resultant samples (stored in model_samples – code not shown) do good job describing both the measured data, as well as the “true” data that was used to draw the measured samples, which is real neat.

gather_draws(model_samples, pred[n]) %>%
  group_by(n) %>%
  summarise(y_pred = quantile(.value, probs = 0.5),
    y_lo = quantile(.value, probs = 0.05),
    y_hi = quantile(.value, probs = 0.95)) %>%
  full_join(tibble(Measured = y_meas) %>%
      mutate(n = 1:n())) %>%
  full_join(tibble(True = y_true) %>%
      mutate(n = 1:n())) %>%
  pivot_longer(Measured:True, names_to = "source") %>%
  ggplot(aes(y_pred, value)) +
  geom_abline(intercept = 0, slope = 1) +
  geom_errorbar(aes(ymin = y_lo, ymax = y_hi)) +
  geom_point(size = 3, shape = 21, fill = "gray", color = "black") +
  facet_wrap(~ source, strip.position = "left") +
  labs(y = NULL,
    x = "Predicted response",
    title = "Measured (relationship + measurement uncertainty) vs True (relationship uncertainty only)") +
  theme_bw() +
  theme(strip.background = element_blank(),
    strip.placement = "outside")

Also notable (code not provided) is that this version of the code does marginally better than a version of the model that does not include uncertainty (using LOOIC from the {loo} package, which is what I would expect.

Thank you @bbbales2!

andrjohns · February 9, 2021, 11:37pm

You can also optimise this a little with the normal_id_glm distribution (linear regression), specified like so:

 target += normal_id_glm_lpdf(y_meas | x, alpha, beta, y_noise + sigma);

You’ll just need to change the declaration of x to matrix[N,1] and beta to vector[1]

Jason_Mercer · February 17, 2021, 4:04am

Just came across that distribution while reading through some other code. Neat solution! Thank you.

Topic		Replies	Views
Non centered parameterization on variance parameter Modeling	31	8459	October 22, 2018
Centered vs non-centered model, different loglikelihoods Modeling specification	6	514	May 15, 2024
Non centered parametrization of lognormal Distribution Modeling specification	4	1611	March 13, 2021
Non-centered parametrisation for hierarchical log_mix (mixed) distribution Modeling rstan , techniques	4	946	February 21, 2022
Non-centered parameterization for likelihood Modeling fitting-issues	2	390	April 25, 2023

Correct log likelihood formulation for non-centered model with response uncertainty

Related topics