Multiple Linear Regression

A multiple linear regression model relates one response variable to several predictor variables through a function that is linear in its unknown coefficients. NIST describes this class of models as linear least squares models: each term multiplies a known predictor by an unknown parameter, at most one parameter has no matching predictor (the intercept), and the terms are added.

For observation $i$,

$y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_p{x}_{ip}+{\epsilon }_i$.

$y_i$ is the response. $x_{ij}$ is the value of predictor $j$. ${\beta }_0$ is the intercept. ${\beta }_j$ are the remaining regression coefficients. ${\epsilon }_i$ is the error term in the model, not an observed residual.

Linear does not mean a straight line in the predictors. The model is linear in the unknown parameters. NIST gives the quadratic example

$y={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+\epsilon$

which is linear in ${\beta }_0,{\beta }_1,{\beta }_2$ even though it is nonlinear in $x$. A model such as $y={\beta }_0+{\beta }_0{\beta }_{1}x+\epsilon$ is linear in $x$ but not linear in the parameters, so it is not a linear least squares model in this sense.

With the intercept stored in the first column of the design matrix, the same model is

$y=X\beta +\epsilon$

with $y\in {\mathbb{R}}^{n\times 1}$, $X\in {\mathbb{R}}^{n\times (p+1)}$, $\beta \in {\mathbb{R}}^{(p+1)\times 1}$, and $\epsilon \in {\mathbb{R}}^{n\times 1}$. This page keeps that intercept convention throughout.

Linear least squares is the standard starting point for relating a quantitative response to several predictors when the scientific model is linear in its parameters. NIST notes that many processes are well described by such models, either because they are approximately linear or because a linear approximation is useful over a limited range, and that the associated theory supports interval estimates when the usual process-modeling assumptions hold.

Direct multiple linear regression becomes difficult in common spectral calibration settings. Wavelengths are often strongly correlated. The number of variables can approach or exceed the number of calibration samples, so $X$ is poorly conditioned or rank deficient. Inverse least squares calibration, a multiple linear regression of a response on spectral predictors, is one of the conventional methods against which principal component regression and partial least squares are compared (Haaland and Thomas 1988).

Those latent variable methods are motivated by collinear and high dimensional predictors. This page does not claim that PCR or PLS is universally better than MLR.

Multiple linear regression specifies the model. Ordinary least squares specifies an estimation criterion for the coefficients. They are related and they are not the same concept.

OLS chooses $\hat{\beta }$ to minimize the sum of squared residuals. NIST states the least squares criterion as the sum of squared deviations between the observed responses and the fitted model values. In matrix form that criterion is

$\hat{\beta }={\mathrm{argmin}}_{\beta }\Vert y-X\beta {\Vert }_2^2$.

Equivalently,

$\hat{\beta }={\mathrm{argmin}}_{\beta }(y-X\beta {)}^{\mathsf{T}}(y-X\beta )$.

The fitted responses are $\hat{y}=X\hat{\beta }$. The residual vector is $e=y-\hat{y}$, with $e_i=y_i-{\hat{y}}_i$. Residuals are observed minus fitted responses. They are not automatically measurement error.

The residual sum of squares is $\mathrm{RSS}=e^{\mathsf{T}}e={\sum }_i(y_i-{\hat{y}}_i{)}^2$. OLS is the estimator that minimizes RSS for a given design matrix and response.

For a new observation the predictors must use the same column order and the same intercept convention,

${\hat{y}}_{\mathrm{new}}=X_{\mathrm{new}}\hat{\beta }$.

Do not refit the model on the new rows. If a preprocessing transformation was estimated from training data, apply those training parameters to the new predictors. Predictor scaling is a modeling choice, not a requirement of OLS. See Normalization & Scaling.

The intercept is included as the first column of $X$. If an intercept is omitted, the column dimension of $X$ is $p$ rather than $p+1$, and the remaining equations keep the same form.

Computing $\hat{\beta }$ requires solving the least squares problem. It does not require normally distributed predictors, and it does not require normally distributed errors. Exact multicollinearity, meaning linear dependence among columns of $X$, removes uniqueness of the closed form inverse solution. Strong but imperfect collinearity can make individual coefficient estimates unstable without making prediction from the fitted values impossible. Correlated predictors do not automatically invalidate every MLR model.

Squared error makes the fit sensitive to unusual observations (NIST). Residual plots can show discrepancies between the fitted model and the data. A residual plot does not by itself prove that every modeling assumption holds. Leverage, Cook distance, and robust regression are left to later resources.

Classical interval estimates and tests use additional assumptions beyond the computing rule in Equation 4, such as correct linear specification, a zero conditional mean for the errors, and, for conventional homoscedastic standard errors, constant error variance and an appropriate independence structure. Those inferential procedures are not developed on this page. OLS is not described here as always BLUE.

In numerical software, least squares problems should be solved with established linear algebra routines rather than by forming $(X^{\mathsf{T}}X{)}^{-1}$ explicitly. The Python and MATLAB examples below use documented implementations. This page does not claim an exact internal solver except where the vendor documentation states it.

UCI Wine: alcohol is regressed on proline and flavanoids with an intercept. This is a two-predictor pedagogical fit, not a claimed chemical model. In-sample $R^2$ is 0.423. The Signature Visual plane is a separate geometric demonstration.

The statistical Python interface is statsmodels.api.OLS. The matrix interface does not add an intercept. Use sm.add_constant when an intercept is required. Fitted coefficients are results.params. Predictions for a new design matrix use results.predict.

The transparent numerical implementation is numpy.linalg.lstsq. It returns a least squares solution minimizing $\Vert y-X\beta {\Vert }_2$. If several minimizers exist, NumPy documents that the one with smallest $\Vert \beta {\Vert }_2$ is returned. The four outputs are the coefficients, the residual sum of squares when it is defined, the rank, and the singular values. This page does not implement np.linalg.inv(X.T @ X) @ X.T @ y as a recommended fit.

MATLAB statistical fitting uses fitlm(X, y). The default linear specification includes an intercept. Coefficient estimates are in mdl.Coefficients. Predictions for new rows use predict. MathWorks documents mdl.Fitted and mdl.Residuals on the LinearModel object. The transparent linear algebra form constructs a leading column of ones and uses backslash, Xdesign \\ y. For a rectangular design matrix, MathWorks documents that backslash returns a least squares solution. That solution is not necessarily the minimum-norm solution if the rank is smaller than the number of columns.

import numpy as npimport statsmodels.api as sm  def design_matrix(X, intercept=True):    X = np.asarray(X, dtype=float)    if X.ndim == 1:        X = X.reshape(-1, 1)    if X.ndim != 2:        raise ValueError("X must be a 1D or 2D array.")    if not np.all(np.isfinite(X)):        raise ValueError("X must contain only finite values.")    if intercept:        return sm.add_constant(X, prepend=True)    return X  def fit_ols(X, y, intercept=True):    y = np.asarray(y, dtype=float).reshape(-1)    if not np.all(np.isfinite(y)):        raise ValueError("y must contain only finite values.")    X_design = design_matrix(X, intercept=intercept)    if X_design.shape[0] != y.shape[0]:        raise ValueError("X and y must have the same number of observations.")    model = sm.OLS(y, X_design)    return model.fit()  def predict_ols(results, X_new, intercept=True):    X_design = design_matrix(X_new, intercept=intercept)    return np.asarray(results.predict(X_design), dtype=float)  def ols_lstsq(X, y, intercept=True):    y = np.asarray(y, dtype=float).reshape(-1)    if not np.all(np.isfinite(y)):        raise ValueError("y must contain only finite values.")    X_design = design_matrix(X, intercept=intercept)    if X_design.shape[0] != y.shape[0]:        raise ValueError("X and y must have the same number of observations.")    beta_hat, residuals, rank, singular_values = np.linalg.lstsq(        X_design,        y,        rcond=None,    )    y_hat = X_design @ beta_hat    e = y - y_hat    return beta_hat, y_hat, e, rank, singular_values 

1. 1.

   NIST/SEMATECH (n.d.). Linear least squares regression. e-Handbook of Statistical Methods, Section 4.1.4.1.

2. 2.

   NIST/SEMATECH (n.d.). Least squares. e-Handbook of Statistical Methods, Section 4.4.3.1.

3. 3.

   Haaland, D. M., & Thomas, E. V. (1988). Partial least-squares methods for spectral analyses. 1. Relation to other quantitative calibration methods and the extraction of qualitative information. Analytical Chemistry, 60(11), 1193-1202.

   doi:10.1021/ac00162a020
4. 4.

   statsmodels Developers (n.d.). statsmodels.regression.linear_model.OLS. statsmodels documentation.

5. 5.

   NumPy Developers (n.d.). numpy.linalg.lstsq. NumPy documentation.

6. 6.

   The MathWorks, Inc. (n.d.). fitlm. MATLAB documentation.

7. 7.

   The MathWorks, Inc. (n.d.). predict (LinearModel). MATLAB documentation.

8. 8.

   The MathWorks, Inc. (n.d.). mldivide (backslash). MATLAB documentation.