Weakly informative priors for logistic regression

1. Motivation: why weakly informative priors in logistic regression?

In logistic regression, nonidentifiability is common. Two main problems:

1. Collinearity
   • Same issue as in linear regression: predictors are nearly linear combinations of each other.
   • This makes the information about some coefficients very weak.
2. Separation
   • A linear combination of predictors perfectly predicts the outcome.
   • For some direction in covariate space, all observations are $y = 0$; for another, all are $y = 1$.
   • When this happens, the maximum likelihood estimate for some coefficients diverges to $+\infty$ or $-\infty$.

Separation is especially common when:

• Predictors are binary (e.g., indicator variables).
• Data sets are not huge.

A common but bad “fix” is:

• Drop predictors until the model becomes identifiable.
• This often means dropping the strongest predictors, which is conceptually backward.

Bayesian logistic regression with weakly informative priors is proposed as a better alternative:

• Use priors that:
  • Are not highly informative,
  • But rule out absurdly huge effects,
  • And stabilize the computation even under separation.

The authors propose independent Student–$t$ priors for the coefficients, with carefully chosen degrees of freedom and scale.

---

2. Basic prior setup: independent Student–t priors

We consider a logistic regression:

• Model: $\Pr(y_i = 1) = \operatorname{logit}^{-1}(X_i \beta)$.
• Coefficients: $\beta = (\beta_1,\dots,\beta_J)$.

Prior assumptions:

• Coefficients are a priori independent:
• $$\beta_j \sim t_\nu(\text{mean } 0,\ \text{scale } s), \quad j = 1,\dots,J.$$
• Here:
  • $\nu$ is the degrees-of-freedom parameter.
  • $s$ is the scale.

Rationale:

• Choose $\nu$ and $s$ so that:
  • Priors are weak, not strongly informative,
  • But they keep coefficients in a reasonable range.
• The $t$ family has heavy tails, which:
  • Prevents over-shrinking truly large effects,
  • Gives robustness to outliers and rare configurations of predictors.

Computation:

• Can be done with full MCMC (e.g., Gibbs + Metropolis) for full Bayesian inference.
• Or approximated by:
  • An EM-style algorithm built on top of the classical IRLS (iteratively reweighted least squares) algorithm for logistic regression.
  • This gives approximate posterior modes and standard errors.

From the point-estimation perspective, the prior:

• Acts as a regularizer, shrinking unstable coefficients toward 0,
• Especially helpful for parameters that would otherwise be very poorly estimated or nonidentified.

---

3. Separation: the problem illustrated

When there is separation:

• The fitted logistic model can perfectly classify some subset of the data.
• Then the maximum likelihood estimate for some coefficients is:
  • Undefined (infinite), or
  • Numerically extreme and meaningless.

Even Bayesian inference with a uniform prior on $\beta$ fails:

• The posterior mode is also at $\pm\infty$ for the separated coefficient.

Example: vote from sex, ethnicity, and income

Model:

• Outcome: probability of Republican vote for president.
• Predictors: sex (female), ethnicity (black), and income (coded as $-2,-1,0,1,2$).
• Fit logistic regressions separately for multiple election years.

Observation:

• In the 1964 data, among 87 Black respondents, none supported the Republican candidate.
• This gives complete separation for the “black” indicator:
  • The maximum likelihood estimate for the “black” coefficient is $-\infty$.
  • Numerical algorithms (like glm in R) may return a large finite number and a huge standard error, but that is just an artifact of stopping after a finite number of iterations.

Profile likelihood for the “black” coefficient:

• If you plot the likelihood as a function of this coefficient (conditional on the other coefficients at their conditional MLEs):
  • It keeps increasing as $\beta_{\text{black}} \to -\infty$.
  • The “best” fit is at $-\infty$, which is clearly unrealistic in the substantive context.

A weakly informative prior:

• Pulls $\beta_{\text{black}}$ back to a large but finite negative value.
• This is more consistent with the idea that the probability of a Black voter choosing Republican is very low but not literally zero.

---

4. Adding a normal prior via augmented data

We now adapt the standard logistic regression algorithm to include normal priors as a first step.

Model:

• Likelihood: $\Pr(y_i = 1) = \operatorname{logit}^{-1}(X_i \beta)$.
• Classical IRLS:
  • At each iteration, approximate the log-likelihood with a quadratic in $\beta$.
  • This generates:
    • Pseudo-data $z_i$,
    • Pseudo-variances $(\sigma_i^z)^2$.
  • Then run weighted least squares to update $\beta$.
  • Iterate until convergence.

Now introduce a normal prior:

• Prior: $\beta_j \sim N(\mu_j,\sigma_j^2)$, for $j = 1,\dots,J$.

This can be folded into IRLS with pseudo-data augmentation:

1. Let:
   • $z$ be the vector of pseudo-data from the likelihood,
   • $\sigma^z$ be the vector of pseudo standard deviations,
   • $X$ be the design matrix for the likelihood.
2. Build augmented pseudo-data:
   • Response: $$z_* = \begin{pmatrix} z \\ \mu \end{pmatrix}.$$
   • Design matrix: $$X_* = \begin{pmatrix} X \\ I_J \end{pmatrix}.$$
   • Weights: $$w_* = \begin{pmatrix} 1 / (\sigma_1^z)^2 \\ \vdots \\ 1 / (\sigma_n^z)^2 \\ 1 / \sigma_1^2 \\ \vdots \\ 1 / \sigma_J^2 \end{pmatrix}.$$
3. Run weighted least squares on $(z_*, X_*, w_*)$ to update $\beta$.

Advantages:

• The augmented regression is always identified, because of the prior rows.
• So the updated estimate $\hat{\beta}$ is finite, even if the likelihood alone would be nonidentifiable due to separation or collinearity.

At convergence:

• $\hat{\beta}$ is the posterior mode under the normal prior.
• The covariance matrix $V_\beta$ is the inverse of the negative Hessian of the log posterior at $\hat{\beta}$.

---

5. Extending to Student–t priors via a normal–inverse–$\chi^2$ mixture

To implement Student–t priors, we use a scale-mixture-of-normals representation.

For each coefficient:

• Represent the $t$ prior as:
  • $\beta_j \mid \sigma_j^2 \sim N(\mu_j,\sigma_j^2)$,
  • $\sigma_j^2 \sim \text{Inv-}\chi^2(\nu_j,s_j^2)$.

Here:

• $\nu_j$ is the degrees of freedom.
• $s_j$ is the scale (roughly the “typical spread”).

This leads to an approximate EM algorithm:

Algorithm structure

1. Initialization
   • Set $\sigma_j^2 = s_j^2$ for each $j$.
   • Pick an initial $\beta$ (e.g., $\beta = 0$ or a rough estimate).
2. Iterative steps (repeat until convergence) Step 1 – IRLS update (approximate likelihood)
   • Given current $\beta$, compute:
     • Pseudo-data $z$,Pseudo-variances $(\sigma_i^z)^2$,
       via the usual logistic IRLS formulas.
     Construct $z_*, X_*, w_*$ from the likelihood and current prior scales $\sigma_j^2$.Run weighted least squares on $(z_*, X_*, w_*)$ to obtain:
     • Updated coefficient vector $\hat{\beta}$,Covariance matrix $V_\beta$ (from the regression).
   Step 2 – Approximate E-step
   • In the log posterior, each prior term involves $(\beta_j – \mu_j)^2 / \sigma_j^2$.Approximate the expectation of $(\beta_j – \mu_j)^2$ under the normal approximation: $$E[(\beta_j – \mu_j)^2 \mid \sigma, y] \approx (\hat{\beta}_j – \mu_j)^2 + (V_\beta)_{jj}.$$This is approximate because the likelihood is only approximated by a normal in $\beta$.
   Step 3 – M-step for $\sigma_j^2$
   • Maximize the approximate expected log posterior with respect to $\sigma_j^2$.The update has a closed form: $$\hat{\sigma}_j^2 = \frac{(\hat{\beta}_j – \mu_j)^2 + (V_\beta)_{jj} + \nu_j s_j^2}{1 + \nu_j}.$$This corresponds to the posterior mode of $\sigma_j^2$ when you have:
     • One “measurement” of squared deviation $(\hat{\beta}j – \mu_j)^2 + (V\beta)_{jj}$,And an inverse-$\chi^2$ prior with parameters $(\nu_j, s_j^2)$.
   Step 4 – Repeat
   • Use the updated $\hat{\sigma}_j^2$ in the prior rows of the next IRLS step.
   • Return to Step 1 and iterate.

After convergence:

• Use the final $\hat{\beta}$ and $V_\beta$ as the approximate posterior mode and covariance.
• The estimated $\hat{\sigma}_j^2$ reflect how much each coefficient was shrunk by the heavy-tailed prior.

---

6. Default weakly informative prior: scaling and Cauchy priors

The goal is a default prior that works reasonably across many applications without tuning.

Step 1: standardize inputs

Distinguish between:

• Inputs: raw variables (age, sex, income, etc.).
• Predictors: functions of inputs used in the model (intercept, age, sex, age × sex, etc.).

Standardization is done at the input level, not the predictor level.

Rules:

1. Binary inputs
   • Recode so that:
     • Mean is 0,
     • The two levels differ by 1.
   • Example: if 10% are black and 90% are not:
     • Define a centered variable $c.\text{black}$ with values $0.9$ (for black) and $-0.1$ (for non-black).
2. Other inputs (continuous)
   • Center so that mean is 0.
   • Rescale so that standard deviation is $0.5$.
   • Then a change of one “unit” in the standardized variable corresponds roughly to a realistic range, and is comparable to a symmetric binary variable with values $\pm 0.5$.

This scaling:

• Depends on the data, but that is acceptable.
• Can be viewed as a rough approximation to a more fully Bayesian treatment where the scale is itself a parameter.

Step 2: choosing $\nu$ and $s$ via a baseline example

Consider the simplest logistic model:

• A single parameter $\theta$ (intercept only).
• Likelihood corresponding to “half success, half failure”:
  • For example, one Bernoulli trial with equal weight on success and failure in a conceptual sense.
• The likelihood as a function of $\theta$ is proportional to $e^{\theta/2} / (1 + e^\theta)$.

This likelihood:

• Is roughly similar to a $t_7$ density with scale $2.5$.
• To be more conservative and allow heavier tails, the authors choose a Cauchy (i.e., $t_1$) distribution with scale $2.5$.

Step 3: default prior prescription

• For each non-intercept coefficient of a standardized input:
  • Use $\beta_j \sim \text{Cauchy}(0, 2.5)$.
• For the intercept:
  • Use $\beta_0 \sim \text{Cauchy}(0, 10)$.

Interpretation:

• A change of 5 on the logit scale is extremely large:
  • It moves probability from about $0.01$ to $0.5$, or from $0.5$ to $0.99$.
• The Cauchy$(0,2.5)$ prior on standardized inputs implies:
  • Most of the prior mass is on effects smaller than 5 in absolute value,
  • But there is non-negligible probability for larger values.
• For the intercept:
  • Cauchy$(0,10)$ says that for an “average” observation:
    • The success probability is probably between something like $10^{-9}$ and $1 – 10^{-9}$.
  • This is weaker (more diffuse) than the $2.5$-scale prior and is often more appropriate for the intercept.

---

7. Examples: bioassay and voting

Bioassay example

• Data: 20 animals, 4 dose levels, binary outcome (death or survival).
• This is a small dataset with no separation, but limited information.

Comparing fits:

• Maximum likelihood:
  • Estimated coefficient for standardized dose: about $10.2 \pm 6.4$.
• With the default priors (Cauchy(0,10) for intercept; Cauchy(0,2.5) for dose after rescaling $x$ to mean 0 and sd 0.5):
  • Estimate reduces to about $5.4 \pm 2.2$.

Interpretation:

• The prior strongly shrinks the effect, because the data are too weak to support an effect as large as $10.2$ on the logit scale.
• The Bayesian fit does not contradict the data:
  • Fit is still consistent with the observed counts within binomial variability for $n=5$ per group.
• The large shift in the estimate reflects:
  • How weak the data are,
  • How implausible extremely large logistic slopes are in typical applications.

Voting from ethnicity (the separation case, revisited)

• For the pre-election survey data:
  • Under flat priors, separation in 1964 drives $\beta_{\text{black}}$ to $-\infty$.
• Using the default Cauchy priors:
  • Coefficients in 1964 become finite and reasonable.
  • Coefficients in other years are almost unchanged relative to the flat-prior estimates.

If you compare Cauchy, $t_7$, and normal priors (all with center 0, scale 2.5 for non-intercepts, 10 for intercept):

• All three:
  • Stabilize the 1964 “black” coefficient,
  • Have little impact on years with rich information.
• Cauchy is slightly heavier-tailed and more conservative.

This shows how the default weakly informative prior works in routine practice:

• Does nothing dramatic when the data are strong.
• Prevents pathologies when the data are weak or separated.

---

8. Default prior vs real prior information

The default Cauchy$(0,2.5)$ prior:

• Does not capture all actual prior knowledge in the voting example.

In reality, we know:

• Black voters are historically less likely to vote Republican than white voters, so the “black” coefficient should be negative.
• Coefficients across years are related, so a hierarchical time-series model over years would be more appropriate.
• The 1964 Republican candidate had particularly conservative positions on civil rights, suggesting an especially negative “black” coefficient that year.

A richer prior could:

• Pool over years,
• Incorporate election-level predictors (e.g., candidate positions),
• Explicitly encode sign expectations (e.g., $\beta_{\text{black}} < 0$ with high probability).

But as a default:

• The Cauchy$(0,2.5)$ prior is simple, automatic, and stabilizing.
• It is weakly informative: it constrains extreme values but does not encode detailed substantive knowledge.
• Additional prior structure can be added if the problem demands it.