R code to simulate and fit a random intercept linear model showing the fixed intercept and average of fitted values credible intervals do not align.

fittedvalues.R

library(brms)
library(bayestestR)

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

## subsample so unequal sample sizes by ID
d <- d[sample(nrow(d), size = nrow(d) * .8), ]

m <- brm(y ~ 1 + (1 | ID),
         data = d, chains = 2, cores = 2,
         iter = 4000, refresh = 0, seed = 1234)
msum <- summary(m)

## fitted values ignoring random effects
fitted.fixed.all <- fitted(m, summary = FALSE, re_formula = NA)

## fitted values including random effects
fitted.random.all <- fitted(m, summary = FALSE)

## fitted values including random effects but only one obs per ID so equal weighting
fitted.random.noduplicates <- fitted(m, summary = FALSE)[, !duplicated(model.frame(m)$ID)]

## fixed effect estimate
msum$fixed

## fitted values no random effect summary (matches fixed effect estimate)
bayestestR::describe_posterior(
  rowMeans(fitted.fixed.all),
  centrality = "mean", ci_method = "eti")

## random effects, unequal weighting (all obs included)
## similar point estimate, smaller CI than fixed
bayestestR::describe_posterior(
  rowMeans(fitted.random.all),
  centrality = "mean", ci_method = "eti")

## random effects, equal weighting (all obs included)
## similar point estimate, smaller CI than fixed
bayestestR::describe_posterior(
  rowMeans(fitted.random.noduplicates),
  centrality = "mean", ci_method = "eti")