Same model, better shape: why centering improves MCMC

Log odds ratios under each parameterization

With 0/1 coding, the log-odds ratio for \(A\) alone (that is, when \(B=0\)) is simply \[ \begin{align*} \text{lOR}_a(B=0) &= (\alpha + \beta_a \cdot 1 + \beta_b \cdot 0 + \beta_{ab} \cdot 0) -(\alpha + \beta_a \cdot 0 + \beta_b \cdot 0 + \beta_{ab} \cdot 0) \\ & = \beta_a \end{align*} \] Analogously, the log-odds ratio for \(B\) alone is \(\beta_b\). And if we want to compare the combination of both \(A=1\) and \(B=1\) to the case where neither is activated, then \[ \begin{align*} \text{lOR}_{ab} &= (\alpha + \beta_a \cdot 1 + \beta_b \cdot 1 + \beta_{ab} \cdot 1) \\ &\quad - (\alpha + \beta_a \cdot 0 + \beta_b \cdot 0 + \beta_{ab} \cdot 0) \\ &= \beta_a + \beta_b + \beta_{ab} \end{align*} \] If instead we center the predictors, defining \(A^* = A - 0.5\) and \(B^* = B - 0.5\), then the log-odds ratio of exposure to \(A\) without exposure to \(B\) relative to exposure to neither is \[ \begin{align*} \text{lOR}_a(B=0) &= (\alpha^* + \gamma_a(0.5) + \gamma_b(-0.5) + \gamma_{ab}(0.5)(-0.5)) \\ &\quad - (\alpha^* + \gamma_a(-0.5) + \gamma_b(-0.5) + \gamma_{ab}(-0.5)(-0.5)) \\ &= (\alpha^* + 0.5\gamma_a - 0.5\gamma_b - 0.25\gamma_{ab}) \\ &\quad - (\alpha^* - 0.5\gamma_a - 0.5\gamma_b + 0.25\gamma_{ab}) \\ &= \gamma_a - 0.5\gamma_{ab}. \end{align*} \] Using the same logic we can show that \[ \text{lOR}_{b} = \gamma_{b} - 0.5 \gamma_{ab} \] and \[ \text{lOR}_{ab} = \gamma_a + \gamma_b. \]

Creating a single data set

Here is the data generation process for a single data set. The outcome \(Y\) is generated using the binary parameterization of \(A\) and \(B\):

s_gen <- function(n = 2000,
                    alpha = -0.8,
                    beta_a = 0.5,
                    beta_b = 0.9,
                    beta_ab = -0.3) {
  
  def <- 
    defData(varname = "A", formula = 0.5, dist = "binary") |>
    defData(varname = "B", formula = 0.5, dist = "binary") |>
    defData(varname = "AB", formula = "A*B", dist = "nonrandom") |>
    defData(varname = "A_c", formula = "A - 0.5", dist = "nonrandom") |>
    defData(varname = "B_c", formula = "B - 0.5", dist = "nonrandom") |>
    defData(varname = "AB_c", formula = "A_c * B_c", dist = "nonrandom") |>
    defData(
      varname = "Y", 
      formula = "..alpha + ..beta_a * A + ..beta_b * B + ..beta_ab * AB",
      dist = "binary", link = "logit"
    )
    
  genData(n, def)
  
}

dd <- s_gen()

The two parameterizations fit the same model

First, here is the frequentist check of both models. The fitted probabilities are identical, even though the coefficients differ.

fit_01 <- glm(Y ~ A * B, data = dd, family = binomial)
fit_c  <- glm(Y ~ A_c * B_c, data = dd, family = binomial)

tidy(fit_01)

## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -1.03     0.0995    -10.4  3.90e-25
## 2 A              0.835    0.135       6.18 6.34e-10
## 3 B              1.14     0.134       8.46 2.73e-17
## 4 A:B           -0.670    0.186      -3.61 3.11e- 4

tidy(fit_c)

## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -0.212    0.0464     -4.57 4.91e- 6
## 2 A_c            0.501    0.0929      5.39 7.05e- 8
## 3 B_c            0.802    0.0929      8.63 6.07e-18
## 4 A_c:B_c       -0.670    0.186      -3.61 3.11e- 4

From the 0/1-coded model, \(\text{lOR}_a = 0.835\), \(\text{lOR}_b = 1.14\), and \(\text{lOR}_{ab} = 0.835 + 1.14 - 0.67 = 1.305.\)

From the centered model, \[ \text{lOR}_a = 0.501 + 0.5*0.670 = 0.836 \] \[ \text{lOR}_b = 0.802 + 0.5*0.670 = 1.137 \] \[ \text{lOR}_{ab} = 0.501 + 0.802 = 1.303 \] So the coefficients themselves change under centering, but the underlying treatment contrasts do not.

Specifying the Bayesian models in JAGS

Now we see that we can recover the same treatment contrasts using two different Bayesian models, though computational performance will be improved with centering.

Here is the JAGS code for each model:

model_01 <- "
model {
  for (i in 1:N) {
    Y[i] ~ dbern(p[i])
    logit(p[i]) <- alpha + beta_a * A[i] + beta_b * B[i] + beta_ab * AB[i]
  }
  
  alpha   ~ dnorm(0, 0.25)
  beta_a  ~ dnorm(0, 0.25)
  beta_b  ~ dnorm(0, 0.25)
  beta_ab ~ dnorm(0, 25)
}
"

model_c <- "
model {
  for (i in 1:N) {
    Y[i] ~ dbern(p[i])
    logit(p[i]) <- alpha + gamma_a * A_c[i] + gamma_b * B_c[i] + gamma_ab * AB_c[i]
  }
  
  alpha   ~ dnorm(0, 0.25)
  gamma_a  ~ dnorm(0, 0.25)
  gamma_b  ~ dnorm(0, 0.25)
  gamma_ab ~ dnorm(0, 25)
}
"

Fitting the models

The function fit_jags fits one of the two models just described:

fit_jags <- function(dat, model_string, centered = FALSE,
                     n_chains = 3, burn = 2000, n_iter = 5000) {
  
  # jdat <- as.list(dat[, .(Y, A, B, AB, A_c, B_c, AB_c)])
  if (centered) {
    jdat <- as.list(dat[, .(Y, A_c, B_c, AB_c)])
    vars <- c("alpha", "gamma_a", "gamma_b", "gamma_ab")
  } else {
    jdat <- as.list(dat[, .(Y, A, B, AB)])
    vars <- c("alpha", "beta_a", "beta_b", "beta_ab")
  }
  jdat$N <- nrow(dat)
  
  mod <- jags.model(
    textConnection(model_string),
    data = jdat,
    n.chains = n_chains,
    quiet = TRUE
  )
  
  update(mod, burn, progress.bar = "none")
  
  samp <- coda.samples(
    mod,
    variable.names = vars,
    n.iter = n_iter,
    progress.bar = "none"
  )
  
  samp
}

Now, we can fit the models, collect the diagnostic data, and take a look at the results:

samp_01 <- fit_jags(dd, model_01, centered = FALSE)
samp_c  <- fit_jags(dd, model_c, centered = TRUE)

diag_tbl <- function(samp, model_name) {
  post <- as_draws_df(samp)
  summ <- summarise_draws(post)
  out <- as.data.table(summ)
  out[, model := model_name]
  out[]
}

diag_01 <- diag_tbl(samp_01, "0/1-coded")
diag_c  <- diag_tbl(samp_c, "centered")

Here are the summary statistics of the posterior distribution as well as the computational diagnostics:

There are a few things to notice here. First, the Bayesian estimates for both the 0/1-coded and centered data are closer to zero than the GLM estimates above. The shrinkage is particularly large for the interaction term, because we placed much more restrictive priors on \(\beta_{ab}\) and \(\gamma_{ab}\). This is what we would expect as the prior is pulling the interaction toward zero.

Second, if we compare the two parameterizations, we see that the R-hat—essentially a measure of whether the chains have converged to the same distribution—is slightly lower for the centered data. There isn’t much to make of the difference here (both are very close to 1), but it does suggest slightly more stable behavior for the centered parameterization.

The biggest impact is on the bulk effective sample size (ESS), which reflects how much independent information the chains contain after accounting for autocorrelation. Even though we ran the same number of iterations, the centered model yields far larger ESS values, indicating much better mixing. The sampler is exploring the posterior much more efficiently under the centered parameterization, and in this case the improvement is quite dramatic. Importantly, these differences have nothing to do with the models themselves since the likelihood is unchanged. Rather, it reflects how easy it is for the sampler to navigate the posterior surface when the data are centered.

A comparison of the trace plots reinforces the stability that centering the data provides. The traces. for the 0/1-coded data (on the left) are a bit more irregular, suggesting less efficient exploration of the posterior. In contrast, the centered parameterization produces tighter, more stable traces with less autocorrelation (on the right), indicating that the chains are mixing more effectively. This aligns with the much larger effective sample sizes observed for the centered model.

Finally, we compare the estimation of the log-odds ratios for the two models, just as we did before with the GLM models, and it is clear that the two Bayesian models also provide the same estimates of the contratsts:

get_lor_summary <- function(samp, model_name) {
  dt <- as.data.table(as_draws_df(samp))
  
  if (model_name == "0/1-coded") {
    dt[, lOR_A := beta_a]
    dt[, lOR_B := beta_b]
    dt[, lOR_AB := beta_a + beta_b + beta_ab]
  } else {
    dt[, lOR_A := gamma_a - 0.5 * gamma_ab]
    dt[, lOR_B := gamma_b - 0.5 * gamma_ab]
    dt[, lOR_AB := gamma_a + gamma_b]
  }
  
  dt[, .(
    mean_A = mean(lOR_A),
    mean_B = mean(lOR_B),
    mean_AB = mean(lOR_AB),
    sd_A = sd(lOR_A),
    sd_B = sd(lOR_B),
    sd_AB = sd(lOR_AB)
  )]
}

lor_01 <- get_lor_summary(samp_01, "0/1-coded")
lor_c  <- get_lor_summary(samp_c,  "centered")