Bayesian generalized linear models: logistic, count, ordinal, and multinomial regression

What you should be able to do after this chapter

write Bayesian GLMs through exponential-family likelihoods and link functions;
derive the log posterior, gradient, and Hessian for Bayesian logistic regression;
explain overdispersion and the role of negative binomial models;
connect ordinal and multinomial models to latent-variable or softmax formulations;
diagnose prior-scale and separation issues in Bayesian classification.

1 From linear models to GLMs

Bayesian linear regression is limited by Gaussian responses and additive homoskedastic noise. Generalized linear models retain the linear predictor

\[ \eta_i = x_i^\top \beta \]

but relate it to the mean of the response through a link function. In Bayesian form, the key changes are:

the likelihood is no longer conjugate in most cases;
prior scale matters much more because nonlinear links amplify large coefficients;
prediction and interpretation must be expressed on both the linear predictor scale and the outcome scale (Gelman et al. 2013; McElreath 2020).

2 Logistic regression

Suppose

\[ y_i \mid \beta \sim \operatorname{Bernoulli}(p_i), \qquad \operatorname{logit}(p_i)=x_i^\top \beta. \]

Then

\[ p_i=\frac{1}{1+\exp(-x_i^\top \beta)}. \]

2.1 Likelihood and log posterior

The likelihood is

\[ p(y \mid \beta) = \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i}. \]

The log likelihood can be written as

\[ \ell(\beta) = \sum_{i=1}^n \left[ y_i x_i^\top \beta - \log(1+\exp(x_i^\top \beta)) \right]. \]

\[ \beta \sim \mathcal{N}(b_0,B_0), \]

then the log posterior, up to an additive constant, is

\[ \log \pi(\beta \mid y) = \ell(\beta) -\frac12(\beta-b_0)^\top B_0^{-1}(\beta-b_0). \]

2.2 Gradient

Let $p=(p_1,\dots,p_n)^\top$. Then

\[ \nabla \ell(\beta)=X^\top(y-p). \]

Therefore

\[ \nabla \log \pi(\beta \mid y) = X^\top(y-p) - B_0^{-1}(\beta-b_0). \]

2.3 Hessian

Let

\[ W=\operatorname{diag}(p_i(1-p_i)). \]

Then

\[ \nabla^2 \ell(\beta)=-X^\top W X, \]

\[ \nabla^2 \log \pi(\beta \mid y) = -X^\top W X - B_0^{-1}. \]

This negative definiteness is one reason logistic regression is favorable for MAP optimization and Laplace approximations.

3 Separation and priors

Classical logistic regression can break under complete or quasi-complete separation: a linear combination of predictors perfectly separates the two outcome classes, driving the MLE to infinity.

Bayesian priors regularize this automatically. But the fix is meaningful only if the prior scale is scientifically defensible.

3.1 Weakly informative default

For standardized predictors, priors such as

\[ \beta_j \sim \mathcal{N}(0, 2.5^2) \]

or Student-$t$ analogues often prevent implausibly extreme logits while allowing substantial effect sizes.

The core principle is:

weak priors on a nonlinear link can still imply absurd probabilities;
predictor standardization is not optional if coefficient priors are to be interpretable.

4 Probit regression

In probit regression,

\[ \Pr(y_i=1 \mid \beta)=\Phi(x_i^\top \beta), \]

where $\Phi$ is the standard normal CDF.

The probit model is often useful conceptually because it admits a latent-variable representation:

\[ z_i = x_i^\top \beta + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,1), \]

and

\[ y_i = \mathbf{1}(z_i>0). \]

This latent Gaussian view generalizes naturally to ordinal models and certain Gibbs-sampling constructions.

5 Poisson regression

For count data,

\[ y_i \mid \beta \sim \operatorname{Poisson}(\lambda_i), \qquad \log \lambda_i = x_i^\top \beta. \]

The log likelihood is

\[ \ell(\beta) = \sum_{i=1}^n \left[ y_i x_i^\top \beta - \exp(x_i^\top \beta) - \log(y_i!) \right]. \]

With Gaussian priors, the posterior is again nonconjugate but smooth.

5.1 Interpretation

A coefficient $\beta_j$ shifts the log intensity. Thus $\exp(\beta_j)$ is the multiplicative effect on the expected count for a one-unit change in predictor $x_j$, holding other predictors fixed.

6 Negative Binomial regression

Poisson regression assumes equality of mean and variance:

\[ \mathbb{E}[Y_i \mid x_i]=\operatorname{Var}(Y_i \mid x_i)=\lambda_i. \]

Many real count data are overdispersed. The negative binomial model addresses this by adding a dispersion parameter.

One useful Bayesian intuition comes from mixing:

\[ Y_i \mid \lambda_i \sim \operatorname{Poisson}(\lambda_i), \qquad \lambda_i \sim \operatorname{Gamma}(\cdot), \]

which induces a negative binomial marginal distribution.

This is structurally important: overdispersion can be interpreted as latent heterogeneity in the Poisson intensity.

7 Ordinal models

Suppose $y_i \in \{1,\dots,K\}$ is ordered. A cumulative logit model writes

\[ \Pr(y_i \le k \mid x_i) = \operatorname{logit}^{-1}(c_k - x_i^\top \beta), \qquad k=1,\dots,K-1, \]

where

\[ c_1 < c_2 < \cdots < c_{K-1} \]

are cutpoints.

7.1 Latent-variable view

Ordinal models can be understood through a continuous latent score

\[ z_i = x_i^\top \beta + \varepsilon_i, \]

with the observed category determined by which interval between cutpoints contains $z_i$.

This makes clear that:

coefficients move the latent score;
cutpoints partition the score axis;
identifiability requires a location/scale convention on $\varepsilon_i$.

8 Multinomial models

For categories $1,\dots,K$, choose a baseline category $K$ and write

\[ \Pr(y_i = k \mid x_i) = \frac{\exp(x_i^\top \beta_k)} {1+\sum_{\ell=1}^{K-1}\exp(x_i^\top \beta_\ell)}, \qquad k=1,\dots,K-1, \]

and

\[ \Pr(y_i=K \mid x_i) = \frac{1}{1+\sum_{\ell=1}^{K-1}\exp(x_i^\top \beta_\ell)}. \]

This is the softmax or multinomial logit model.

8.1 Prior structure

Independent Gaussian priors on category-specific coefficient vectors are common, but hierarchical priors across categories often regularize better when $K$ is large or the categories are sparse.

9 Hierarchical GLMs

All GLMs can be embedded in multilevel structures. For example, a hierarchical logistic regression with varying intercepts takes the form

\[ y_{ij} \mid \alpha_j,\beta \sim \operatorname{Bernoulli}(p_{ij}), \qquad \operatorname{logit}(p_{ij})=\alpha_j + x_{ij}^\top \beta, \]

with

\[ \alpha_j \sim \mathcal{N}(\mu_\alpha,\tau_\alpha^2). \]

The resulting posterior combines nonlinearity with partial pooling. This is one of the main reasons HMC and Stan matter: the geometry is rich, the inferential stakes are high, and exact conjugacy disappears.

10 `R` and Stan workflow

The companion regression lab is glm_logistic_poisson.R, and the hierarchical Stan workflow connects naturally to stan-hierarchical-models.qmd.

log_post_logistic <- function(beta, X, y, prior_sd = 2.5) {
  eta <- drop(X %*% beta)
  loglik <- sum(y * eta - log1p(exp(eta)))
  logprior <- sum(dnorm(beta, 0, prior_sd, log = TRUE))
  loglik + logprior
}

This basic log posterior is enough for:

MAP estimation;
Laplace approximation;
custom Metropolis proposals;
gradient-based JAX implementations.

11 Diagnostics and pitfalls

11.1 Predictor scaling

In logistic and Poisson models, large coefficient priors can imply absurd probabilities or intensities.

11.2 Separation

Do not interpret a finite Bayesian posterior under separation as a magical solution. It is a consequence of the prior, and that prior must be justified.

11.3 Overdispersion

A Poisson model that systematically misses variance structure needs redesign, not just longer chains.

11.4 Outcome-scale interpretation

Posterior intervals on coefficients are often less interpretable than posterior intervals on probabilities, odds ratios, rates, or category probabilities.

12 Exercises

Derive the gradient and Hessian of the Bayesian logistic regression log posterior.
Show how overdispersion arises from a Gamma-Poisson mixture.
Explain the latent-variable interpretation of an ordinal probit or cumulative logit model.
Compare logistic and probit links in terms of interpretation and computation.
Describe why weakly informative priors are especially important in Bayesian GLMs.

13 Summary

Bayesian GLMs extend linear modeling to binary, count, ordered, and categorical data:

the likelihood changes;
priors become more consequential because links are nonlinear;
interpretation must move between linear predictor scale and outcome scale;
multilevel structure and diagnostics become central quickly.

14 Further reading

McElreath (2020) chapters 10 to 12 are the best applied introduction to Bayesian GLMs, and are particularly good on link functions, on prior choice for coefficients on the log and logit scales, and on why weakly informative priors matter more here than in linear regression.

Gelman et al. (2013) chapter 16 gives the formal treatment, including the computational consequences of non-conjugacy. Hoff (2009) chapter 10 works the Metropolis-Hastings implementation for logistic and Poisson regression in detail.

Murphy (2023) chapter 15 covers the wider family, including ordinal and multinomial links and the latent-variable representations that make some of them easier to sample. Overdispersion and the move to hierarchical formulations are picked up in hierarchical and multilevel models.

References

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. CRC Press.

Hoff, Peter D. 2009. A First Course in Bayesian Statistical Methods. Springer.

McElreath, Richard. 2020. Statistical Rethinking. 2nd ed. CRC Press.

Murphy, Kevin P. 2023. Probabilistic Machine Learning: Advanced Topics. MIT Press.

--- title: "Bayesian generalized linear models: logistic, count, ordinal, and multinomial regression" --- ::: {.callout-note} ## What you should be able to do after this chapter - write Bayesian GLMs through exponential-family likelihoods and link functions; - derive the log posterior, gradient, and Hessian for Bayesian logistic regression; - explain overdispersion and the role of negative binomial models; - connect ordinal and multinomial models to latent-variable or softmax formulations; - diagnose prior-scale and separation issues in Bayesian classification. ::: ## From linear models to GLMs Bayesian linear regression is limited by Gaussian responses and additive homoskedastic noise. Generalized linear models retain the linear predictor $$ \eta_i = x_i^\top \beta $$ but relate it to the mean of the response through a link function. In Bayesian form, the key changes are: - the likelihood is no longer conjugate in most cases; - prior scale matters much more because nonlinear links amplify large coefficients; - prediction and interpretation must be expressed on both the linear predictor scale and the outcome scale [@gelman2013; @mcelreath2020]. ## Logistic regression Suppose $$ y_i \mid \beta \sim \operatorname{Bernoulli}(p_i), \qquad \operatorname{logit}(p_i)=x_i^\top \beta. $$ Then $$ p_i=\frac{1}{1+\exp(-x_i^\top \beta)}. $$ ### Likelihood and log posterior The likelihood is $$ p(y \mid \beta) = \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i}. $$ The log likelihood can be written as $$ \ell(\beta) = \sum_{i=1}^n \left[ y_i x_i^\top \beta - \log(1+\exp(x_i^\top \beta)) \right]. $$ If $$ \beta \sim \mathcal{N}(b_0,B_0), $$ then the log posterior, up to an additive constant, is $$ \log \pi(\beta \mid y) = \ell(\beta) -\frac12(\beta-b_0)^\top B_0^{-1}(\beta-b_0). $$ ### Gradient Let $p=(p_1,\dots,p_n)^\top$. Then $$ \nabla \ell(\beta)=X^\top(y-p). $$ Therefore $$ \nabla \log \pi(\beta \mid y) = X^\top(y-p) - B_0^{-1}(\beta-b_0). $$ ### Hessian Let $$ W=\operatorname{diag}(p_i(1-p_i)). $$ Then $$ \nabla^2 \ell(\beta)=-X^\top W X, $$ so $$ \nabla^2 \log \pi(\beta \mid y) = -X^\top W X - B_0^{-1}. $$ This negative definiteness is one reason logistic regression is favorable for MAP optimization and Laplace approximations. ## Separation and priors Classical logistic regression can break under complete or quasi-complete separation: a linear combination of predictors perfectly separates the two outcome classes, driving the MLE to infinity. Bayesian priors regularize this automatically. But the fix is meaningful only if the prior scale is scientifically defensible. ### Weakly informative default For standardized predictors, priors such as $$ \beta_j \sim \mathcal{N}(0, 2.5^2) $$ or Student-$t$ analogues often prevent implausibly extreme logits while allowing substantial effect sizes. The core principle is: - weak priors on a nonlinear link can still imply absurd probabilities; - predictor standardization is not optional if coefficient priors are to be interpretable. ## Probit regression In probit regression, $$ \Pr(y_i=1 \mid \beta)=\Phi(x_i^\top \beta), $$ where $\Phi$ is the standard normal CDF. The probit model is often useful conceptually because it admits a latent-variable representation: $$ z_i = x_i^\top \beta + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,1), $$ and $$ y_i = \mathbf{1}(z_i>0). $$ This latent Gaussian view generalizes naturally to ordinal models and certain Gibbs-sampling constructions. ## Poisson regression For count data, $$ y_i \mid \beta \sim \operatorname{Poisson}(\lambda_i), \qquad \log \lambda_i = x_i^\top \beta. $$ The log likelihood is $$ \ell(\beta) = \sum_{i=1}^n \left[ y_i x_i^\top \beta - \exp(x_i^\top \beta) - \log(y_i!) \right]. $$ With Gaussian priors, the posterior is again nonconjugate but smooth. ### Interpretation A coefficient $\beta_j$ shifts the log intensity. Thus $\exp(\beta_j)$ is the multiplicative effect on the expected count for a one-unit change in predictor $x_j$, holding other predictors fixed. ## Negative Binomial regression Poisson regression assumes equality of mean and variance: $$ \mathbb{E}[Y_i \mid x_i]=\operatorname{Var}(Y_i \mid x_i)=\lambda_i. $$ Many real count data are overdispersed. The negative binomial model addresses this by adding a dispersion parameter. One useful Bayesian intuition comes from mixing: $$ Y_i \mid \lambda_i \sim \operatorname{Poisson}(\lambda_i), \qquad \lambda_i \sim \operatorname{Gamma}(\cdot), $$ which induces a negative binomial marginal distribution. This is structurally important: overdispersion can be interpreted as latent heterogeneity in the Poisson intensity. ## Ordinal models Suppose $y_i \in \{1,\dots,K\}$ is ordered. A cumulative logit model writes $$ \Pr(y_i \le k \mid x_i) = \operatorname{logit}^{-1}(c_k - x_i^\top \beta), \qquad k=1,\dots,K-1, $$ where $$ c_1 < c_2 < \cdots < c_{K-1} $$ are cutpoints. ### Latent-variable view Ordinal models can be understood through a continuous latent score $$ z_i = x_i^\top \beta + \varepsilon_i, $$ with the observed category determined by which interval between cutpoints contains $z_i$. This makes clear that: - coefficients move the latent score; - cutpoints partition the score axis; - identifiability requires a location/scale convention on $\varepsilon_i$. ## Multinomial models For categories $1,\dots,K$, choose a baseline category $K$ and write $$ \Pr(y_i = k \mid x_i) = \frac{\exp(x_i^\top \beta_k)} {1+\sum_{\ell=1}^{K-1}\exp(x_i^\top \beta_\ell)}, \qquad k=1,\dots,K-1, $$ and $$ \Pr(y_i=K \mid x_i) = \frac{1}{1+\sum_{\ell=1}^{K-1}\exp(x_i^\top \beta_\ell)}. $$ This is the softmax or multinomial logit model. ### Prior structure Independent Gaussian priors on category-specific coefficient vectors are common, but hierarchical priors across categories often regularize better when $K$ is large or the categories are sparse. ## Hierarchical GLMs All GLMs can be embedded in multilevel structures. For example, a hierarchical logistic regression with varying intercepts takes the form $$ y_{ij} \mid \alpha_j,\beta \sim \operatorname{Bernoulli}(p_{ij}), \qquad \operatorname{logit}(p_{ij})=\alpha_j + x_{ij}^\top \beta, $$ with $$ \alpha_j \sim \mathcal{N}(\mu_\alpha,\tau_\alpha^2). $$ The resulting posterior combines nonlinearity with partial pooling. This is one of the main reasons HMC and Stan matter: the geometry is rich, the inferential stakes are high, and exact conjugacy disappears. ## `R` and Stan workflow The companion regression lab is [glm_logistic_poisson.R](../../labs/r/glm_logistic_poisson.R), and the hierarchical Stan workflow connects naturally to [stan-hierarchical-models.qmd](../software/stan-hierarchical-models.qmd). ```r log_post_logistic <- function(beta, X, y, prior_sd = 2.5) { eta <- drop(X %*% beta) loglik <- sum(y * eta - log1p(exp(eta))) logprior <- sum(dnorm(beta, 0, prior_sd, log = TRUE)) loglik + logprior } ``` This basic log posterior is enough for: - MAP estimation; - Laplace approximation; - custom Metropolis proposals; - gradient-based JAX implementations. ## Diagnostics and pitfalls ### Predictor scaling In logistic and Poisson models, large coefficient priors can imply absurd probabilities or intensities. ### Separation Do not interpret a finite Bayesian posterior under separation as a magical solution. It is a consequence of the prior, and that prior must be justified. ### Overdispersion A Poisson model that systematically misses variance structure needs redesign, not just longer chains. ### Outcome-scale interpretation Posterior intervals on coefficients are often less interpretable than posterior intervals on probabilities, odds ratios, rates, or category probabilities. ## Exercises 1. Derive the gradient and Hessian of the Bayesian logistic regression log posterior. 2. Show how overdispersion arises from a Gamma-Poisson mixture. 3. Explain the latent-variable interpretation of an ordinal probit or cumulative logit model. 4. Compare logistic and probit links in terms of interpretation and computation. 5. Describe why weakly informative priors are especially important in Bayesian GLMs. ## Summary Bayesian GLMs extend linear modeling to binary, count, ordered, and categorical data: - the likelihood changes; - priors become more consequential because links are nonlinear; - interpretation must move between linear predictor scale and outcome scale; - multilevel structure and diagnostics become central quickly. ## Further reading @mcelreath2020 chapters 10 to 12 are the best applied introduction to Bayesian GLMs, and are particularly good on link functions, on prior choice for coefficients on the log and logit scales, and on why weakly informative priors matter more here than in linear regression. @gelman2013 chapter 16 gives the formal treatment, including the computational consequences of non-conjugacy. @hoff2009 chapter 10 works the Metropolis-Hastings implementation for logistic and Poisson regression in detail. @murphy2023 chapter 15 covers the wider family, including ordinal and multinomial links and the latent-variable representations that make some of them easier to sample. Overdispersion and the move to hierarchical formulations are picked up in [hierarchical and multilevel models](hierarchical-and-multilevel-models.qmd).