# Bafumi-Gelman correction for correlation between fixed and random effects

**URL:** <https://discourse.mc-stan.org/t/bafumi-gelman-correction-for-correlation-between-fixed-and-random-effects/39185>\
**Category:** Modeling\
**Created:** [March 30, 2025, 9:25pm UTC](https://discourse.mc-stan.org/t/bafumi-gelman-correction-for-correlation-between-fixed-and-random-effects/39185 "2025-03-30T21:25:39Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Caio\_Simon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/caio_simon/32/19634_2.png) [@Caio\_Simon](https://discourse.mc-stan.org/u/Caio_Simon)\
**Post date:** [March 30, 2025, 9:25pm UTC](https://discourse.mc-stan.org/t/bafumi-gelman-correction-for-correlation-between-fixed-and-random-effects/39185/1 "2025-03-30T21:25:39Z")

</div>

Hello,

I was reading a [paper](https://sites.stat.columbia.edu/gelman/research/unpublished/Bafumi_Gelman_Midwest06.pdf#page=3.10) by Bafumi and Gelman where they showcase a simple method for adjusting for correlations between predictors and group-level effects. I tried to implement this in a model I am currently working on for school. The objective is basically to perform MrP to predict cycling in small areas in Europe. However, when I apply the simple correction, the model takes incredibly long to fit (around 3 hours on 20 thousand rows; when I do it without the correction the time goes to around 40 minutes to 1 hour), Here is the code:

```no-highlight
data{
  int<lower = 0> N; // Number of rows
  array[N] int<lower = 0, upper = 1> y; // Outcome variable (cycling)
  array[N] int<lower = 1, upper = 7> age_id; // Age ID
  array[N] int<lower = 1, upper = 246> nuts_id; // Unique nuts region identifier
  int<lower=0> K; // Number of fixed effects predictors
  matrix[N, K] x; // Fixed effects matrix
  matrix[246, K] fe_adj; // [N_regions, N_fixed effects] median fixed effects by region
  matrix[7, K] fe_adj_age; // [N_category, N_fixed effects] median fixed effects by age category
}
parameters{
  real alpha; // Global intercept
  
  // Random effects: Creating scale and intercept parameters
  real<lower = 0> gamma_age_scale;
  vector[7] gamma_age;
  
  real<lower = 0> gamma_region_scale;
  vector[246] gamma_region;
  
  // Fixed effects coefficients
  vector[K] beta;
  vector[K] adj_beta;
  vector[K] adj_beta_age;
}
transformed parameters{
  vector[N] mu;
  
  // Non-centred parametrisation
  vector[7] gamma_age_star = gamma_age * gamma_age_scale + fe_adj_age * adj_beta_age;
  vector[246] gamma_region_star = gamma_region * gamma_region_scale + fe_adj * adj_beta;
  
  mu = alpha + gamma_age_star[age_id] + gamma_region_star[nuts_id] + x * beta;

}
model{
  alpha ~ normal(0, 5);
  
  // Standard random effects
  to_vector(gamma_age) ~ normal(0, 1);
  gamma_age_scale ~ normal(0, 5);
  
  to_vector(gamma_region) ~ normal(0, 1);
  gamma_region_scale ~ normal(0, 5);
  
  beta ~ normal(0, 5);
  adj_beta ~ normal(0, 3);
  adj_beta_age ~ normal(0, 3);
  
  y ~ bernoulli_logit(mu);
}
generated quantities{
  vector[N] log_lik; // log likelihood for model checking and comparison later
    for (n in 1:N) {
        log_lik[n] = bernoulli_lpmf(y[n] | inv_logit(mu[n]));
    }
}

```

I adjust for the fixed effects correlation using the fe\_adj \* fe\_beta line, where fe\_adj is a matrix where for every unique region I put the group-level medians of the fixed effects. Same goes for fe\_adj\_age.

Am I doing something wrong? If so, please tell me .

Kind regards,  
Caio

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [March 31, 2025, 9:22pm UTC](https://discourse.mc-stan.org/t/bafumi-gelman-correction-for-correlation-between-fixed-and-random-effects/39185/2 "2025-03-31T21:22:20Z")

</div>

I don’t know the paper you mention, but we can ask @andrewgelman.

If you declare `mu` last, you can put the declaration nd definition on the same line. Usually we can make these non-centered parameterizations one-liners, but the adjustment you’re adding won’t support that.

You can use `bernoulli_logit_lpmf` in the generated quantities block—it’ll be easier for readers since it follows the original model definition.

---

<div class="post-metadata">

**Author:** ![andrewgelman](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrewgelman/32/1577_2.png) [@andrewgelman](https://discourse.mc-stan.org/u/andrewgelman)\
**Post date:** [March 31, 2025, 9:30pm UTC](https://discourse.mc-stan.org/t/bafumi-gelman-correction-for-correlation-between-fixed-and-random-effects/39185/3 "2025-03-31T21:30:22Z")

</div>

Heres the paper: [https://sites.stat.columbia.edu/gelman/research/unpublished/Bafumi\_Gelman\_Midwest06.pdf](https://sites.stat.columbia.edu/gelman/research/unpublished/Bafumi_Gelman_Midwest06.pdf)

---

<div class="post-metadata">

**Author:** ![StaffanBetner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/staffanbetner/32/11056_2.png) [@StaffanBetner](https://discourse.mc-stan.org/u/StaffanBetner)\
**Post date:** [April 1, 2025, 8:16am UTC](https://discourse.mc-stan.org/t/bafumi-gelman-correction-for-correlation-between-fixed-and-random-effects/39185/4 "2025-04-01T08:16:39Z")

</div>

Just wanted to add a note on the terminology and related approaches. The Bafumi-Gelman correction you’ve implemented is a specific application of what’s more generally known as the **Mundlak-Chamberlain device** or a **Correlated Random Effects (CRE)** model. It’s designed precisely to handle the correlation between predictors and group-level effects that you mentioned.

@ruben has written an excellent post about [this](https://rubenarslan.github.io/posts/2023-03-11-latent-group-mean-centering-revisited), framing it as **latent group-mean centering** , particularly in the context of brms. A key distinction in the “latent” approach he discusses is that the group means (or similar summaries) of the predictors are estimated as parameters within the model itself, rather than being pre-calculated and passed in as data (like your fe\_adj medians).