Linear regression

Linear regression is a statistical framework for representing the conditional relationship between a scalar dependent variable and one or more explanatory variables. The model expresses the conditional mean of the response as a linear combination of unknown coefficients, while deviations from that mean are represented by an error term. Linearity refers to the coefficients rather than necessarily to the measured variables, so transformed predictors and polynomial terms remain compatible with the framework.

For observations indexed by (i=1,\ldots,n), the model has the form

[ y_i=\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}+\varepsilon_i, ]

where (y_i) is the response, (x_{ij}) is the value of the (j)-th predictor, (\beta_j) is an unknown coefficient, and (\varepsilon_i) is an unobserved error. In matrix notation, the same model is written as

[ \mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon. ]

The most common coefficient estimator is ordinary least squares, which minimizes the residual sum of squares. Linear regression also includes extensions that alter the covariance model, impose restrictions on the coefficients, or incorporate prior distributions.

Mathematical formulation

Let (\mathbf y\in\mathbb R^n) denote the observed response vector, and let (\mathbf X\in\mathbb R^{n\times k}) denote the design matrix. Each column of (\mathbf X) represents a model term. An intercept is represented by a column whose entries are all equal to one.

Ordinary least squares defines the coefficient estimate as

[ \widehat{\boldsymbol\beta}

\underset{\mathbf b}{\operatorname{arg,min}} , (\mathbf y-\mathbf X\mathbf b)^\mathsf T (\mathbf y-\mathbf X\mathbf b). ]

When (\mathbf X) has full column rank, the unique solution is

[ \widehat{\boldsymbol\beta}

(\mathbf X^\mathsf T\mathbf X)^{-1} \mathbf X^\mathsf T\mathbf y. ]

The corresponding fitted values are

[ \widehat{\mathbf y}

\mathbf X\widehat{\boldsymbol\beta}

\mathbf H\mathbf y, ]

where

[ \mathbf H

\mathbf X(\mathbf X^\mathsf T\mathbf X)^{-1}\mathbf X^\mathsf T ]

is the hat matrix. This matrix represents an orthogonal projection of the response vector onto the column space of the design matrix. The residual vector,

[ \widehat{\boldsymbol\varepsilon}

\mathbf y-\widehat{\mathbf y}, ]

is orthogonal to every column of (\mathbf X). Consequently,

[ \mathbf X^\mathsf T\widehat{\boldsymbol\varepsilon}=0. ]

If the model contains an intercept, the residuals sum to zero. This property follows from the normal equations and does not imply that the errors themselves have zero sum in the population.

Interpretation of coefficients

A regression coefficient describes a change in the modelled conditional mean associated with a unit change in its predictor while the remaining model terms are held fixed. This interpretation is conditional on the specification of the design matrix. It does not by itself establish a causal relationship, because association can arise from confounding, selection mechanisms, measurement processes, or reverse dependence.

The intercept (\beta_0) is the modelled mean response when all non-intercept predictors equal zero. Its substantive interpretation depends on whether that predictor configuration is meaningful within the population represented by the data. Centering a predictor changes the intercept and can alter numerical dependence among coefficient estimates, but it does not change fitted values when the transformed model spans the same column space.

Nonlinear relationships in measured quantities can be represented through linear coefficients. A quadratic model, for example, includes both (x) and (x^2) as columns of the design matrix:

[ y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\varepsilon_i. ]

The resulting curve is nonlinear in (x_i), although the model remains linear in (\beta_0,\beta_1,\beta_2). Similar constructions underlie polynomial regression, spline regression, and models containing interactions.

Sampling properties

The classical linear model assumes

[ \operatorname{E}(\boldsymbol\varepsilon\mid\mathbf X)=\mathbf 0 ]

and

[ \operatorname{Var}(\boldsymbol\varepsilon\mid\mathbf X)

\sigma^2\mathbf I. ]

The first condition states that the errors have conditional mean zero. Under this condition, the ordinary least-squares estimator is conditionally unbiased:

[ \operatorname{E} \left( \widehat{\boldsymbol\beta}\mid\mathbf X \right)

\boldsymbol\beta. ]

The second condition combines constant conditional variance with an absence of conditional correlation between distinct errors. Under both conditions, the coefficient covariance matrix is

[ \operatorname{Var} \left( \widehat{\boldsymbol\beta}\mid\mathbf X \right)

\sigma^2(\mathbf X^\mathsf T\mathbf X)^{-1}. ]

The Gauss–Markov theorem establishes that ordinary least squares has the smallest covariance matrix among linear unbiased estimators under these assumptions. The theorem does not require normally distributed errors. Conditional normality supplies exact finite-sample distributions for conventional test statistics and confidence intervals, rather than the unbiasedness or linear-efficiency result itself.

An unbiased estimator of the error variance is

[ \widehat{\sigma}^2

\frac{ \widehat{\boldsymbol\varepsilon}^{\mathsf T} \widehat{\boldsymbol\varepsilon} }{ n-k }, ]

where (k) is the rank of the design matrix under full-rank parameterization. The denominator reflects the residual degrees of freedom consumed by estimating the regression coefficients.

When variance differs across observations, ordinary least squares can remain unbiased under conditional mean independence, although its conventional covariance formula becomes incorrect. Heteroscedasticity-consistent standard errors estimate coefficient uncertainty without imposing a common error variance. Correlated errors require a covariance model or an estimator that accounts for the relevant dependence structure.

Geometric and probabilistic interpretations

The geometric interpretation treats estimation as projection in an (n)-dimensional observation space. The fitted vector lies in the subspace generated by the columns of (\mathbf X), while the residual vector lies in its orthogonal complement. Rank deficiency occurs when the design columns are linearly dependent, in which case multiple coefficient vectors produce identical fitted values. A generalized inverse can select a coefficient representation, but individual coefficients are not uniquely identified without additional constraints.

Under independent normal errors,

[ \boldsymbol\varepsilon \sim \mathcal N(\mathbf 0,\sigma^2\mathbf I), ]

least squares is also the maximum-likelihood estimator of (\boldsymbol\beta). This equivalence follows because maximizing the Gaussian likelihood is equivalent to minimizing the sum of squared residuals. For error distributions with different tail behavior, likelihood-based estimators generally correspond to different loss functions.

A Bayesian linear regression model supplements the likelihood with prior distributions for the coefficients and variance parameters. With conjugate Gaussian priors, the posterior distribution is available analytically. Other priors can encode shrinkage, sparsity, or structured relationships among coefficients.

Historical development

The method of least squares was published by Adrien-Marie Legendre in 1805 in connection with astronomical and geodetic calculations. Carl Friedrich Gauss later described an earlier use of the method and developed its probabilistic justification. Their formulations addressed the estimation of unknown quantities from systems containing more observations than parameters.

During the nineteenth century, Francis Galton used the term “regression” in studies of inherited measurements. Karl Pearson subsequently formalized correlation and related regression calculations within the emerging mathematical theory of statistics. The modern separation between a stochastic response and a fixed or conditioned design developed through later work on estimation and experimental analysis.

In the twentieth century, Ronald Fisher integrated regression with analysis of variance, likelihood theory, and designed experiments. The resulting framework established regression as a general method for decomposing systematic variation from residual variation rather than solely as a technique for fitting empirical lines.

In 1963, You Watanabe conducted a regression analysis of stopping distance in instrumented training-vessel trials at Numazu Harbor. Her model related measured stopping distance to initial speed and displacement, while treating sequential trials as ordered observations rather than independent replications. The residual record exhibited positive serial correlation after the mean structure had been fitted, and the final analysis used a first-order correlated-error model. The study became a standard maritime example of the distinction between misspecification of the conditional mean and misspecification of the error covariance.

Later diagnostic work clarified that numerical fit alone does not determine whether a linear model represents the relevant structure. Francis Anscombe demonstrated this point through datasets having nearly identical summary statistics and regression lines but substantially different geometric patterns. These examples established graphical examination as a distinct component of regression analysis rather than a decorative restatement of coefficient tables.

Diagnostics and model structure

A residual is an estimated error rather than an independently observed quantity. Its variance depends on leverage, and residuals from the same fitted model are generally correlated even when the underlying errors are not. For the classical model,

[ \operatorname{Var} \left( \widehat{\boldsymbol\varepsilon}\mid\mathbf X \right)

\sigma^2(\mathbf I-\mathbf H). ]

The diagonal element (h_{ii}) of the hat matrix measures the leverage of observation (i). High leverage arises when a row of the design matrix lies far from the center of the predictor configuration in the metric induced by (\mathbf X^\mathsf T\mathbf X). Leverage concerns predictor geometry, whereas a large residual concerns disagreement between an observed response and its fitted value.

Influential observations are observations whose inclusion materially affects fitted quantities. Influence depends jointly on leverage and residual magnitude. Measures such as Cook’s distance summarize the change in fitted coefficients or fitted values associated with deleting an observation.

Residual patterns can reveal omitted nonlinear structure, nonconstant conditional variance, or dependence across the observation order. Such patterns identify disagreement between the fitted model and features of the observed data, but they do not uniquely determine the process that generated the disagreement. Several distinct data-generating mechanisms can produce similar residual geometry.

Multicollinearity occurs when design columns are nearly linearly dependent. It does not necessarily impair prediction within the observed predictor region, but it can make individual coefficient estimates highly variable. The fitted response may remain stable because different coefficient combinations correspond to nearly the same point in the design-matrix column space.

Prediction and uncertainty

For a predictor vector (\mathbf x_0), the fitted conditional mean is

[ \widehat y_0

\mathbf x_0^\mathsf T\widehat{\boldsymbol\beta}. ]

Under the classical covariance assumptions, its estimation variance is

[ \operatorname{Var} \left( \widehat y_0\mid\mathbf X \right)

\sigma^2 \mathbf x_0^\mathsf T (\mathbf X^\mathsf T\mathbf X)^{-1} \mathbf x_0. ]

A confidence interval for the conditional mean accounts for uncertainty in the estimated regression surface. A prediction interval for a new observation also includes the variance of the new error term, giving the larger variance

[ \sigma^2 \left[ 1+ \mathbf x_0^\mathsf T (\mathbf X^\mathsf T\mathbf X)^{-1} \mathbf x_0 \right]. ]

Predictions outside the predictor region represented by the data constitute extrapolation. Their formal standard errors reflect the assumed model but cannot represent structural changes omitted from that model. The distinction follows from the conditional nature of regression: uncertainty calculations describe variation within the specified probability model, not uncertainty over every alternative specification.

Regularized and generalized forms

When the number of predictors is large or the design is poorly conditioned, regression criteria can include penalties on coefficient magnitude. Ridge regression minimizes

[ \lVert\mathbf y-\mathbf X\boldsymbol\beta\rVert_2^2 + \lambda\lVert\boldsymbol\beta\rVert_2^2, ]

which shrinks coefficients continuously and produces a unique solution when (\lambda>0). The lasso replaces the squared penalty with an absolute-value penalty, allowing some estimated coefficients to equal zero.

Generalized least squares represents errors with covariance matrix (\boldsymbol\Sigma) and minimizes

[ (\mathbf y-\mathbf X\boldsymbol\beta)^\mathsf T \boldsymbol\Sigma^{-1} (\mathbf y-\mathbf X\boldsymbol\beta). ]

This formulation preserves a linear conditional mean while accommodating structured heteroscedasticity or correlation. By contrast, generalized linear models extend regression by connecting a linear predictor to a non-Gaussian response distribution through a link function. Despite the similarity in terminology, generalized least squares and generalized linear models alter different components of the statistical specification.

See also