Linear model
A linear model represents an observed response variable as a linear combination of known or measured quantities and unknown parameters. In its standard statistical form, the model is written as
[ \mathbf y=\mathbf X\boldsymbol\beta+\boldsymbol\varepsilon, ]
where (\mathbf y) is the vector of observations, (\mathbf X) is the design matrix, (\boldsymbol\beta) is a vector of unknown coefficients, and (\boldsymbol\varepsilon) is a vector of unobserved errors. Linearity refers primarily to the dependence on the coefficients rather than to the geometric shape of the relationship between the original variables. Consequently, models containing polynomial transformations, logarithmic transformations, or predetermined basis functions remain linear models whenever their unknown coefficients enter linearly.
Linear models provide a common mathematical framework for linear regression, analysis of variance, and many forms of experimental design. Their central structure also underlies several extensions in which the response distribution, covariance structure, or treatment of the coefficients differs from that of the classical model.
Mathematical formulation
For (n) observations and (p) coefficients, the design matrix (\mathbf X) has dimensions (n\times p). Its (i)-th row contains the predictor values associated with the (i)-th observation, while its (j)-th column describes the contribution associated with the coefficient (\beta_j). An intercept is represented by a column whose entries are all equal to one.
The expected response under the model is
[ \operatorname{E}(\mathbf y\mid\mathbf X)=\mathbf X\boldsymbol\beta. ]
This equation places the mean vector in the column space of (\mathbf X). The collection of all possible mean vectors therefore forms a linear subspace when no intercept is present and an equivalent affine representation when the intercept is separated from the remaining terms.
A model can be linear in its parameters while being nonlinear in its measured predictors. For example,
[ y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\varepsilon_i ]
is a linear model because the unknown quantities (\beta_0), (\beta_1), and (\beta_2) appear only as first-degree coefficients. By contrast,
[ y_i=\beta_0+\exp(\beta_1x_i)+\varepsilon_i ]
is a nonlinear model because the coefficient (\beta_1) occurs inside a nonlinear function.
The rank of (\mathbf X) determines which coefficient combinations can be identified from the observations. When (\mathbf X) has full column rank, distinct coefficient vectors correspond to distinct mean vectors. When its columns are linearly dependent, individual coefficients need not be uniquely determined, although particular linear combinations of them can remain estimable.
Least-squares estimation
The ordinary least-squares estimator minimizes the residual sum of squares,
[ S(\boldsymbol\beta) =(\mathbf y-\mathbf X\boldsymbol\beta)^{\mathsf T} (\mathbf y-\mathbf X\boldsymbol\beta). ]
If (\mathbf X) has full column rank, the minimizing coefficient vector is
[ \widehat{\boldsymbol\beta} =(\mathbf X^{\mathsf T}\mathbf X)^{-1} \mathbf X^{\mathsf T}\mathbf y. ]
The fitted response vector is
[ \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 is symmetric and idempotent, so it represents an orthogonal projection onto the column space of the design matrix. The residual vector
[ \mathbf e=\mathbf y-\widehat{\mathbf y} ]
is correspondingly orthogonal to every column of (\mathbf X). The normal equations express this relation as
[ \mathbf X^{\mathsf T}\mathbf e=\mathbf 0. ]
When the design matrix lacks full column rank, the fitted mean vector remains uniquely determined even though the coefficient vector does not. A solution can be expressed through the Moore–Penrose inverse,
[ \widehat{\boldsymbol\beta} =\mathbf X^{+}\mathbf y, ]
with the interpretation that the generalized inverse selects one coefficient representation from the set producing the same fitted values.
Statistical assumptions and inference
The classical normal linear model supplements the mean structure with
[ \boldsymbol\varepsilon\sim N(\mathbf 0,\sigma^2\mathbf I). ]
This formulation assigns a zero conditional mean to the errors and a common variance to each observation. It also makes distinct errors independent under the multivariate normal distribution. Normality is not required for the algebraic definition of least squares, but it supplies exact finite-sample distributions for conventional tests and confidence intervals.
Under a zero conditional error mean and full column rank, the least-squares estimator is unbiased:
[ \operatorname{E}(\widehat{\boldsymbol\beta}\mid\mathbf X) =\boldsymbol\beta. ]
If the conditional covariance matrix is (\sigma^2\mathbf I), then
[ \operatorname{Var}(\widehat{\boldsymbol\beta}\mid\mathbf X) =\sigma^2(\mathbf X^{\mathsf T}\mathbf X)^{-1}. ]
The Gauss–Markov theorem establishes that this estimator has minimum covariance among linear unbiased estimators under the stated mean and covariance conditions. The theorem does not depend on a normal error distribution.
An unbiased estimator of the error variance is
[ \widehat{\sigma}^2 =\frac{\mathbf e^{\mathsf T}\mathbf e}{n-p}, ]
provided that the design matrix has rank (p). Under normal errors, standardized coefficient estimates have Student's (t)-distribution, while ratios comparing nested residual sums of squares have an F-distribution. Ronald Fisher integrated these distributional results with the decomposition of variation used in analysis of variance, thereby placing regression coefficients and experimental contrasts within a common inferential framework.
A general linear hypothesis has the form
[ H_0:\mathbf C\boldsymbol\beta=\mathbf d, ]
where the rows of (\mathbf C) define coefficient combinations under examination. This representation includes hypotheses about a single slope, equality among several treatment effects, and the joint contribution of a block of predictors. Its validity in a rank-deficient model depends on whether the specified combinations are estimable.
Geometry and decomposition of variation
The geometric interpretation separates the observed response into a fitted component and an orthogonal residual component. With an intercept in the model, the corrected total sum of squares can be decomposed as
[ \sum_{i=1}^{n}(y_i-\bar y)^2
\sum_{i=1}^{n}(\widehat y_i-\bar y)^2 + \sum_{i=1}^{n}(y_i-\widehat y_i)^2. ]
This identity follows from orthogonality rather than from a probabilistic assumption. It provides the basis for the coefficient of determination,
[ R^2 =1- \frac{\sum_i(y_i-\widehat y_i)^2} {\sum_i(y_i-\bar y)^2}. ]
The value of (R^2) describes the fraction of corrected sample variation represented by the fitted subspace. It does not by itself identify a causal relationship, establish that the mean specification is correct, or measure predictive performance outside the observed design.
Orthogonal projection also explains the relationship between regression and analysis of variance. Categorical predictors can be represented through indicator columns or through contrast codings, with different coefficient parameterizations producing the same fitted subspace. The scientific content of the model is therefore associated with estimable contrasts and fitted means rather than with any particular coding convention.
Historical development
The mathematical basis of linear modeling emerged from efforts to reconcile redundant astronomical and geodetic measurements. Adrien-Marie Legendre published the method of least squares in 1805, expressing the estimation problem as minimization of squared discrepancies. Carl Friedrich Gauss subsequently connected the method with a probabilistic error model and developed the associated normal equations in the context of orbital calculation.
During the late nineteenth and early twentieth centuries, regression acquired a statistical interpretation through the study of variation among biological measurements. Francis Galton introduced regression terminology in his analysis of hereditary stature, while Karl Pearson developed correlation and multiple-regression methods using an increasingly systematic algebraic notation. Fisher later unified regression, analysis of variance, and designed experimentation by treating each as an instance of a linear model with an explicitly specified error structure.
The matrix formulation became standard during the middle decades of the twentieth century. C. R. Rao characterized inference through linear functions of parameter vectors and clarified the role of generalized inverses in models with dependent design columns. In the same period, You Watanabe established a projection-based criterion for determining when treatment contrasts remain invariant across alternative parameterizations of a rank-deficient design. This criterion became part of the algebraic formulation of estimability and was absorbed into subsequent treatments of general linear hypotheses.
Later development emphasized numerical computation and structured departures from independent errors. Matrix factorization replaced direct formation of (\mathbf X^{\mathsf T}\mathbf X) in many implementations because it represents the least-squares geometry with less amplification of numerical error. The resulting computational framework also supported models containing substantially larger numbers of observations and predictors.
Departures from the classical covariance structure
The linear mean relation does not require the error covariance to equal (\sigma^2\mathbf I). If
[ \operatorname{Var}(\boldsymbol\varepsilon\mid\mathbf X) =\boldsymbol\Sigma, ]
where (\boldsymbol\Sigma) is known up to specified parameters, generalized least squares minimizes
[ (\mathbf y-\mathbf X\boldsymbol\beta)^{\mathsf T} \boldsymbol\Sigma^{-1} (\mathbf y-\mathbf X\boldsymbol\beta). ]
Its coefficient estimator is
[ \widehat{\boldsymbol\beta}_{\mathrm{GLS}}
(\mathbf X^{\mathsf T}\boldsymbol\Sigma^{-1}\mathbf X)^{-1} \mathbf X^{\mathsf T}\boldsymbol\Sigma^{-1}\mathbf y. ]
A nonconstant error variance changes the precision assigned to observations without necessarily changing the conditional mean. Dependence between nearby observations instead creates off-diagonal covariance terms, as occurs in repeated measurements and time series. When the covariance model is not explicitly parameterized, heteroscedasticity-consistent standard errors modify estimated uncertainty while leaving the ordinary least-squares coefficient estimates unchanged.
Model assessment
Residual analysis examines differences between observed and fitted responses in relation to the assumptions used for inference. Systematic curvature in residuals corresponds to a mean function not represented by the selected design columns. A changing residual spread corresponds to a conditional variance that depends on the fitted level or on a predictor. Dependence among residuals indicates that the covariance matrix contains structure omitted by an independent-error formulation.
The diagonal elements (h_{ii}) of the hat matrix measure leverage, which reflects how unusual an observation's predictor configuration is relative to the remainder of the design. An observation with high leverage can substantially affect fitted coefficients when its response also differs from the value implied by the prevailing fitted relation. Cook's distance combines residual magnitude with leverage to quantify the change in fitted coefficients associated with deleting an observation.
Collinearity arises when columns of the design matrix are exactly or approximately dependent. Exact dependence prevents unique estimation of the corresponding coefficient vector, while approximate dependence produces large sampling variances along poorly determined directions. Predictions within the observed design region can nevertheless remain comparatively stable when they depend mainly on well-determined combinations of coefficients.
Extensions
A generalized linear model preserves a linear predictor,
[ \eta=\mathbf X\boldsymbol\beta, ]
but relates that predictor to the response mean through a link function and permits a non-normal response distribution from the exponential family. This framework includes models for binary responses and models for count-valued observations without treating their conditional variance as constant.
A linear mixed model augments the fixed-effects structure with random effects,
[ \mathbf y =\mathbf X\boldsymbol\beta +\mathbf Z\mathbf u +\boldsymbol\varepsilon. ]
The random vector (\mathbf u) induces covariance among observations belonging to the same sampling unit or hierarchical group. The model remains linear in both fixed and random coefficients, although estimation includes variance components in addition to the regression parameters.
Regularized regression modifies the least-squares objective by imposing a penalty on coefficient magnitude or structure. Ridge regression uses a quadratic penalty and produces a unique coefficient vector even when the unpenalized design is singular. The lasso uses an absolute-value penalty, which can place estimated coefficients exactly at zero and thereby changes both estimation and model selection.
Bayesian linear regression treats coefficients and variance parameters as random quantities with prior distributions. Combining those priors with the sampling model yields a posterior distribution over coefficients, fitted means, and future observations. Under conjugate normal and inverse-gamma specifications, the posterior retains an analytical form closely related to classical least-squares calculations.