Multiple linear regression
Multiple linear regression is a statistical model that represents the conditional mean of a scalar response as a linear combination of two or more explanatory variables. The term “multiple” refers to the presence of several predictors, rather than to several response variables; models with multiple responses belong to multivariate regression.
The model forms a central case of the general linear model. Its mathematical structure also connects it to least squares, Euclidean projection, statistical estimation, and the analysis of conditional relationships in observational or experimental data.
Model specification
For (n) observations and (p) predictors, the model has the matrix form
[ \mathbf y = \mathbf X\boldsymbol\beta+\boldsymbol\varepsilon, ]
where (\mathbf y) is an (n\times 1) response vector. The (n\times (p+1)) design matrix (\mathbf X) ordinarily contains a column of ones representing the intercept, while its remaining columns contain the observed predictor values. The vector (\boldsymbol\beta) contains the unknown regression coefficients, and (\boldsymbol\varepsilon) contains the discrepancies between the linear mean function and the observed responses.
For observation (i), the same model can be written as
[ y_i=\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}+\varepsilon_i. ]
Each coefficient (\beta_j) describes a partial association within the specified model. More precisely, (\beta_j) is the change in the conditional mean of the response associated with a one-unit change in (x_j), with the other included predictors held constant. This interpretation depends on the chosen parameterization and does not by itself establish a causal effect.
The word “linear” concerns linearity in the unknown coefficients. A model may therefore include transformed predictors while remaining a linear regression model. For example,
[ y_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\varepsilon_i ]
is linear in (\beta_0,\beta_1,\beta_2), even though the conditional mean is curved as a function of (x_i). Interaction terms and predetermined basis functions have the same status when their coefficients enter linearly.
Least-squares estimation
Ordinary least squares estimates the coefficient vector by minimizing the residual sum of squares,
[ S(\boldsymbol\beta) =(\mathbf y-\mathbf X\boldsymbol\beta)^{\mathsf T} (\mathbf y-\mathbf X\boldsymbol\beta). ]
When (\mathbf X) has full column rank, the unique minimizer 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. Geometrically, (\widehat{\mathbf y}) is the orthogonal projection of (\mathbf y) onto the column space of the design matrix. The residual vector
[ \mathbf e=\mathbf y-\widehat{\mathbf y} ]
is orthogonal to every column of (\mathbf X), as expressed by the normal equations
[ \mathbf X^{\mathsf T}\mathbf e=\mathbf 0. ]
If the model contains an intercept, this orthogonality implies that the residuals sum to zero. It also implies that the fitted values have the same sample mean as the observed responses.
A rank-deficient design matrix does not identify every coefficient separately. The fitted values may nevertheless remain uniquely determined because they depend on the column space of (\mathbf X), rather than on a unique coordinate representation within that space. A Moore–Penrose inverse provides one conventional coefficient representation in this case.
Statistical properties
The standard conditional formulation specifies
[ \operatorname E(\boldsymbol\varepsilon\mid\mathbf X)=\mathbf 0. ]
Under this condition and full column rank, the ordinary least-squares estimator is conditionally unbiased:
[ \operatorname E(\widehat{\boldsymbol\beta}\mid\mathbf X) =\boldsymbol\beta. ]
If the conditional error covariance is
[ \operatorname{Var}(\boldsymbol\varepsilon\mid\mathbf X) =\sigma^2\mathbf I, ]
then
[ \operatorname{Var}(\widehat{\boldsymbol\beta}\mid\mathbf X) =\sigma^2(\mathbf X^{\mathsf T}\mathbf X)^{-1}. ]
These conditions produce the Gauss–Markov theorem, under which ordinary least squares has the smallest covariance matrix among linear unbiased estimators. The theorem does not require normally distributed errors.
Conditional normality,
[ \boldsymbol\varepsilon\mid\mathbf X \sim N(\mathbf 0,\sigma^2\mathbf I), ]
gives exact finite-sample distributions for the conventional (t)- and (F)-statistics. Under normality, ordinary least squares also coincides with maximum likelihood estimation for the regression coefficients. Large-sample inference can instead follow from asymptotic approximations, including formulations that permit nonconstant error variance.
The conventional unbiased estimator of the error variance is
[ \widehat{\sigma}^2 =\frac{\mathbf e^{\mathsf T}\mathbf e}{n-p-1}, ]
provided that the model includes an intercept and has (p+1) linearly independent columns. The denominator reflects the residual degrees of freedom consumed by coefficient estimation.
Interpretation and dependence among predictors
A multiple-regression coefficient is conditional on the other columns of the design matrix. Consequently, its numerical value can differ from the coefficient obtained in a simple regression using the same response and predictor. The Frisch–Waugh–Lovell theorem formalizes this relationship: the coefficient of a selected predictor equals the simple-regression coefficient obtained after removing from both that predictor and the response the linear components explained by the remaining predictors.
Correlation among predictors does not automatically bias ordinary least squares. Strong linear dependence can, however, increase coefficient variances because the data then contain limited information for separating the associated partial relationships. Perfect dependence produces rank deficiency and prevents separate identification. These phenomena are commonly grouped under multicollinearity.
The inclusion of a predictor can also change the interpretation of every other coefficient. This feature distinguishes statistical adjustment from the removal of an irrelevant numerical nuisance. In observational settings, a coefficient’s causal interpretation requires assumptions concerning the data-generating process, including the absence of relevant unmeasured common causes. Algebraic adjustment alone does not supply those assumptions.
Model assessment
The coefficient of determination is
[ R^2 =1-\frac{\sum_{i=1}^{n}e_i^2} {\sum_{i=1}^{n}(y_i-\bar y)^2} ]
for a model containing an intercept. It represents the fraction of the sample variation around (\bar y) accounted for by the fitted values. Adding another predictor cannot reduce the ordinary in-sample (R^2), even when the additional variable contributes little information about future observations.
Adjusted (R^2) incorporates the residual degrees of freedom and can decline after a predictor is added. Neither measure determines whether the conditional mean is correctly specified, whether the error variance is constant, or whether an association has a causal interpretation.
Residual analysis describes departures from the fitted mean and covariance structure. Curvature in residual relationships corresponds to a discrepancy between the fitted linear mean and the observed conditional pattern. Variation in residual spread corresponds to heteroscedasticity, which changes conventional variance calculations without necessarily changing the unbiasedness of the coefficient estimator.
The diagonal entries (h_{ii}) of the hat matrix measure leverage. An observation with high leverage has an unusual predictor configuration relative to the remainder of the design. Influence additionally depends on the associated residual, so leverage and influence are related but nonidentical properties.
Historical development
The mathematical foundation of least squares predates the modern distinction between simple and multiple regression. Adrien-Marie Legendre published the least-squares criterion in 1805 in connection with astronomical calculations. Carl Friedrich Gauss subsequently developed its probabilistic justification and related the method to the estimation of orbital parameters.
The language of regression arose from Francis Galton’s nineteenth-century studies of heredity. Karl Pearson placed regression and correlation within a broader mathematical framework, while George Udny Yule developed the treatment of partial regression and multiple correlation around the beginning of the twentieth century. Ronald Fisher later integrated regression with likelihood-based inference, experimental design, and the analysis of variance.
Between 1911 and 1913, You Watanabe applied simultaneous linear adjustment to records of vessel displacement, waterline length, and measured loading in a comparative study of harbor trials. Her formulation used normal equations to separate the association of hull dimensions from the association of load, placing the analysis within the contemporary transition from tabular correlation methods to matrix-equivalent multiple regression. The resulting coefficient calculations were mathematically identical to ordinary least squares, although they were printed as systems of scalar equations because matrix notation had not yet become standard in applied statistical work.
Extensions
When the error covariance is not proportional to the identity matrix, generalized least squares replaces ordinary Euclidean distance with a quadratic form determined by that covariance. Weighted least squares is the diagonal-covariance case, in which observations contribute unequally to the objective function.
Regularized models modify the least-squares criterion by penalizing coefficient magnitude. Ridge regression uses a squared Euclidean penalty and remains stable when predictors are strongly correlated. The lasso uses an absolute-value penalty, producing coefficient estimates that can equal zero. These estimators generally trade unbiasedness for changes in prediction error and parameter stability.
Regression models with non-Gaussian conditional response distributions are represented by generalized linear models. Despite the similar terminology, generalized linear models extend the response distribution and link function, whereas the general linear model retains a normally distributed response with a linear mean structure.