why do average fitted values have less uncertainty than the intercept in a Bayesian random effects model?

Question

Consider a linear, random intercept model:

$$ y_{ij} \sim \mathcal{N}(\mu_i, \sigma_e^2) $$

where

$$ \mu_i \sim \mathcal{N}(\gamma_0, \sigma_{mu}^2) $$

Estimated in a Bayesian framework using MCMC sampling, I would expect that extracting fitted values for each participant in sample for each MCMC draw of the posterior would give me:

$$ \mu_{i,draw} $$

If I were to average across the $i$ participants, I would get samples from a posterior distribution like:

$$ \bar{\mu}_{draw} $$

which I would expect to have the same value and same uncertainty as $\gamma_0$.

In most software, the model is separated as:

$$ \mu_i = \gamma_0 + u_i $$

where

$$ \gamma_0 \sim \mathcal{N}(0, \sigma_{\mu}^2) $$

When I actually try this, my observation is that the posterior draws for $\gamma_0$ and the average of the fitted values $\mu_i$ have very similar average values, but $\gamma_0$ consistently has a wider uncertainty interval than the uncertainty interval for the average of $\mu_i$, the fitted values for each participant.

Now why I care about this. In the linear model case, I could just use the parameter estimates for $\gamma_0$. However, when fitting say a logistic mixed effects model, the estimate for $\gamma_0$ is on the log odds scale and back transformed to a probability is not the same as getting the averaged probability taking into account the random intercept.

My best guess is that in taking the fitted values:

$$ \hat{\mu_i} = \gamma_0 + u_i $$

and then averaging across all participants and then summarizing that posterior distribution is somehow missing a source of variance, but I've been over this a dozen times and tried to look for articles or books discussing it and have come up empty, so am hoping for some insight here.

Answer 1

In that bayesian model, $\gamma_0$ posterior depends on $\vec\mu$, while $\vec\mu$ posterior depends on data (you may think of it as in a DAG). $\gamma_0$ is conditionally independent on data (conditioned on $\vec\mu$).

In particular, drafted $\gamma_0$ will vary around $\bar\mu$, and this variation is independent on data, so it just sums up to $\bar\mu$'s variance.

You may think it this way: information streams up from data to $\vec\mu$ (low variance), and then to $\gamma_0$ (higher variance), because you set a low informative prior. If you set a prior more informative than data, $\gamma_0$ would have been less variable than $\bar\mu$.

Answer 2

Ok this is the model you're fitting in Stan:

// generated with brms 2.15.0
functions {
}
data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  // data for group-level effects of ID 1
  int<lower=1> N_1;  // number of grouping levels
  int<lower=1> M_1;  // number of coefficients per level
  int<lower=1> J_1[N];  // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_1;
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> sigma;  // residual SD
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  vector[N_1] z_1[M_1];  // standardized group-level effects
}
transformed parameters {
  vector[N_1] r_1_1;  // actual group-level effects
  r_1_1 = (sd_1[1] * (z_1[1]));
}
model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = Intercept + rep_vector(0.0, N);
    
    for (n in 1:N) {
      // add more terms to the linear predictor
      mu[n] += r_1_1[J_1[n]] * Z_1_1[n];
    }
    
    target += normal_lpdf(Y | mu, sigma);
  }
  // priors including constants
  target += student_t_lpdf(Intercept | 3, 48.4, 8.3);
  target += student_t_lpdf(sigma | 3, 0, 8.3)
    - 1 * student_t_lccdf(0 | 3, 0, 8.3);
  target += student_t_lpdf(sd_1 | 3, 0, 8.3)
    - 1 * student_t_lccdf(0 | 3, 0, 8.3);
  target += std_normal_lpdf(z_1[1]);
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
}

I am not a fluent Stan user but I am not sure the model you described above matches this Stan code. Anyway, back to the math you showed us above, $Var(\gamma_{0})$ can't be greater than $Var(\mu)$ because the variance of the former would be additional to the variance of the latter. So maybe double check your code or revise the math. The math you described above can be fitted in JAGS this way:

library(R2jags)
#### Simulate some data ####
nID <- 100 ## number of people
k <- 4 ## number of observations per person

## make it reproducible
set.seed(1234)

## generate random intercepts
mu <- rnorm(nID, mean = 50, sd = 6)

## generate random observations -- 4 per person
y <- rnorm(nID * k, mean = rep(mu, each = k), sd = 6)

## put it all together
d <- data.frame(y = y, ID = rep(seq_len(nID), each = k))

# Jags code to fit the model to the simulated data
model_code <- "
model
{
  # Likelihood
  for (n in 1:N) {
    y[n] ~ dnorm(mu[n], sigma_e^-2)
    mu[n] ~ dnorm(gamma_0, sigma_mu^-2)
    
  }
  # Priors
  gamma_0 ~ dnorm(0.0,0.01)
  sigma_e ~ dunif(0, 100)
  sigma_mu ~ dunif(0, 100)
}
"

# Set up the data
model_data <- list(N = T, y = d$y)

# Choose the parameters to watch
model_parameters <- c("mu", "gamma_0")

# Run the model
model_run <- jags(
  data = model_data,
  parameters.to.save = model_parameters,
  model.file = textConnection(model_code))

print(model_run)

And here are the results which clearly show that $\mu$ has wider CI compared to the $\gamma$:

> print(model_run)
Inference for Bugs model at "6", fit using jags,
 3 chains, each with 2000 iterations (first 1000 discarded)
 n.sims = 3000 iterations saved
         mu.vect sd.vect    2.5%    25%    50%    75%  97.5%  Rhat n.eff
gamma_0    1.589   9.825 -17.581 -4.926  1.775  7.975 21.384 1.001  2000
mu        23.318  35.147 -50.320  1.920 23.662 45.349 94.786 1.004  2900