Multiple Linear Regression

Multiple Linear Regression

Multiple linear regression generalizes simple linear regression from one predictor to $p-1$ predictors per observation. With observations $\{({\mathbf{x}}_i,y_i){\}}_{i=1}^n$ where each ${\mathbf{x}}_i\in {\mathbb{R}}^{p-1}$ holds the predictor values for observation $i$, the model is
$y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_{p-1}{x}_{i,p-1}+{\epsilon }_i,$
or equivalently $y_i={\mathbf{x}}_i^{\top }\mathbf{\beta }+{\epsilon }_i$ with ${\mathbf{x}}_i=(1,x_{i1},\dots ,x_{i,p-1}{)}^{\top }$ absorbing the intercept and $\mathbf{\beta }=({\beta }_0,{\beta }_1,\dots ,{\beta }_{p-1}{)}^{\top }\in {\mathbb{R}}^p$. Stacking all $n$ observations into matrix form,
$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{\epsilon },$
where $\mathbf{y}\in {\mathbb{R}}^n$, $\mathbf{X}\in {\mathbb{R}}^{n\times p}$ is the design matrix whose row $i$ is ${\mathbf{x}}_i^{\top }$, and $\mathbf{\epsilon }\in {\mathbb{R}}^n$ is the vector of noise terms.

The OLS estimator is the closed-form minimizer of $\Vert \mathbf{y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2$:
$\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}.$
The inverse exists when the columns of $\mathbf{X}$ are linearly independent, which requires $p\le n$ at minimum. The geometric interpretation is that $\hat{\mathbf{y}}=\mathbf{X}\hat{\mathbf{\beta }}$ is the orthogonal projection of $\mathbf{y}$ onto the column space of $\mathbf{X}$, with the residual vector $\mathbf{y}-\hat{\mathbf{y}}$ perpendicular to every predictor.

Multicollinearity ${\mathbf{X}}^{\top }\mathbf{X}$ is close to singular and the OLS estimator becomes numerically unstable: small changes in $\mathbf{y}$ produce large swings in $\hat{\mathbf{\beta }}$, individual coefficient standard errors blow up, and individual coefficients can flip sign run-to-run while the joint fit remains stable. The standard remedies are dropping correlated predictors, pre-processing with PCA, or moving to ridge regression which adds $\lambda \mathbf{I}$ to the diagonal of ${\mathbf{X}}^{\top }\mathbf{X}$ and stabilizes the inverse.

When predictors are highly correlated,

Intuition

Multiple linear regression is the same idea as simple regression, scaled up: find the coefficient vector that makes the linear combination of predictors closest to the response in squared-error sense. The matrix formula $\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$ is the multidimensional version of the slope formula from simple regression, with $({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}$ playing the role of $1/\widehat{\mathrm{Var}}(x)$ and ${\mathbf{X}}^{\top }\mathbf{y}$ playing the role of $\widehat{\mathrm{Cov}}(x,y)$.

The geometric picture is even cleaner. The fitted vector $\hat{\mathbf{y}}$ lies in the column space of $\mathbf{X}$ (it is a linear combination of the predictors), and the residual vector $\mathbf{y}-\hat{\mathbf{y}}$ is the part of $\mathbf{y}$ that the predictors cannot explain. The OLS criterion picks the unique $\hat{\mathbf{y}}$ in the column space that minimizes the residual length, which by elementary geometry is the orthogonal projection of $\mathbf{y}$ onto the column space. The orthogonality condition ${\mathbf{X}}^{\top }(\mathbf{y}-\mathbf{X}\hat{\mathbf{\beta }})=0$ is exactly the matrix equation that produces the closed-form solution.

Properties

(OLS Solution) The closed-form estimator is $\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$ when $\mathbf{X}$ has full column rank.

The derivation comes from the orthogonality of the residual vector to the column space of $\mathbf{X}$. Setting ${\mathbf{X}}^{\top }(\mathbf{y}-\mathbf{X}\mathbf{\beta })=0$ gives the normal equations ${\mathbf{X}}^{\top }\mathbf{X}\mathbf{\beta }={\mathbf{X}}^{\top }\mathbf{y}$, and inverting the Gram matrix ${\mathbf{X}}^{\top }\mathbf{X}$ produces the solution. The Gram matrix is invertible exactly when the predictor columns are linearly independent, which requires $\mathrm{rank}(\mathbf{X})=p$ and consequently $n\ge p$. When $n<p$ or columns are collinear, the inverse fails to exist and the problem is ill-posed without regularization.

(Hat Matrix) The fitted values are a linear projection of the response.

Defining the hat matrix $\mathbf{H}=\mathbf{X}({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }$, the fitted values are $\hat{\mathbf{y}}=\mathbf{Hy}$, justifying the name (it puts the “hat” on $\mathbf{y}$). The hat matrix is the orthogonal projector onto the column space of $\mathbf{X}$: it is symmetric, idempotent (${\mathbf{H}}^2=\mathbf{H}$), and its trace equals $p$ (the number of fitted parameters, equal to the rank of the projection). The diagonal entries $h_{ii}$ are the leverages of each observation, measuring how much observation $i$‘s own response influences its own fitted value, and large leverages identify points whose removal would significantly change the fit.

(Gauss-Markov Theorem) OLS is the BLUE under uncorrelated, constant-variance noise.

If $\mathbb{E}[\mathbf{\epsilon }]=0$ and $\mathrm{Cov}(\mathbf{\epsilon })={\sigma }^2\mathbf{I}$, then among all linear unbiased estimators of $\mathbf{\beta }$, the OLS estimator has the smallest variance — it is the best linear unbiased estimator (BLUE). The result requires no distributional assumption beyond the first two moments. When the noise covariance is not scalar, the BLUE becomes the generalized least squares estimator ${\hat{\mathbf{\beta }}}_{\text{GLS}}=({\mathbf{X}}^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{y}$, which weights observations by the inverse noise covariance.

(Residual Variance) ${\hat{\sigma }}^2=\Vert \mathbf{y}-\mathbf{X}\hat{\mathbf{\beta }}{\Vert }^2/(n-p)$ is the unbiased estimator of ${\sigma }^2$.

The denominator is the degrees of freedom remaining after fitting $p$ parameters. The MLE under Gaussian noise is the corresponding biased estimator that divides by $n$; dividing by $n-p$ corrects for the lost degrees of freedom and produces an unbiased estimator. This is the variance estimator that goes into standard errors and $t$-statistics.

(Adjusted $R^2$) Penalizes added predictors that do not improve fit enough.

The ordinary $R^2=1-\mathrm{SSE}/\mathrm{SST}$ can only increase as predictors are added, since the optimization has more freedom. The adjusted version
$R_{\text{adj}}^2=1-\frac{\mathrm{SSE}/(n-p)}{\mathrm{SST}/(n-1)}=1-\frac{n-1}{n-p}(1-R^2)$
divides each sum of squares by its degrees of freedom, so adding a predictor that does not reduce SSE proportionally enough actually decreases $R_{\text{adj}}^2$. This makes it a more honest metric for comparing models of different sizes.

(Inference Under Normality) The coefficient vector is jointly Gaussian.

Under $\mathbf{\epsilon }\sim \mathcal{N}(0,{\sigma }^2\mathbf{I})$, the estimator is exactly Gaussian: $\hat{\mathbf{\beta }}\sim \mathcal{N}(\mathbf{\beta },{\sigma }^2({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1})$. Individual coefficients can be tested via $t_j={\hat{\beta }}_j/\mathrm{SE}({\hat{\beta }}_j)\sim t_{n-p}$ under $H_0:{\beta }_j=0$, and groups of coefficients via $F$-tests on nested models. Without Gaussianity, the same statistics are approximately Normal by the central limit theorem when $n$ is large.

(Spurious Regression) Regressing one trending variable on another can produce false significance.

When two unrelated variables both trend over time (or both have unit-root behavior), an OLS regression of one on the other often produces a high $R^2$ and small $p$-values even though the variables share no real dependence — the common trend is mistaken for a relationship. This is why blindly regressing time-series quantities is dangerous; the standard remedies are differencing the series, including a trend term explicitly, or testing for cointegration. See 2.4 Time Series Regression for the time-series-specific treatment.

(Iterative Computation by Gradient Descent) When $n$ or $p$ is too large for the closed form, gradient descent solves the same optimization.

For very large $n$ (millions of rows) or very large $p$ (thousands of features), the $O(n{p}^2+p^3)$ cost of forming and inverting ${\mathbf{X}}^{\top }\mathbf{X}$ becomes prohibitive. Iterative methods such as gradient descent
${\mathbf{\beta }}^{(t+1)}={\mathbf{\beta }}^{(t)}-\alpha {\mathbf{X}}^{\top }(\mathbf{X}{\mathbf{\beta }}^{(t)}-\mathbf{y})$
converge to the same OLS minimum at $O(np)$ per step and avoid the matrix inversion entirely. This is the algorithm of choice when fitting linear models on big data, and it is the same machinery (with a different loss) that drives logistic regression and most of deep learning.

Examples

Predicting house price from multiple features

$y_i$ is sale price; predictors are square footage, number of bedrooms, age, and distance to downtown ($p=5$ including intercept). Fit OLS on $n=1000$ houses.

Reading the fit $\hat{\mathbf{\beta }}$ has five entries, each interpreted as the marginal effect of a one-unit change in that predictor holding all others fixed. So ${\hat{\beta }}_{\text{bedrooms}}$ is the predicted price increase per added bedroom, conditional on square footage, age, and distance staying constant. The standard error of each coefficient depends on ${\mathbf{X}}^{\top }\mathbf{X}$; if bedrooms and square footage are highly correlated, both standard errors will be large because the data cannot cleanly attribute price increases to one or the other. The adjusted $R^2$ then tells you what fraction of price variance is explained, fairly accounting for the larger model.

The coefficient vector

Polynomial regression as a special case

Fit $y_i={\beta }_0+{\beta }_1{x}_i+{\beta }_2{x}_i^2+{\beta }_3{x}_i^3+{\epsilon }_i$ to a curved scatter plot.

Why this is still linear regression function of $x$ is non-linear, the parameters $({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3)$ enter linearly. Setting up the design matrix with columns $(1,x_i,x_i^2,x_i^3)$ recovers a multiple linear regression problem with $p=4$. The OLS solution and all the inferential machinery apply unchanged. This is the standard trick for fitting non-linear shapes with linear regression: pick a basis (polynomials, splines, sinusoids) and treat the basis evaluations as new predictors.

Although the

Geometric derivation of the normal equations

Show that $\hat{\mathbf{y}}$ being the orthogonal projection of $\mathbf{y}$ onto $\mathrm{span}(\mathbf{X})$ implies $\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$.

Solution $\mathrm{span}(\mathbf{X})$ such that $\mathbf{y}-\hat{\mathbf{y}}$ is perpendicular to every column of $\mathbf{X}$, i.e., ${\mathbf{X}}^{\top }(\mathbf{y}-\hat{\mathbf{y}})=0$. Substituting $\hat{\mathbf{y}}=\mathbf{X}\hat{\mathbf{\beta }}$ gives the normal equations ${\mathbf{X}}^{\top }\mathbf{X}\hat{\mathbf{\beta }}={\mathbf{X}}^{\top }\mathbf{y}$, and inverting the Gram matrix produces $\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}$. The geometric story is therefore not separate from the algebraic one: the orthogonality condition is the matrix equation that defines the OLS solution.

The orthogonal projection is the unique vector in

See Also

• 1.2. Linear Regression: the simple-regression starting point.
• 1.4. Nonlinear Regression: when the response depends on parameters non-linearly.
• 1.5. Regularization: ridge and lasso for ill-conditioned or high-dimensional designs.
• 1.6. Cross-Validation: model-selection machinery for choosing predictor sets.
• 2.4 Time Series Regression: applying MLR to time-indexed data with correlated errors.
• 3.8.1 Bayesian Linear Regression: putting a prior on $\mathbf{\beta }$ and reporting a posterior.