Bayesian Robust Regression With Heavy-Tailed Errors

Why Use Heavy-Tailed Regression Errors?

In Gaussian regression, an observation far from its fitted value can exert substantial influence on the coefficient estimates. That may be appropriate when unusually large errors are genuinely implausible. When extreme residuals are possible, however, a Gaussian model can make the analysis sensitive to a small number of observations. Bayesian regression with Student-t errors offers a model in which large errors are less surprising, while retaining a convenient conditional Gaussian structure.

The key construction assigns each observation a latent precision. Conditional on that precision, the error is Gaussian; integrating out the precision gives a Student-t error with heavier tails. The latent precisions also yield simple conditional distributions that support a Gibbs sampler.

Model: For observations \(i=1,\ldots,n\), let \(y_i\) be the outcome and let \(x_i\) be a \(p\)-vector of regressors. Write \(y_i=x_i^\mathsf{T}\beta+\epsilon_i\). Conditional on \(\lambda_i\) and \(\sigma^2\), assume \(\epsilon_i\mid\lambda_i,\sigma^2\sim N(0,\sigma^2/\lambda_i)\), with independent \(\lambda_i\sim\operatorname{Gamma}(\nu/2,\nu/2)\). Here the Gamma distribution is parameterized by shape and rate, and \(\nu>0\) is the degrees of freedom. The \(\lambda_i\) are independent across observations.

In this specification, \(\sigma^2\) is the scale squared of the marginal Student-t errors, not generally their variance. For \(\nu>2\), the marginal error variance is \(\nu\sigma^2/(\nu-2)\); for \(\nu\leq 2\), it is infinite. The degrees of freedom \(\nu\) control tail thickness: smaller values produce heavier tails. In what follows, \(\nu\) is treated as fixed.

From a Normal–Gamma Mixture to Student-t Errors

Let \(e_i=y_i-x_i^\mathsf{T}\beta\). Conditional on \(\lambda_i\), the error density is normal, and the prior density for \(\lambda_i\) is Gamma. Multiplying these densities gives a kernel in \(\lambda_i\) proportional to \(\lambda_i^{(\nu+1)/2-1}\exp\{-\lambda_i(\nu+e_i^2/\sigma^2)/2\}\). Integrating this Gamma-shaped kernel over positive \(\lambda_i\) yields the marginal density.

p(e_i\mid\sigma^2,\nu) = \frac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2)\sqrt{\nu\pi\sigma^2}} \left(1+\frac{e_i^2}{\nu\sigma^2}\right)^{-(\nu+1)/2}.

This is the density of a Student-t distribution with \(\nu\) degrees of freedom, location zero, and scale squared \(\sigma^2\). The mixture representation is more than a way to express the marginal density: it allows the model to adapt the conditional precision of each observation to its residual.

Worked Example: Deriving a Latent-Precision Update

Suppose \(\nu=4\), \(\sigma^2=4\), and an observation has residual \(e_i=8\). Under the model above, derive the conditional distribution of its latent precision.

The likelihood as a function of \(\lambda_i\) is proportional to \(\lambda_i^{1/2}\exp\{-\lambda_i e_i^2/(2\sigma^2)\}\). Multiplying by the Gamma prior kernel gives \(\lambda_i^{(\nu+1)/2-1}\exp\{-\lambda_i(\nu+e_i^2/\sigma^2)/2\}\). Therefore the conditional distribution is Gamma with shape \((\nu+1)/2\) and rate \((\nu+e_i^2/\sigma^2)/2\). Substituting the values,

\lambda_i\mid e_i,\sigma^2,\nu \sim \operatorname{Gamma}\left(\frac{5}{2},\frac{4+64/4}{2}\right) =\operatorname{Gamma}(2.5,10).

Its conditional mean is \(2.5/10=0.25\). Since the shape exceeds one, its expected reciprocal is rate divided by shape minus one: \(E(1/\lambda_i\mid e_i,\sigma^2,\nu)=10/(2.5-1)=6.667\). Thus the conditional variance evaluated at the mean precision is \(4/0.25=16\), while the conditional expected error variance is \(4E(1/\lambda_i\mid\cdot)=4\{10/(2.5-1)\}=80/3\), or approximately \(26.67\). These are different quantities: in general, \(E(1/\lambda_i)\) is not \(1/E(\lambda_i)\).

Conditional Distributions for Bayesian Regression

To specify a conjugate prior for the regression coefficients and error scale, use \(\beta\mid\sigma^2\sim N(b_0,\sigma^2 V_0)\) and \(\sigma^2\sim\operatorname{InverseGamma}(a_0,d_0)\), independently of the latent precisions. Here \(V_0\) is a positive-definite \(p\times p\) matrix, \(a_0,d_0>0\), and the inverse-Gamma density is proportional to \((\sigma^2)^{-a_0-1}\exp(-d_0/\sigma^2)\). Let \(X\) be the \(n\times p\) design matrix with rows \(x_i^\mathsf{T}\), \(y\) the outcome vector, and \(\Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n)\).

Conditional updates: Under this model and these priors, \[ \beta\mid\sigma^2,\lambda,y\sim N(m,\sigma^2 A^{-1}),\qquad A=V_0^{-1}+X^\mathsf{T}\Lambda X,\qquad m=A^{-1}(V_0^{-1}b_0+X^\mathsf{T}\Lambda y). \] For each observation, \[ \lambda_i\mid\beta,\sigma^2,y_i\sim \operatorname{Gamma}\left(\frac{\nu+1}{2},\frac{\nu+(y_i-x_i^\mathsf{T}\beta)^2/\sigma^2}{2}\right). \] Finally, \[ \sigma^2\mid\beta,\lambda,y\sim \operatorname{InverseGamma}\left( a_0+\frac{n+p}{2},\ d_0+\frac{1}{2}\left[(y-X\beta)^\mathsf{T}\Lambda(y-X\beta)+(\beta-b_0)^\mathsf{T}V_0^{-1}(\beta-b_0)\right]\right). \]

The coefficient update follows by combining the weighted Gaussian likelihood with the Gaussian prior. Completing the square in \(\beta\) gives precision matrix \(A/\sigma^2\) and mean \(m\). For the scale update, the likelihood contributes a power \((\sigma^2)^{-n/2}\), and the conditional prior for the \(p\)-dimensional coefficient vector contributes \((\sigma^2)^{-p/2}\). Their combined power gives the inverse-Gamma shape \(a_0+(n+p)/2\); the weighted residual sum of squares and prior quadratic form give its scale parameter.

The resulting Gibbs sampler cycles through these conditional distributions: sample \(\beta\) given the current scale and precisions, sample each \(\lambda_i\) given the current coefficients and scale, then sample \(\sigma^2\) given the current coefficients and precisions. These are conditional updates within the stated joint model; they do not mean that the coefficients can be estimated by treating one fixed set of latent precisions as known.

Worked Example: Calculating a Coefficient Update

Consider a regression through the origin with one coefficient, \(y_i=x_i\beta+\epsilon_i\). Suppose \(x=(1,2)^\mathsf{T}\), \(y=(1,2.4)^\mathsf{T}\), the current latent precisions are \((1,0.25)\), and the current \(\sigma^2\) is 1. Use the prior \(\beta\mid\sigma^2\sim N(0,4\sigma^2)\), so \(V_0=4\).

The posterior conditional precision, apart from the common factor \(1/\sigma^2\), is \(V_0^{-1}+\sum_i\lambda_i x_i^2\). Its value is \(0.25+1(1^2)+0.25(2^2)=2.25\). The corresponding numerator for the conditional mean is \(V_0^{-1}b_0+\sum_i\lambda_i x_i y_i=0+1(1)(1)+0.25(2)(2.4)=2.2\). Thus

\beta\mid\sigma^2,\lambda,y \sim N\left(\frac{2.2}{2.25},\frac{1}{2.25}\right) =N(0.9778,0.4444).

The second observation contributes less to the data precision because its current \(\lambda_i\) is 0.25. The update still incorporates the prior and uncertainty about \(\beta\); it is not simply an unweighted regression fit.

How Heavy Tails Downweight Outlying Observations

The latent-precision update makes the robust mechanism explicit. A large squared residual increases the conditional rate of \(\lambda_i\), reducing its conditional mean. The corresponding observation then contributes less to \(X^\mathsf{T}\Lambda X\) and \(X^\mathsf{T}\Lambda y\) in the coefficient update. Small residuals tend to be associated with precisions closer to one, the prior mean of \(\lambda_i\).

This is an adaptive weighting mechanism, not a rule that permanently removes observations. Each \(\lambda_i\) is uncertain and is updated along with the coefficients and scale. Posterior coefficient summaries average over that uncertainty, so the final Bayesian estimates are not generally identical to a single weighted least-squares fit.

Worked Example: A Large Residual and Its Precision

Use \(\nu=4\) and \(\sigma^2=1\). Compare an observation with residual \(e_i=0.5\) to one with residual \(e_i=6\). Their conditional precision distributions have the same shape, \((4+1)/2=2.5\), but different rates.

For the smaller residual, the rate is \((4+0.25)/2=2.125\), so \(E(\lambda_i\mid e_i)=2.5/2.125=1.176\). For the larger residual, the rate is \((4+36)/2=20\), so \(E(\lambda_i\mid e_i)=2.5/20=0.1250\). The large-residual observation therefore has much lower expected conditional precision.

The corresponding expected conditional variances are not found by dividing by these means. Since \(E(1/\lambda_i)=\text{rate}/(\text{shape}-1)\), they are \(1\{2.125/(2.5-1)\}=1.417\) and \(1\{20/(2.5-1)\}=13.33\), respectively. This calculation illustrates how a large residual is compatible with greater conditional error variance in the mixture representation.

Worked Example: Comparing a Gaussian Fit With a Downweighted Update

Suppose a hypothetical data set has three points \((x_i,y_i)=(0,0),(1,1),(2,10)\), and the regression includes an intercept and slope. The third point has a much larger outcome than the first two. For comparison, first calculate the ordinary least-squares line, then calculate a weighted least-squares line using current precision values \((1,1,0.05)\). The latter is a conditional likelihood calculation for those fixed precisions, without a coefficient prior.

For ordinary least squares, \(\bar{x}=1\) and \(\bar{y}=11/3\). The slope is \[ \frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})}{\sum_i(x_i-\bar{x})^2} =\frac{10}{2}=5, \] and the intercept is \(\bar{y}-5\bar{x}=-4/3\). The fitted line is therefore approximately \(y=-1.333+5x\).

For the weighted calculation, \(\sum\lambda_i=2.05\), \(\sum\lambda_i x_i=1.1\), \(\sum\lambda_i y_i=1.5\), \(\sum\lambda_i x_i^2=1.2\), and \(\sum\lambda_i x_i y_i=2.0\). The weighted slope is \[ \frac{(2.05)(2.0)-(1.1)(1.5)}{(2.05)(1.2)-(1.1)^2} =\frac{2.45}{1.25}=1.96. \] The intercept is \(\{1.5-(1.96)(1.1)\}/2.05=-0.3200\). The weighted line, \(y=-0.3200+1.96x\), is much closer to the first two points because the third point has lower precision. In a full Bayesian analysis, posterior uncertainty in the precisions and the coefficient prior also affect the result; this fixed-weight comparison isolates the role of downweighting.

Common Pitfalls and Interpretation

Confusing scale and variance. The parameter \(\sigma^2\) is the Student-t scale squared. Only when \(\nu>2\) does the marginal variance exist, and then it is \(\nu\sigma^2/(\nu-2)\).
Taking the reciprocal of an expected precision. In general, \(E(1/\lambda_i)\ne 1/E(\lambda_i)\). For a Gamma variable with shape \(a>1\) and rate \(b\), \(E(1/\lambda_i)=b/(a-1)\), whereas \(1/E(\lambda_i)=b/a\).
Calling a conditional weighted fit the posterior estimate. The weights are latent, the coefficients have a prior, and Bayesian summaries account for posterior uncertainty across all unknowns.
Treating robustness as automatic protection from every modeling problem. Student-t errors address heavy-tailed outcome errors under the specified regression model. They do not correct omitted-variable bias, an unsuitable mean function, or errors in the recorded data.

Key takeaway: Student-t regression represents each error as conditionally Gaussian with a latent precision. Large residuals tend to produce smaller precisions, reducing their contribution to coefficient updates, while posterior inference still propagates uncertainty about the precisions, coefficients, and error scale.

Check Your Understanding

Work through these questions using the model and parameterizations defined above.

For a Student-t error with \(\nu=5\) and \(\sigma^2=3\), what is its marginal variance?
With \(\nu=6\), \(\sigma^2=2\), and residual \(e_i=4\), derive the shape and rate of the conditional Gamma distribution for \(\lambda_i\).
Why does \(E(1/\lambda_i)\) generally differ from \(1/E(\lambda_i)\), and what condition is needed for the former to be finite under a Gamma distribution?
In the coefficient update, where do the latent precisions appear, and how does a small precision affect an observation’s contribution?
Why is a weighted least-squares line using one set of latent precisions not, by itself, a complete Bayesian posterior summary?