Consider a linear, random intercept model:

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


$$ \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 $$


$$ \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.

  • $\begingroup$ Incidentally, I can share specific R code if helpful, but it seems to be software agnostic as I tried it through a couple different packages and ended up in the same place. $\endgroup$
    – Joshua
    Commented Oct 26, 2021 at 0:07
  • $\begingroup$ It would be much better if you have shared your code. $\endgroup$
    – Amin Shn
    Commented Oct 26, 2021 at 8:53
  • $\begingroup$ @AminShn sure here is a link to a minimally reproducible R examples using simulated data gist.github.com/JWiley/31975ca598abe4ab249a1cb19b3ea53d $\endgroup$
    – Joshua
    Commented Oct 26, 2021 at 9:01
  • $\begingroup$ "but $\gamma_0$ consistently has a wider uncertainty interval than the average of $\mu_i$". Is this definitely the right way around? $\endgroup$
    – Eoin
    Commented Oct 26, 2021 at 11:37
  • $\begingroup$ @Eoin sorry that was not well, written, I have edited it to: "γ0 consistently has a wider uncertainty interval than the uncertainty interval for the average of μi," $\endgroup$
    – Joshua
    Commented Oct 26, 2021 at 22:29

2 Answers 2


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$.

  • $\begingroup$ I had always thought of higher level parameters as more certainty since "noise" would have been removed lower down, but they are further from data so the information stream analogy makes sense. $\endgroup$
    – Joshua
    Commented Oct 26, 2021 at 23:00

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:

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

## make it reproducible

## 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 <- "
  # 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))


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
  • $\begingroup$ I think one difference is that I am interested in mubar, not mu itself. mu itself definitely has more variance, but mubar does not (in the examples I've worked with at least). $\endgroup$
    – Joshua
    Commented Oct 26, 2021 at 22:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.