import numpy as np
import statsmodels.api as sm
rng = np.random.default_rng(0)
x = rng.uniform(0, 10, 200)
y = 1 + 2 * x + rng.normal(0, 1, 200)
fit = sm.OLS(y, sm.add_constant(x)).fit()
print(fit.params)[1.04711649 1.97469853]Linear Regression
Simple and multiple linear regression: model, OLS estimators, their properties, dummy variables and interactions.
1 Regression and the Linear Model
Given observations \((y_i,x_{i1},\dots,x_{ip})\), \(i=1,\dots,n\), of a scalar response (dependent variable) \(y\) and covariates (explanatory variables) \(x=(x_1,\dots,x_p)\), Regression Analysis studies the conditional expectation function
\[ \mu(x)=\mathbb{E}[Y\mid X=x]. \]
This function can have any shape, so it cannot be estimated without further assumptions. Linear regression assumes that \(\mu\) is affine in \(x\),
\[ \mu(x)=\beta_0+\beta_1x_1+\dots+\beta_px_p, \]
which allows closed-form parameter estimates by ordinary least squares (OLS), see Least Squares Estimation. With \(p=1\) it is called simple linear regression (SLR), and with \(p>1\) multiple linear regression (MLR). Both differ from multivariate linear regression, where several correlated responses are predicted. Regression models the conditional distribution of the response given the covariates rather than the joint distribution of all variables, which is the subject of multivariate analysis.
2 The Model
For each statistical unit \(i=1,\dots,n\),
\[ y_i=\beta_0+\beta_1x_{i1}+\dots+\beta_px_{ip}+\varepsilon_i =\mathbf{x}_i^\top\boldsymbol{\beta}+\varepsilon_i, \]
where \(\beta_0\) is the intercept, \(\beta_1,\dots,\beta_p\) are the slopes (the effect on \(\mu\) of a unit change in \(x_j\), other covariates held fixed) and \(\varepsilon_i\) is the error term that adds noise to the relationship. In matrix form
\[ \mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}, \]
where
\[ \mathbf{y}=\begin{bmatrix}y_1\\\vdots\\y_n\end{bmatrix},\quad \mathbf{X}=\begin{bmatrix}1&x_{11}&\cdots&x_{1p}\\ \vdots&\vdots&&\vdots\\1&x_{n1}&\cdots&x_{np}\end{bmatrix},\quad \boldsymbol{\beta}=\begin{bmatrix}\beta_0\\\vdots\\\beta_p\end{bmatrix},\quad \boldsymbol{\varepsilon}=\begin{bmatrix}\varepsilon_1\\\vdots\\\varepsilon_n \end{bmatrix}. \]
The standard assumptions are
- Linearity and weak exogeneity: \(\mathbb{E}[\varepsilon\mid X]=0\),
- Constant variance (homoscedasticity): \(\mathrm{Var}(\varepsilon_i)=\sigma^2\),
- Independence (or uncorrelatedness) of the errors,
- No perfect multicollinearity: \(\mathbf{X}\) has full column rank.
For exact finite-sample inference one also assumes \(\boldsymbol{\varepsilon}\sim\mathcal{N}(0,\sigma^2I_n)\), see Normal Distribution.
3 Simple Linear Regression
For \(p=1\) the model is \(y_i=\beta_0+\beta_1x_i+\varepsilon_i\).
The ordinary least squares estimators minimise the residual sum of squares
\[ \mathrm{RSS}(\beta_0,\beta_1)=\sum_{i=1}^n(y_i-\beta_0-\beta_1x_i)^2. \]
Setting the partial derivatives to zero gives the normal equations
\[ \sum_i(y_i-\beta_0-\beta_1x_i)=0,\qquad \sum_ix_i(y_i-\beta_0-\beta_1x_i)=0, \]
with solution
\[ \hat\beta_1=\frac{\sum_i(x_i-\bar x)(y_i-\bar y)}{\sum_i(x_i-\bar x)^2} =\frac{S_{xy}}{S_{xx}},\qquad \hat\beta_0=\bar y-\hat\beta_1\bar x, \]
where \(\bar x,\bar y\) are the sample means. The fitted line passes through \((\bar x,\bar y)\), and the slope is the sample covariance of \(x,y\) over the sample variance of \(x\).
4 OLS in General
For general \(p\), minimising \(\mathrm{RSS}(\boldsymbol\beta)=\|\mathbf{y}-\mathbf{X}\boldsymbol\beta\|^2\) and setting the gradient \(-2\mathbf{X}^\top(\mathbf{y}-\mathbf{X}\boldsymbol\beta)\) to zero gives the normal equations \(\mathbf{X}^\top\mathbf{X}\boldsymbol\beta= \mathbf{X}^\top\mathbf{y}\), hence
\[ \hat{\boldsymbol\beta}=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}. \]
The matrix \(\mathbf{X}^\top\mathbf{X}\) is invertible exactly when there is no perfect multicollinearity. Geometrically, the fitted values \(\hat{\mathbf{y}}=\mathbf{X}\hat{\boldsymbol\beta}=H\mathbf{y}\) with \(H=\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\) are the orthogonal projection of \(\mathbf{y}\) onto the column space of \(\mathbf{X}\) (see Orthogonality and Projections), and the residuals \(\hat{\boldsymbol\varepsilon}=\mathbf{y}-\hat{\mathbf{y}}\) are orthogonal to every column of \(\mathbf{X}\). An unbiased estimator of the error variance is
\[ \hat\sigma^2=\frac{\mathrm{RSS}}{n-p-1}. \]
5 Properties of the Estimators
Under assumptions 1-4, \(\hat{\boldsymbol\beta}\) is unbiased, \(\mathbb{E}[\hat{\boldsymbol\beta}]=\boldsymbol\beta\), with covariance matrix
\[ \mathrm{Var}(\hat{\boldsymbol\beta})=\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}, \]
and it has the smallest covariance among all linear unbiased estimators (it is the BLUE).
Proof. Since \(\hat{\boldsymbol\beta}=\boldsymbol\beta+(\mathbf{X}^\top\mathbf{X})^{-1} \mathbf{X}^\top\boldsymbol\varepsilon\), taking expectations gives unbiasedness and taking covariances gives the variance formula. For another linear unbiased estimator \(\tilde{\boldsymbol\beta}=C\mathbf{y}\) write \(C=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top+D\). Unbiasedness forces \(D\mathbf{X}=0\), so \(\mathrm{Var}(\tilde{\boldsymbol\beta})=\sigma^2CC^\top =\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}+\sigma^2DD^\top\), and \(DD^\top\) is positive semi-definite.
For SLR this gives
\[ \mathrm{Var}(\hat\beta_1)=\frac{\sigma^2}{S_{xx}},\qquad \mathrm{Var}(\hat\beta_0)=\sigma^2\left(\frac1n+\frac{\bar x^2}{S_{xx}}\right). \]
If in addition the errors are normal, then \(\hat{\boldsymbol\beta}\sim \mathcal{N}(\boldsymbol\beta,\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1})\) (it is also the maximum likelihood estimator), and \((n-p-1)\hat\sigma^2/\sigma^2\sim\chi^2_{n-p-1}\) independently of \(\hat{\boldsymbol\beta}\). Each standardised coefficient \((\hat\beta_j-\beta_j)/\widehat{\mathrm{se}}(\hat\beta_j)\) then follows a [[students-t-distribution|Student \(t\) distribution]] with \(n-p-1\) degrees of freedom, giving confidence intervals and tests, and joint hypotheses use the F-Distribution.
6 Multiple Regression: Categorical Predictors and Interactions
In MLR each slope \(\beta_j\) is interpreted holding the other covariates fixed, which is why correlated covariates make estimates unstable (large \(\mathrm{Var}(\hat\beta_j)\)).
Dummy variables. A categorical covariate with \(L\) levels is encoded by \(L-1\) indicator variables \(d_\ell=\mathbb{1}\{\text{level}=\ell\}\), with one level kept as the baseline. Using all \(L\) indicators together with an intercept would give perfect multicollinearity. The coefficient on \(d_\ell\) is the shift in the mean response relative to the baseline level.
Interaction terms. To let the effect of one covariate depend on another, add their product,
\[ y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_1x_2+\varepsilon, \]
so that the slope in \(x_1\) is \(\beta_1+\beta_3x_2\). Interacting a dummy with a continuous covariate gives a separate slope for each category.
Simulate \(y=1+2x+\varepsilon\) and recover the coefficients with OLS.
Residual diagnostics (residuals against fitted values, normal Q-Q plot) check the constant variance, linearity and normality assumptions.
Choosing which covariates to include is the subject of Model Selection. Responses that are binary or counts require the generalisation to the Generalized Linear Models, for instance Logistic Regression and Poisson Regression.
7 References
- Hastie, T., Tibshirani, R. and Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed. Springer, Chapter 3.
- Wooldridge, J. M. (2019). Introductory Econometrics: A Modern Approach, 7th ed. Cengage.