1. Why unequal variances and correlations matter

The standard normal linear regression model

$y \mid \beta, \sigma, X \sim N(X\beta, \sigma^2 I)$

assumes:

  1. Linearity: $E(y \mid \theta, X)$ is a linear function of the predictors in $X$.
  2. Normality: Errors are normally distributed.
  3. Independence and equal variance: Errors are independent and all have the same variance $\sigma^2$.

In real data, none of these assumptions is exactly true. The key question is not “Are they exactly true?” but:

Does the mismatch between assumptions and reality damage the inferences we care about?

To answer this, one can:

  • Use posterior predictive checks to see how well the model reproduces important features of the data.
  • Fit alternative models (e.g., allow unequal variances or correlations) and see how sensitive the substantive conclusions are to these assumptions.

2. Diagnosing and handling nonlinearity

If the true relationship between $y$ and $X$ is nonlinear, the linear model’s mean structure $E(y\mid X, \theta) = X\beta$ is misspecified.

How to detect nonlinearity

  • Plot residuals vs explanatory variables.
  • If there’s a pattern (e.g., curved shape), the linear form may be inadequate.

How to fix it

  • Include transformed versions of predictors (e.g., squares, logs).
  • Example:
    • In medicine, improvement might depend on age in a nonlinear way: improvement increases for younger patients, then decreases for older ones.
    • Adding a quadratic term in age can capture this.

3. Diagnosing and handling nonnormality

Sometimes nonnormality is structurally obvious:

  • $y$ is discrete or bounded.
  • Residual plots (e.g., like those in the incumbency example) suggest heavy tails or skewness.

If nonnormality is serious:

  1. Apply a transformation to $y$ (e.g., log, logit).
  2. Or use a generalized linear model that matches the data type (binary, counts, etc.), instead of normal regression.

4. Diagnosing and handling unequal variances (heteroscedasticity)

Unequal error variances can be seen by:

  • Plotting absolute or squared residuals versus predictors.
  • If residual spread systematically varies with some covariate, variances are unequal.

Causes and remedies

  1. Missing covariates
    Example:
    • Agricultural yield regressed on soil and fertilizer may show unequal variances if precipitation is missing.
    • Adding precipitation as an explanatory variable can reduce heteroscedasticity.
  2. Structural reasons
    • When units are aggregates:
      • Example: hospital-level outcomes where each data point is an average of patients.
      • Variance is roughly proportional to $1/\text{sample size}$ in each hospital.
    • When measurement precision differs:
      • Example: two technicians with different skill levels produce measurements with different variances.

In such cases, unequal variances are intrinsic to the data generation process and must be explicitly modeled, not just “fixed” by adding more predictors.

5. Diagnosing and handling correlations in errors

If the errors $y_i – X_i\beta$ and $y_j – X_j\beta$ are correlated, independence is violated.

How to detect correlations

  • Look at correlation of residuals with respect to a relevant structure:
    • Time → autocorrelation plots if data are collected sequentially.
    • Cluster or family → residuals grouped by family, school, etc.

Why it matters

If correlation is ignored:

  • Posterior inferences about parameters look too precise because the data contain less information than nnn independent observations.
  • Predictions for new units are inaccurate if correlation structure is relevant.

Example: siblings’ heights

  • Even after controlling for age and sex, siblings’ heights remain correlated.
  • If one also controls for parents’ heights (two extra columns in $X$), the residual correlation among siblings becomes smaller, but usually still positive.
  • This shows:
    • Adding more predictors can reduce dependence, but some genuine correlation remains when there are unmodeled shared factors (e.g. genetics, shared environment).

Conclusion: sometimes adding covariates simplifies the covariance structure, but nonzero correlations usually remain when units are inherently related.

6. General model with unequal variances and correlations

To model both unequal variances and correlated errors, the model is generalized to:

$y \sim N(X\beta, \Sigma_y)$

where $\Sigma_y$ is an arbitrary symmetric positive definite $n \times n$ covariance matrix (not just $\sigma^2 I$).

This is more flexible, but:

  • Modeling and estimation become more complicated.
  • One must either specify $\Sigma_y$​ or put an informative prior on it.

7. Bayesian regression when the covariance matrix is known

First, consider the simple case where $\Sigma_y$​ is known.

Trick: linear transformation to standardize variance

Let $\Sigma_y^{1/2}$​ be a Cholesky factor of $\Sigma_y$​, so that

$\Sigma_y = \Sigma_y^{1/2} (\Sigma_y^{1/2})^T$

Multiply both sides of the regression equation by $\Sigma_y^{-1/2}$​:

$\Sigma_y^{-1/2} y \mid \theta, X \sim N(\Sigma_y^{-1/2} X\beta, I)$

Now we have a regression with identity covariance, so all the standard formulas apply if we replace:

  • $y$ with $y\tilde{y} = \Sigma_y^{-1/2} y$
  • $X$ with $X\tilde{X} = \Sigma_y^{-1/2} X$

Under a uniform prior for $\beta$, the posterior is:

$\hat{\beta} = (X^T \Sigma_y^{-1} X)^{-1} X^T \Sigma_y^{-1}$

$V_\beta = (X^T \Sigma_y^{-1} X)^{-1}$

and

$\beta \mid y \sim N(\hat{\beta}, V_\beta)$

(These are generalized least squares formulas.)

Prediction with known covariance

For new data $\tilde{y}$​ with predictors $\tilde{X}$, one must specify the joint covariance of old and new observations:

$\begin{pmatrix} y \\ \tilde{y} \end{pmatrix} \sim N \left( \begin{pmatrix} X\beta \\ \tilde{X}\beta \end{pmatrix}, \begin{pmatrix} \Sigma_y & \Sigma_{y,\tilde{y}} \\ \Sigma_{\tilde{y},y} & \Sigma_{\tilde{y}} \end{pmatrix} \right)$

Given $(y, \beta, \Sigma_y)$, the predictive mean and variance of $\tilde{y}$​ are:

  • Mean:
    • $E(\tilde{y} \mid y, \beta, \Sigma_y) = \tilde{X}\beta + \Sigma_{\tilde{y},y} \Sigma_y^{-1} (y – X\beta)$
  • Variance:
    • $\text{Var}(\tilde{y} \mid y, \beta, \Sigma_y) = \Sigma_{\tilde{y}} – \Sigma_{\tilde{y},y} \Sigma_y^{-1} \Sigma_{y,\tilde{y}}$

These follow from standard multivariate normal conditioning formulas.

8. Bayesian regression when $\Sigma_y$​ is unknown (general case)

Now suppose $\Sigma_y$​ is unknown.

Assume:

  • Prior on $\beta$ is uniform given $\Sigma_y$​:
    • $p(\beta \mid \Sigma_y) \propto 1$
  • There is some prior $p(\Sigma_y)$ on the covariance.

The joint posterior is:

$p(\beta, \Sigma_y \mid y) \propto p(\Sigma_y) \, N(y \mid X\beta, \Sigma_y)$

One can integrate out $\beta$ to get the marginal posterior of $\Sigma_y$​.

This leads to:

$p(\Sigma_y \mid y) \propto p(\Sigma_y)\, |\Sigma_y|^{-1/2} |V_\beta|^{1/2} \exp\left\{-\frac{1}{2}(y – X\hat{\beta})^T \Sigma_y^{-1} (y – X\hat{\beta})\right\}$

where $\hat{\beta}$​ and $V_\beta$​ depend on $\Sigma_y$​.

Problems:

  • This is easy to evaluate given $\Sigma_y$​, but sampling from it is difficult, especially since $\hat{\beta}$​ and $V_\beta$​ themselves depend on $\Sigma_y$​.
  • Specifying a reasonable prior on a full covariance matrix is also nontrivial.

So, instead of the most general case, it is common to work with restricted forms of $\Sigma_y$​ that are easier to model.

9. Special case: covariance known up to a scalar factor

Assume:

$\Sigma_y = Q_y \sigma^2$

where $Q_y$ is known and $\sigma^2$ is unknown.

Use a standard noninformative prior:

$p(\beta, \sigma^2) \propto \sigma^{-2}$

Define transformed data:

  • $y\tilde{y} = Q_y^{-1/2} y$
  • $X\tilde{X} = Q_y^{-1/2} X$

Then everything reduces to ordinary regression with variance $\sigma^2 I$, and:

$\hat{\beta} = (X^T Q_y^{-1} X)^{-1} X^T Q_y^{-1} y$

$V_\beta = (X^T Q_y^{-1} X)^{-1}$

and the usual scaled inverse-$\chi^2$ posterior for $\sigma^2$ with:

$s^2 = \frac{1}{n – k} (y – X\hat{\beta})^T Q_y^{-1} (y – X\hat{\beta})$

This is a direct generalization of the formula for ordinary regression (where $Q_y = I$).

10. Weighted linear regression

A very important special case is weighted regression, where:

  • $\Sigma_y$​ is diagonal:
    • $\Sigma_{ii} = \frac{\sigma^2}{w_i}$
  • Here $w_i$ are known weights (e.g., proportional to sample size or measurement precision).

Then:

$Q_y^{-1} = \text{diag}(w_1, \dots, w_n)$

Intuitively:

  • An observation with larger weight $w_i$​ has smaller variance and gets more influence in estimation.
  • This can be viewed as treating $w_i$ as an additional covariate that affects variance but not the mean.

11. Parametric models for unequal variances

More generally, each variance can depend on some weight and an unknown parameter:

$\Sigma_{ii} = \sigma^2 v(w_i, \phi)$

For example:

  • $v(w_i, \phi) = w_i^{-\phi}$
    • $\phi = 0$ → equal variances
    • $\phi = 1$ → variances proportional to $1 / w_i$​, like weighted regression
  • Or $v(w_i, \phi) = (1 – \phi) + \phi / w_i$ with similar behavior.

A simple noninformative prior is:

  • $\phi \sim \text{Uniform}(0,1)$

To avoid dependence on scale, the weights $w_i$are usually rescaled so that:

$\prod_{i=1}^n w_i = 1$

Given $\phi$, this again becomes weighted regression with

$Q_y = \text{diag}(v(w_i, \phi))$

Then one can derive a numerical posterior for $\phi$, often via evaluation on a grid:

$p(\phi \mid y) \propto p(\phi)\, |V_\beta|^{1/2} s^{-(n – k)}$

where $\hat{\beta}, V_\beta, s^2$ all depend on $\phi$.

Sampling proceeds via:

  1. Draw $\phi$ from its posterior (using grid and normalization).
  2. Given $\phi$, draw $\sigma^2$.
  3. Given $(\phi, \sigma^2)$, draw $\beta$.

12. Several unknown variance parameters (variance components)

Now consider the case where different groups of observations have different variances.

Suppose:

  • The $n$ observations are divided into $I$ groups (“batches”), with:
    • group $i$ containing $n_i$​ observations,
    • each group having its own variance $\sigma_i^2$​.
  • Then $\Sigma_y$​ is diagonal, with:
    • $\Sigma_{jj} = \sigma_i^2$​ for observation $j$ in group $i$.

Assume a noninformative prior:

$p(\beta, \sigma_1^2, \dots, \sigma_I^2) \propto \prod_{i=1}^I \sigma_i^{-2}$

This corresponds to treating each group’s variance like a separate experiment.

Posterior is:

$p(\beta, \sigma_1^2, \dots, \sigma_I^2 \mid y) \propto \left[\prod_{i=1}^I \sigma_i^{-(n_i + 2)}\right] \exp\left\{-\frac{1}{2} (y – X\beta)^T \Sigma_y^{-1}(y – X\beta)\right\}$

The conditional posterior of $\beta$ given $\Sigma_y$​ is again a generalized least squares normal distribution.

However, there is no simple closed form for the marginal posterior of the variance components, so iterative methods (e.g., EM, Gibbs, or Metropolis) are typically required.

Example: incumbency advantage with two variances

In the congressional incumbency example:

  • It is plausible that elections with incumbents running have different variability than open-seat elections.
  • Let:
    • $\sigma_1^2$​ = variance for elections with an incumbent running.
    • $\sigma_2^2$​ = variance for open-seat elections.

The model:

$y \sim N(X\beta, \Sigma_y)$

with:

  • $\Sigma_{ii} = \sigma_1^2$​ for incumbents,
  • $\Sigma_{ii} = \sigma_2^2$​ for open seats.

Fitting this model (using EM to find a mode and then Gibbs sampling) yields:

  • Estimates of β (incumbency effect) very similar to the equal-variance model.
  • Variance estimates differ:
    • Open-seat elections typically have larger variance.
    • This mainly affects the predictive distribution (uncertainty for different types of races), not the estimated average incumbency advantage.

Over time, the posterior medians of $\sigma_1$​ and $\sigma_2$​ can be plotted, showing differences between incumbents and open seats and year-to-year fluctuation. The pattern often looks noisy because some years have few open seats, suggesting a hierarchical model over time (smoothing variances across years) would be preferable.

13. General variance modeling structure

All the variance models discussed can be summarized as:

  • Mean model:
    • $E(y \mid X, \theta) = X\beta$
  • Variance model (on log scale):
    • $\log(\text{var}(y \mid X, \theta)) = W \phi$

where:

  • $W$ is a design matrix for variance (e.g., indicator variables for groups, functions of weights),
  • $\phi$ is a vector of variance parameters.

Given such a structure:

  • One can use iterative simulation methods (Gibbs, Metropolis, etc.) to sample from the joint posterior of $(\beta, \phi)$.
  • This yields flexible models that can handle heteroscedasticity and simple correlation structures.

One-sentence wrap-up

Unequal variances and correlated errors are common in practice; they can be diagnosed using residuals and often partially addressed by adding informative covariates, but when needed, Bayesian regression can explicitly model a general covariance structure $\Sigma_y$​, possibly with multiple variance components and parametric variance functions, combining generalized least squares formulas with prior distributions and iterative posterior simulation to obtain realistic inferences and predictions.