Correlation and regression analysis

Intuition Beginner

Two variables are correlated when they tend to vary together. Taller people tend to weigh more. Countries with higher education spending tend to have higher GDP. Hours of study tend to be associated with exam scores. Correlation measures the strength and direction of the linear relationship between two quantitative variables.

The most common measure of correlation is the Pearson correlation coefficient $r$, which ranges from $-1$ to $+1$. A correlation of $+1$ means a perfect positive linear relationship: as one variable increases, the other increases proportionally, and all data points fall exactly on a line with positive slope. A correlation of $-1$ means a perfect negative linear relationship. A correlation near 0 means no linear relationship: knowing the value of one variable tells you nothing about the other.

Correlation does not imply causation. This is perhaps the most repeated warning in all of statistics. Ice cream sales and drowning deaths are positively correlated, but eating ice cream does not cause drowning. Both increase in hot weather. The lurking variable (temperature) explains the correlation. Confusing correlation with causation is one of the most common and most dangerous statistical fallacies.

Regression analysis extends correlation by fitting a line (or curve) to the data. The regression line describes how the dependent variable $y$ (the response) changes as the independent variable $x$ (the predictor) changes. The most common form is simple linear regression, which fits the model $y=a+bx$ to the data by choosing $a$ (intercept) and $b$ (slope) to minimise the sum of squared residuals. This is the method of least squares.

A residual is the vertical distance from each data point to the regression line: $e_i=y_i-(a+b{x}_i)$. The least squares line minimises the total of the squared residuals. It is the line that makes the residuals as small as possible in the sense of sum of squared vertical distances. The least squares method was developed independently by Gauss and Legendre in the early nineteenth century for astronomical calculations.

The slope $b$ of the regression line has a direct interpretation: for each one-unit increase in $x$, the model predicts a $b$-unit change in $y$. If the regression of blood pressure on age has slope 0.8, the model predicts that blood pressure increases by 0.8 mm Hg for each additional year of age. The intercept $a$ is the predicted value of $y$ when $x=0$. The intercept is meaningful only if $x=0$ is within or near the range of the data.

The coefficient of determination $R^2$ measures the fraction of the variability in $y$ that is explained by the regression on $x$. $R^2=1$ means the regression line fits the data perfectly (all points on the line). $R^2=0$ means the regression line explains none of the variability. In simple linear regression, $R^2=r^2$: the square of the correlation coefficient gives the proportion of variance explained.

Regression makes predictions, but extrapolation beyond the range of the data is risky. A regression model fitted to adults aged 20-60 may predict nonsensical values for children or centenarians. The linear relationship that holds within the observed range may not extend outside it. Extrapolation is one of the most common sources of prediction error.

Multiple regression extends the idea to several predictors simultaneously. Instead of fitting a line in two dimensions, you fit a plane (with two predictors) or a hyperplane (with more) through the data. The model $y=b_0+b_1{x}_1+b_2{x}_2+\dots +b_k{x}_k$ predicts $y$ from $k$ predictors. Each coefficient $b_j$ represents the change in $y$ associated with a one-unit increase in $x_j$, holding all other predictors constant. This "holding constant" interpretation is what makes multiple regression so powerful: it approximates the effect of a controlled experiment using observational data.

The "holding constant" interpretation has limits. In observational data, the predictors are often correlated with each other (multicollinearity), which makes it difficult to isolate the effect of any single predictor. When two predictors are highly correlated, the data contain little information about their separate effects, and the regression coefficients become unstable (small changes in the data produce large changes in the coefficients). Multicollinearity does not violate any regression assumption, but it makes the coefficients hard to interpret and inflates their standard errors.

Regression diagnostics are checks on the assumptions of the regression model. The most important assumptions are linearity (the relationship is linear), independence (the residuals are independent), homoscedasticity (the residuals have constant variance), and normality (the residuals are normally distributed). Diagnostic plots (residuals versus fitted values, normal Q-Q plots, scale-location plots, leverage plots) reveal violations of these assumptions. A good regression analysis always includes a diagnostic check: fitting the model is only half the work, and checking it is the other half.

Visual Beginner

Concept	Formula	Interpretation
Correlation $r$	$r=\frac{S_{XY}}{\sqrt{S_{XX}\cdot S_{YY}}}$	Strength and direction of linear relationship
Slope $b$	$b=r\frac{s_y}{s_x}$	Change in $y$ per unit change in $x$
Intercept $a$	$a=\bar{y}-b\bar{x}$	Predicted $y$ when $x=0$
$R^2$	$R^2=r^2$ (simple regression)	Fraction of variance in $y$ explained by $x$

The scatter plot is the fundamental visual tool for studying the relationship between two quantitative variables. The pattern of points reveals the direction (positive or negative), form (linear or curved), and strength (tight or loose) of the relationship, as well as any unusual observations (outliers or influential points).

Worked example Beginner

A statistics instructor collects data on study hours ($x$) and exam scores ($y$) from 8 students:

Student	$x$ (hours)	$y$ (score)
1	2	55
2	3	62
3	4	68
4	5	73
5	6	78
6	7	82
7	8	88
8	9	94

Computing the summaries: $\bar{x}=5.5$, $\bar{y}=75$, $s_x=2.45$, $s_y=13.1$.

The correlation is $r=0.994$, indicating an extremely strong positive linear relationship.

The slope is $b=r\cdot s_y/s_x=0.994\times 13.1/2.45=5.31$.

The intercept is $a=\bar{y}-b\bar{x}=75-5.31(5.5)=75-29.2=45.8$.

The regression line is $\hat{y}=45.8+5.31x$.

Interpretation: each additional hour of study is associated with a predicted increase of 5.31 points on the exam.

$R^2=0.994^2=0.988$. About 98.8% of the variability in exam scores is explained by study hours.

For a student who studies 6.5 hours, the predicted score is $\hat{y}=45.8+5.31(6.5)=80.3$.

The residuals (observed minus predicted) provide a check on the model. Student 2 studied 3 hours and scored 62; the prediction is $45.8+5.31(3)=61.7$, so the residual is $62-61.7=0.3$. Student 7 studied 8 hours and scored 88; the prediction is $45.8+5.31(8)=88.3$, so the residual is $88-88.3=-0.3$. The residuals are all small (the largest in absolute value is about 1.5 points), confirming that the linear model fits well.

This example is unusually clean. Real data rarely produce correlations of 0.994 or $R^2$ values of 0.988. In practice, social science data typically yields $R^2$ values between 0.1 and 0.5, and even well-fitting physical science models rarely exceed 0.95. The high $R^2$ here reflects the fact that study hours is a strong predictor of exam performance and the data were constructed to be nearly linear.

A word of caution about this regression: it does not prove that studying causes higher exam scores. It is plausible that students who study more are also more motivated, better prepared, or more knowledgeable about effective study techniques. The regression controls for none of these confounders. The coefficient 5.31 is the association between study hours and exam scores, not the causal effect of an additional hour of study. Establishing causation would require randomised assignment of study hours, which is impractical.

Check your understanding Beginner

Formal definition Intermediate+

The linear regression model

The simple linear regression model is:

$Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i,i=1,\dots ,n$

where ${\beta }_0$ is the intercept, ${\beta }_1$ is the slope, and ${ϵ}_i$ are independent random errors with $E[{ϵ}_i]=0$ and $\text{Var}({ϵ}_i)={\sigma }^2$.

Ordinary least squares estimators

The OLS estimators minimise the residual sum of squares $\text{RSS}={\sum }_{i=1}^n(Y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{X}_i{)}^2$.

${\hat{\beta }}_1=\frac{{\sum }_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{{\sum }_{i=1}^n(X_i-\bar{X}{)}^2}=\frac{S_{XY}}{S_{XX}}$

${\hat{\beta }}_0=\bar{Y}-{\hat{\beta }}_1\bar{X}$

Properties of OLS estimators

Under the regression model assumptions, the OLS estimators have the following properties.

Unbiasedness. $E[{\hat{\beta }}_0]={\beta }_0$ and $E[{\hat{\beta }}_1]={\beta }_1$.

Variances. $\text{Var}({\hat{\beta }}_1)={\sigma }^2/S_{XX}$ and $\text{Var}({\hat{\beta }}_0)={\sigma }^2(1/n+{\bar{X}}^2/S_{XX})$, where $S_{XX}=\sum (X_i-\bar{X}{)}^2$.

Gauss-Markov theorem. Under the assumptions of linearity, independence, homoscedasticity ($\text{Var}({ϵ}_i)={\sigma }^2$ for all $i$), and zero mean of errors, the OLS estimators are BLUE: Best Linear Unbiased Estimators. "Best" means they have the smallest variance among all linear unbiased estimators.

The coefficient of determination

The total sum of squares decomposes as $\text{TSS}=\text{RSS}+\text{RegSS}$:

$\sum (Y_i-\bar{Y}{)}^2=\sum (Y_i-{\hat{Y}}_i{)}^2+\sum ({\hat{Y}}_i-\bar{Y}{)}^2$

$R^2=1-\frac{\text{RSS}}{\text{TSS}}=\frac{\text{RegSS}}{\text{TSS}}$

Inference for regression coefficients

Under the additional assumption that ${ϵ}_i\sim N(0,{\sigma }^2)$ (normal errors), the OLS estimators are normally distributed:

${\hat{\beta }}_1\sim N\left({\beta }_1,\frac{{\sigma }^2}{S_{XX}}\right),{\hat{\beta }}_0\sim N\left({\beta }_0,{\sigma }^2\left(\frac{1}{n}+\frac{{\bar{X}}^2}{S_{XX}}\right)\right)$

Replacing ${\sigma }^2$ with the unbiased estimate ${\hat{\sigma }}^2=\text{RSS}/(n-2)$ gives t-distributed test statistics:

$t=\frac{{\hat{\beta }}_1-{\beta }_{1,0}}{\text{SE}({\hat{\beta }}_1)}\sim t_{n-2}$

Multiple regression

The multiple regression model is:

$Y_i={\beta }_0+{\beta }_1{X}_{i1}+{\beta }_2{X}_{i2}+\cdots +{\beta }_p{X}_{ip}+{ϵ}_i$

In matrix notation: $\mathbf{Y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$, where $\mathbf{Y}$ is $n\times 1$, $\mathbf{X}$ is $n\times (p+1)$, $\mathbf{\beta }$ is $(p+1)\times 1$, and $\mathbf{ϵ}$ is $n\times 1$.

The OLS estimator is $\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{Y}$, assuming ${\mathbf{X}}^{\top }\mathbf{X}$ is invertible.

The fitted values are $\hat{\mathbf{Y}}=\mathbf{X}\hat{\mathbf{\beta }}=\mathbf{HY}$, where $\mathbf{H}=\mathbf{X}({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }$ is the hat matrix (the projection matrix onto the column space of $\mathbf{X}$).

$\text{Var}(\hat{\mathbf{\beta }})={\sigma }^2({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}$.

Adjusted R-squared and model comparison

Adding predictors always increases $R^2$, even if the predictors are unrelated to the response. The adjusted $R^2$ penalises model complexity:

$R_{\text{adj}}^2=1-\frac{\text{RSS}/(n-p-1)}{\text{TSS}/(n-1)}$

The adjustment replaces the total sum of squares with the mean square error, which accounts for the degrees of freedom used by the model. Unlike $R^2$, the adjusted $R^2$ can decrease when a predictor is added, if the predictor does not improve the fit enough to compensate for the lost degree of freedom.

Prediction intervals

A confidence interval for the mean response at $x=x_0$ is narrower than a prediction interval for an individual response. The confidence interval for $E[Y|X=x_0]$ accounts for the uncertainty in estimating the regression line. The prediction interval for a new observation at $x_0$ must also account for the variability of the individual observation around the line.

The prediction interval is: ${\hat{y}}_0\pm t_{\alpha /2,n-2}\cdot s\sqrt{1+1/n+(x_0-\bar{x}{)}^2/S_{XX}}$. The "1" under the square root represents the individual variability, and the remaining terms represent the uncertainty in the estimated mean. The prediction interval is always wider than the confidence interval, and it widens as $x_0$ moves away from $\bar{x}$, reflecting the increased uncertainty of extrapolation.

The adjusted $R^2$ accounts for the number of predictors: $R_{\text{adj}}^2=1-\frac{\text{RSS}/(n-p-1)}{\text{TSS}/(n-1)}$.

Diagnostics

The validity of regression inference depends on several assumptions: linearity, independence of errors, constant variance (homoscedasticity), and normality of errors. Diagnostic tools include residual plots (plotting residuals against fitted values or predictors to detect non-linearity and heteroscedasticity), normal probability plots of residuals, leverage values $h_{ii}$ (diagonal elements of the hat matrix), and Cook's distance $D_i=\frac{(\hat{\mathbf{Y}}-{\hat{\mathbf{Y}}}_{(i)}{)}^{\top }(\hat{\mathbf{Y}}-{\hat{\mathbf{Y}}}_{(i)})}{p{\hat{\sigma }}^2}$ measuring the influence of observation $i$ on the fitted values.

Key theorem with proof Intermediate+

The Gauss-Markov theorem

Theorem. Under the linear model $Y_i={\mathbf{x}}_i^{\top }\mathbf{\beta }+{ϵ}_i$ with $E[{ϵ}_i]=0$, $\text{Var}({ϵ}_i)={\sigma }^2$, and $\text{Cov}({ϵ}_i,{ϵ}_j)=0$ for $i\ne j$, the OLS estimator $\hat{\mathbf{\beta }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{Y}$ is BLUE: among all linear unbiased estimators of $\mathbf{\beta }$, it has the smallest variance.

Proof. Let $\tilde{\mathbf{\beta }}=\mathbf{CY}$ be any linear unbiased estimator, where $\mathbf{C}$ is a $(p+1)\times n$ matrix. Unbiasedness requires $E[\mathbf{CY}]=\mathbf{CX}\mathbf{\beta }=\mathbf{\beta }$ for all $\mathbf{\beta }$, so $\mathbf{CX}=\mathbf{I}$.

Write $\mathbf{C}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }+\mathbf{D}$, where $\mathbf{D}=\mathbf{C}-({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }$. Then $\mathbf{CX}=\mathbf{I}+\mathbf{DX}=\mathbf{I}$ requires $\mathbf{DX}=0$.

The covariance matrix of $\tilde{\mathbf{\beta }}$ is:

$\text{Var}(\tilde{\mathbf{\beta }})={\sigma }^2\mathbf{C}{\mathbf{C}}^{\top }={\sigma }^2\left[({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }+\mathbf{D}\right]\left[\mathbf{X}({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}+{\mathbf{D}}^{\top }\right]$

$={\sigma }^2\left[({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}+\mathbf{D}{\mathbf{D}}^{\top }\right]$

using $\mathbf{DX}=0$. Since $\mathbf{D}{\mathbf{D}}^{\top }$ is positive semi-definite, $\text{Var}(\tilde{\mathbf{\beta }})\ge {\sigma }^2({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}=\text{Var}(\hat{\mathbf{\beta }})$ in the matrix sense. $\square$

Maximum likelihood estimation under normality

Theorem. Under $Y_i={\beta }_0+{\beta }_1{X}_i+{ϵ}_i$ with ${ϵ}_i\sim N(0,{\sigma }^2)$, the MLE of $({\beta }_0,{\beta }_1,{\sigma }^2)$ coincides with OLS for the regression coefficients and ${\hat{\sigma }}_{\text{MLE}}^2=\text{RSS}/n$.

Proof sketch. The log-likelihood is:

$\ell ({\beta }_0,{\beta }_1,{\sigma }^2)=-\frac{n}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}{\sum }_{i=1}^n(Y_i-{\beta }_0-{\beta }_1{X}_i{)}^2$

Maximising over $({\beta }_0,{\beta }_1)$ is equivalent to minimising $\sum (Y_i-{\beta }_0-{\beta }_1{X}_i{)}^2$, which is the OLS criterion. The maximising values of ${\sigma }^2$ given the OLS coefficients is ${\hat{\sigma }}^2=\text{RSS}/n$.

Exercises Intermediate+

Advanced results Master

The geometry of least squares

The OLS solution has an elegant geometric interpretation. The vector of fitted values $\hat{\mathbf{Y}}$ is the orthogonal projection of the response vector $\mathbf{Y}$ onto the column space of the design matrix $\mathbf{X}$. The hat matrix $\mathbf{H}=\mathbf{X}({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }$ is the projection matrix: it is symmetric ($\mathbf{H}={\mathbf{H}}^{\top }$), idempotent (${\mathbf{H}}^2=\mathbf{H}$), and has trace equal to $p+1$ (the number of parameters).

The residual vector $\mathbf{e}=(\mathbf{I}-\mathbf{H})\mathbf{Y}$ is orthogonal to the column space of $\mathbf{X}$. This orthogonality is the geometric expression of the normal equations: the residuals are uncorrelated with every column of $\mathbf{X}$.

The eigenvalues of $\mathbf{H}$ are all 0 or 1. The rank of $\mathbf{H}$ equals the number of linearly independent columns of $\mathbf{X}$, which is the number of parameters in the model. This geometric perspective connects regression to the theory of orthogonal projections in linear algebra.

Generalised least squares and weighted regression

When the errors have non-constant variance (heteroscedasticity) or are correlated, ordinary least squares is no longer optimal. Generalised least squares (GLS) accounts for the covariance structure of the errors. If $\text{Var}(\mathbf{ϵ})={\sigma }^2\mathbf{V}$, where $\mathbf{V}$ is a known positive-definite matrix, the GLS estimator is ${\hat{\mathbf{\beta }}}_{\text{GLS}}=({\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{y}$. The GLS estimator is the BLUE under the general Gauss-Markov theorem.

Weighted least squares (WLS) is the special case where $\mathbf{V}$ is diagonal (errors are uncorrelated but have different variances). The weight for observation $i$ is $w_i=1/{\sigma }_i^2$, giving more weight to observations with smaller variance (more precise measurements). WLS is commonly used when the response variable is an average of varying numbers of observations, when the variability increases with the level of the response, or when the data are from a stratified sample.

Feasible GLS (FGLS) replaces the unknown $\mathbf{V}$ with an estimate $\hat{\mathbf{V}}$ obtained from the OLS residuals. The iterative procedure (estimate $\mathbf{V}$, compute GLS estimates, re-estimate $\mathbf{V}$, repeat) converges to the GLS estimator under mild conditions. FGLS is widely used in econometrics for models with serial correlation (Cochrane-Orcutt procedure) or panel data with group-specific variances.

Weighted least squares and generalised least squares

When the errors have non-constant variance (heteroscedasticity) or are correlated, OLS is no longer efficient. Weighted least squares (WLS) handles heteroscedasticity by minimising $\sum w_i(Y_i-{\mathbf{x}}_i^{\top }\mathbf{\beta }{)}^2$ where $w_i=1/\text{Var}({ϵ}_i)$. Observations with smaller variance receive larger weights, contributing more to the fit.

Generalised least squares (GLS) handles correlated errors by transforming the model. If $\text{Var}(\mathbf{ϵ})={\sigma }^2\mathbf{V}$ where $\mathbf{V}$ is a known positive-definite matrix, the GLS estimator is ${\hat{\mathbf{\beta }}}_{\text{GLS}}=({\mathbf{X}}^{\top }{\mathbf{V}}^{-1}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }{\mathbf{V}}^{-1}\mathbf{Y}$. This is the BLUE under the general covariance structure.

Feasible GLS (FGLS) replaces $\mathbf{V}$ with an estimate $\hat{\mathbf{V}}$ when the covariance structure is unknown but can be estimated from the data. Iteratively reweighted least squares (IRLS) alternates between estimating the weights and fitting the model, converging to the GLS solution.

Regularisation: ridge regression and the lasso

When the design matrix has collinear columns or the number of predictors approaches the sample size, the OLS estimator becomes unstable (high variance). Regularisation methods address this by adding a penalty to the least squares criterion.

Ridge regression (Tikhonov regularisation) minimises $\Vert \mathbf{Y}-\mathbf{X}\mathbf{\beta }{\Vert }^2+\lambda \Vert \mathbf{\beta }{\Vert }^2$, where $\lambda >0$ is a tuning parameter. The ridge estimator is ${\hat{\mathbf{\beta }}}_{\text{ridge}}=({\mathbf{X}}^{\top }\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{Y}$. The penalty shrinks the coefficients toward zero, trading bias for reduced variance. Ridge regression is equivalent to a Bayesian posterior mean under a normal prior on the coefficients.

The lasso (least absolute shrinkage and selection operator) minimises $\Vert \mathbf{Y}-\mathbf{X}\mathbf{\beta }{\Vert }^2+\lambda \sum |{\beta }_j|$. The lasso penalty produces sparse solutions: some coefficients are shrunk exactly to zero, performing automatic variable selection. The lasso is computationally harder than ridge (no closed-form solution) but produces more interpretable models.

Elastic net combines ridge and lasso penalties: $\Vert \mathbf{Y}-\mathbf{X}\mathbf{\beta }{\Vert }^2+{\lambda }_1\Vert \mathbf{\beta }{\Vert }^2+{\lambda }_2\sum |{\beta }_j|$. This handles correlated predictors better than the lasso alone.

Generalised linear models

The linear regression model assumes that the response is continuous and normally distributed. Generalised linear models (GLMs) extend regression to non-normal responses by specifying a link function $g$ that connects the mean of $Y$ to the linear predictor: $g(E[Y])=\mathbf{X}\mathbf{\beta }$.

Logistic regression models binary responses using the logit link: $\log (p/(1-p))=\mathbf{X}\mathbf{\beta }$, where $p=P(Y=1)$. The response follows a Bernoulli distribution. Logistic regression is the workhorse of classification in statistics.

Poisson regression models count data using the log link: $\log (\lambda )=\mathbf{X}\mathbf{\beta }$, where $\lambda =E[Y]$. The response follows a Poisson distribution. Poisson regression is used for modelling rates and frequencies.

GLMs are fitted by iteratively reweighted least squares, which solves the maximum likelihood equations. The theory of GLMs, developed by Nelder and Wedderburn in 1972, unifies regression, logistic regression, Poisson regression, and many other models under a single framework.

Nonparametric regression

Parametric regression assumes a specific functional form (linear, polynomial). Nonparametric regression lets the data determine the shape of the relationship. Kernel regression estimates the conditional mean $E[Y|X=x]$ as a weighted average of nearby observations, with weights determined by a kernel function. Local polynomial regression (loess) fits low-degree polynomials locally using weighted least squares.

Splines fit piecewise polynomials with smoothness constraints. Regression splines use a fixed set of knots; smoothing splines penalise roughness by minimising $\sum (Y_i-f(x_i){)}^2+\lambda \int f^{\prime \prime }(x{)}^{2}dx$, trading fit against smoothness. The effective degrees of freedom are controlled by $\lambda$.

Causal inference and regression

Regression can estimate causal effects under strict conditions: random assignment of the treatment, no unmeasured confounders, and correct model specification. In observational studies, regression adjustments can reduce confounding bias but cannot eliminate it if unmeasured confounders exist.

The potential outcomes framework (Rubin causal model) defines the causal effect of a treatment as $Y_i(1)-Y_i(0)$, the difference between potential outcomes under treatment and control. Since only one potential outcome is observed for each unit, individual causal effects are not identified. Regression estimates the average treatment effect under the assumption of ignorable treatment assignment.

Instrumental variable regression addresses endogeneity (correlation between the predictor and the error term) by using an instrument $Z$ that affects $Y$ only through $X$. Two-stage least squares first regresses $X$ on $Z$ to obtain predicted values $\hat{X}$, then regresses $Y$ on $\hat{X}$. The instrumental variable approach identifies the causal effect of $X$ on $Y$ even in the presence of unmeasured confounders, provided the instrument satisfies the exclusion restriction.

Model selection and information criteria

Choosing between regression models with different sets of predictors requires balancing goodness of fit against model complexity. Adding predictors always improves $R^2$, but the improvement may be due to overfitting rather than genuine predictive power.

Akaike's Information Criterion (AIC) trades off fit against complexity: $\text{AIC}=-2\log L+2p$, where $L$ is the maximised likelihood and $p$ is the number of parameters. Lower AIC indicates a better model. AIC estimates the expected Kullback-Leibler divergence between the model and the true distribution.

The Bayesian Information Criterion (BIC) imposes a stronger penalty for complexity: $\text{BIC}=-2\log L+p\log n$. For large $n$, BIC selects simpler models than AIC. BIC is consistent (it selects the true model with probability approaching 1 as $n\to \infty$) while AIC is efficient (it selects the model with the best predictive performance). The choice between AIC and BIC reflects a philosophical choice between prediction and truth.

Mallows' $C_p$ statistic provides another model selection criterion: $C_p={\text{RSS}}_p/{\hat{\sigma }}^2-n+2p$. A good model has $C_p\approx p$. The $C_p$ criterion is equivalent to AIC for linear models with normal errors.

Cross-validation provides a direct estimate of predictive performance by fitting the model on a training set and evaluating on a held-out test set. K-fold cross-validation divides the data into $k$ folds, fits on $k-1$ folds, and evaluates on the held-out fold, repeating for each fold. Leave-one-out cross-validation (LOOCV) is the special case $k=n$, which is computationally expensive but provides an nearly unbiased estimate of prediction error.

Connections Master

• Descriptive statistics 26.01.01. The regression line is a summary of the bivariate relationship between $X$ and $Y$, analogous to how the mean summarises a univariate distribution. The correlation coefficient is a standardised covariance.

• Sampling distributions 26.04.01. The regression coefficients are statistics with sampling distributions. The normality of these distributions (under normal errors) enables hypothesis tests and confidence intervals for the coefficients.

• Hypothesis testing 26.05.01. The t-test for the slope and the F-test for the overall model are applications of the hypothesis testing framework to regression parameters.

• Experimental design 26.09.01. ANOVA is a special case of regression with categorical predictors. The F-test in ANOVA is equivalent to testing whether the regression coefficients for the group indicators are all zero.

• Bayesian statistics 26.07.01. Bayesian regression places prior distributions on the regression coefficients and computes the posterior distribution. Ridge regression corresponds to a normal prior on the coefficients.

• Linear algebra 01.01.09. OLS regression is a projection onto the column space of the design matrix. The hat matrix, leverage values, and residual analysis are all grounded in the theory of orthogonal projections.

• Machine learning and data science. Regression is the foundation of supervised learning. Linear regression, logistic regression, and regularised regression are among the most widely used predictive models. The bias-variance trade-off in regularised regression connects to the broader theory of model selection.

• Econometrics. Instrumental variables, panel data regression, and time series regression are extensions of the basic regression model developed to handle the specific challenges of economic data (endogeneity, serial correlation, unobserved heterogeneity).

• Nonparametric methods 26.08.01. Nonparametric regression relaxes the linearity assumption, using kernel smoothing, splines, or local polynomials to estimate the conditional mean. The bias-variance trade-off in choosing the bandwidth or smoothing parameter is analogous to model selection in parametric regression.

• Statistical literacy 26.10.01. Misinterpretation of regression results (confusing correlation with causation, extrapolating beyond the data, ignoring confounders) is one of the most common statistical errors. Understanding the assumptions and limitations of regression is essential for responsible data analysis.

Historical and philosophical context Master

Galton and regression toward the mean

Francis Galton discovered the concept of regression in 1886 while studying the relationship between the heights of parents and their children. Galton observed that tall parents tended to have children who were shorter than themselves (though still above average) and short parents tended to have children who were taller. He called this phenomenon "regression toward mediocrity" and later "regression toward the mean."

Galton's insight was that the relationship between parent and child heights was not a line of slope 1 (perfect inheritance) but a line with slope less than 1 (partial regression). This statistical phenomenon, not a biological force, is a consequence of the correlation being less than 1. Extreme values are partly due to the underlying tendency and partly due to random variation. The random component tends to be less extreme in subsequent measurements, producing regression toward the mean.

Regression toward the mean is one of the most commonly misinterpreted statistical phenomena. A sports team that performs exceptionally well in one season tends to perform worse the next season, not because of a decline in skill but because the exceptional performance was partly due to luck. A medical patient who is extremely sick at one visit tends to improve at the next, not because of treatment but because the extreme measurement was partly due to random fluctuation. Failing to account for regression toward the mean leads to spurious conclusions about the effectiveness of interventions.

Pearson and the correlation coefficient

Karl Pearson developed the correlation coefficient and the mathematical framework of regression in a series of papers beginning in 1896. Pearson generalised Galton's work from the bivariate normal distribution to arbitrary distributions and provided the formulas for the sample correlation coefficient and the regression line that are still used today. Pearson's correlation coefficient became the standard measure of association for quantitative variables and remains one of the most widely used statistics in science.

Pearson's contributions to regression extended beyond the correlation coefficient. He developed the method of moments for estimating parameters, introduced the chi-square test for goodness of fit, and established the first academic department of statistics at University College London. His textbook The Grammar of Science (1892) articulated a positivist philosophy in which correlation replaced causation as the fundamental concept, arguing that science should describe relationships rather than explain mechanisms.

Gauss and the method of least squares

The method of least squares was developed independently by Carl Friedrich Gauss and Adrien-Marie Legendre in the early nineteenth century. Legendre published first (1805), but Gauss claimed he had been using the method since 1795 for astronomical calculations. The priority dispute between Gauss and Legendre was one of many such disputes in the history of mathematics.

Gauss's contribution went beyond the method itself. He showed that the least squares estimator is the maximum likelihood estimator under normal errors, connecting the algebraic method to a probabilistic model. He also proved that the normal distribution is the unique error distribution for which the sample mean is the maximum likelihood estimator of the location parameter, providing a theoretical justification for both the normal distribution and the method of least squares.

The method of least squares was initially used for astronomical calculations: determining the orbits of comets and asteroids from noisy positional observations. The problem was to find the curve (orbit) that best fit a set of observed positions, where "best" meant minimising the sum of squared residuals. Legendre's 1805 publication Nouvelles methodes pour la determination des orbites des cometes introduced the method with the clear statement that "of all the principles that can be proposed for this purpose, there is none more general, more exact, or easier to apply, than that which consists of rendering the sum of the squares of the errors a minimum."

Gauss's 1809 work Theoria motus corporum coelestium showed that least squares is optimal under the assumption of normally distributed errors, and his 1823 work established the Gauss-Markov theorem (though in a weaker form than the modern version), showing that least squares is the best linear unbiased estimator even without the normality assumption. These results placed the method of least squares on a rigorous mathematical foundation.

Fisher and the development of regression inference

Fisher's 1922 paper "The Goodness of Fit of Regression Formulae, and Its Distribution" established the distributional theory for regression coefficients, including the t-distribution for testing the significance of the slope and the F-distribution for comparing nested models. Fisher also introduced the analysis of variance as a decomposition of the total sum of squares into components attributable to the regression and to the residual variation.

Fisher's work transformed regression from a method of curve fitting into a complete framework for statistical inference. The combination of point estimation (least squares), interval estimation (confidence intervals for coefficients), and hypothesis testing (t-tests and F-tests) provided the tools that scientists needed to draw conclusions from regression analyses.

The expansion to multiple regression and beyond

Multiple regression emerged naturally from the matrix formulation of least squares. The development of electronic computers in the 1950s and 1960s made it practical to invert large matrices, enabling the routine use of multiple regression with many predictors. The first statistical software packages (BMD, SPSS, SAS) included regression as a core procedure.

The development of regularisation methods in the 1970s and 1980s (ridge regression by Hoerl and Kennard in 1970, the lasso by Tibshirani in 1996) addressed the problems of multicollinearity and high-dimensional predictors. The machine learning era brought further innovations, including support vector regression, random forests, gradient boosting, and neural networks, all of which can be viewed as extensions of the basic regression idea to more flexible functional forms.

The philosophical significance of regression

Regression raises deep philosophical questions about the nature of prediction, explanation, and causation. The distinction between prediction (forecasting outcomes) and explanation (understanding mechanisms) maps onto the distinction between correlational and causal analysis. Regression can do both, but the requirements are different.

For prediction, all that matters is accuracy: does the model forecast well on new data? Regularised regression, cross-validation, and ensemble methods are designed to optimise prediction. For explanation, the model must correctly represent the underlying mechanism, which requires understanding the causal structure. The same model that predicts well may explain poorly if it captures spurious correlations rather than genuine causal relationships.

The tension between prediction and explanation is one of the defining issues of modern data science. Machine learning prioritises prediction; traditional statistics prioritises explanation. Both are legitimate goals, and the choice between them depends on the purpose of the analysis. A regression model used to guide medical treatment must be explanatory (the coefficients must represent causal effects). A regression model used to forecast sales need only be predictive.

Regression and the philosophy of science

The history of regression reflects the broader history of the philosophy of science. Galton's work was motivated by eugenics, a programme that confounded correlation with causation and assumed that biological inheritance determined social outcomes. Pearson's positivism replaced causal explanation with statistical description, arguing that science could only measure correlations, not discover causes. Fisher's randomisation framework provided the foundation for causal inference through experimental design, but his advocacy of randomised experiments implicitly conceded that observational data could not establish causation.

The modern potential outcomes framework (Rubin, 1974; Holland, 1986) makes the causal assumptions explicit and testable. It shows that regression can estimate causal effects, but only under conditions that are often not met in practice. The assumption of "no unmeasured confounders" is untestable from the data alone and requires substantive knowledge of the domain. This is why regression analysis in medicine, social science, and policy requires not just statistical expertise but also domain expertise.

Regression in the age of machine learning

Machine learning has both extended and challenged traditional regression analysis. On one hand, machine learning algorithms (decision trees, random forests, neural networks) can capture nonlinear relationships that linear regression misses. On the other hand, these algorithms often sacrifice interpretability for predictive accuracy. The resulting "black box" models may predict well but provide little insight into the mechanisms generating the data.

The explainable AI (XAI) movement seeks to recover interpretability by developing methods (SHAP values, LIME, counterfactual explanations) that approximate the behaviour of complex models with simpler, more interpretable ones. SHAP values, in particular, decompose the prediction of any model into additive contributions from each feature, generalising the coefficients of a linear regression to arbitrary models. The connection to regression is not coincidental: the goal of attributing predictions to features is fundamentally the same goal that motivated Galton and Pearson.

The mathematics of least squares

The method of least squares has a beautiful geometric interpretation. In matrix notation, the linear regression model is $\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{ϵ}$, where $\mathbf{y}$ is the $n\times 1$ response vector, $\mathbf{X}$ is the $n\times p$ design matrix, $\mathbf{\beta }$ is the $p\times 1$ coefficient vector, and $\mathbf{ϵ}$ is the error vector. The least squares estimator $\hat{\mathbf{\beta }}=({\mathbf{X}}^{\prime }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\prime }\mathbf{y}$ is the orthogonal projection of $\mathbf{y}$ onto the column space of $\mathbf{X}$.

The hat matrix $\mathbf{H}=\mathbf{X}({\mathbf{X}}^{\prime }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\prime }$ maps $\mathbf{y}$ to $\hat{\mathbf{y}}$. It is symmetric and idempotent (${\mathbf{H}}^2=\mathbf{H}$), properties shared by all projection matrices. The diagonal elements $h_{ii}$ of the hat matrix are the leverages: they measure how much influence observation $i$ has on its own fitted value. Observations with high leverage ($h_{ii}$ close to 1) can have a disproportionate influence on the regression line.

The Gauss-Markov theorem states that the least squares estimator is the best linear unbiased estimator (BLUE): among all unbiased estimators that are linear functions of $\mathbf{y}$, the least squares estimator has the smallest variance. This theorem does not require normality; it requires only that the errors have zero mean, constant variance, and zero correlation. The normality assumption is needed only for hypothesis tests and confidence intervals, not for point estimation.

Regression diagnostics and model checking

No regression analysis is complete without checking whether the model assumptions are satisfied. The four main assumptions are linearity (the relationship between predictors and response is linear), independence (the errors are independent), homoscedasticity (the errors have constant variance), and normality (the errors are normally distributed). These are conveniently remembered by the acronym LINE.

Residual plots are the primary diagnostic tool. A plot of residuals versus fitted values should show no pattern: a funnel shape indicates heteroscedasticity, a curved pattern indicates nonlinearity, and clusters indicate lack of independence. A normal Q-Q plot of the residuals should approximate a straight line: systematic deviations indicate non-normality. A scale-location plot (square root of absolute residuals versus fitted values) amplifies heteroscedasticity. A residuals versus leverage plot identifies influential observations.

Influential observations are those that substantially change the regression coefficients when removed. Cook's distance measures the overall influence of observation $i$ on all fitted values. An observation with $D_i>4/n$ or $D_i>1$ (conventions vary) is considered influential. DFBETAS measures the influence of observation $i$ on each individual coefficient. Identifying influential observations is not a license to delete them; it is a prompt to investigate whether they are genuine data points or recording errors.

Regression and the reproducibility crisis

The misuse of regression has contributed to the reproducibility crisis. Stepwise regression (automatically selecting predictors based on p-values) inflates effect sizes and produces models that do not replicate. The garden of forking paths (Gelman and Loken, 2014) is particularly acute in regression: there are many defensible choices of predictors, transformations, and model specifications, and the researcher may unconsciously choose the one that gives the strongest result.

The solution is pre-registration of the regression model (specifying predictors and functional form before seeing the data), validation on independent datasets, and emphasis on effect sizes and confidence intervals rather than p-values. These practices are becoming standard in fields that rely heavily on regression analysis.

Logistic regression and generalised linear models

Not all response variables are continuous. When the response is binary (success/failure, yes/no, diseased/healthy), linear regression is inappropriate because it can predict values outside the range [0, 1]. Logistic regression models the log-odds of the response as a linear function of the predictors: $\log (p/(1-p))={\beta }_0+{\beta }_1{x}_1+\dots +{\beta }_p{x}_p$.

The coefficients in logistic regression have an odds ratio interpretation: $e^{{\beta }_j}$ is the multiplicative change in the odds of success for a one-unit increase in $x_j$. Logistic regression is estimated by maximum likelihood, not by least squares, because the errors are not normally distributed. The deviance (analogous to the residual sum of squares) measures the goodness of fit.

Generalised linear models (GLMs) provide a unified framework for regression with different types of response variables. A GLM specifies a distribution for the response (normal, binomial, Poisson, gamma), a linear predictor $\eta ={\beta }_0+{\beta }_1{x}_1+\dots +{\beta }_p{x}_p$, and a link function $g$ connecting the mean of the response to the linear predictor: $g(\mu )=\eta$. Linear regression is a GLM with normal distribution and identity link. Logistic regression is a GLM with binomial distribution and logit link. Poisson regression (for count data) is a GLM with Poisson distribution and log link. GLMs are estimated by iteratively reweighted least squares (IRLS), which converges to the maximum likelihood estimates.

Regression and the problem of overfitting

Overfitting occurs when a regression model captures noise in the training data rather than the true underlying relationship. An overfitted model performs well on the data used to fit it but poorly on new data. The risk of overfitting increases with the number of predictors relative to the sample size: with $n$ observations and $p=n-1$ predictors, any linear regression will fit the training data perfectly ($R^2=1$), regardless of whether there is a genuine relationship.

The bias-variance trade-off captures the tension between underfitting (too few predictors, high bias) and overfitting (too many predictors, high variance). As the model complexity increases, bias decreases but variance increases. The optimal model complexity minimises the total prediction error, which is the sum of bias squared and variance. Regularisation methods (ridge, lasso) navigate this trade-off by shrinking the coefficients toward zero, reducing variance at the cost of some bias.

Cross-validation estimates the prediction error by partitioning the data into training and validation sets. K-fold cross-validation divides the data into $k$ equal parts, fits the model on $k-1$ parts, and evaluates on the held-out part, repeating for each part. The average prediction error across the $k$ folds estimates the out-of-sample prediction error. This estimate is used to select the optimal model complexity (number of predictors, regularisation parameter, polynomial degree).

Regularisation: ridge regression and the lasso

Ridge regression adds a penalty on the sum of squared coefficients to the least squares criterion: ${\min }_{\mathbf{\beta }}\sum (Y_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta }{)}^2+\lambda \sum {\beta }_j^2$. The penalty shrinks the coefficients toward zero, reducing their variance at the cost of introducing a small bias. Ridge regression is particularly effective when the predictors are highly correlated (multicollinearity), where OLS coefficients are unstable.

The lasso (Least Absolute Shrinkage and Selection Operator) uses an absolute value penalty: ${\min }_{\mathbf{\beta }}\sum (Y_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta }{)}^2+\lambda \sum |{\beta }_j|$. Unlike ridge, the lasso can set coefficients exactly to zero, performing automatic variable selection. The lasso produces sparse models (with fewer nonzero coefficients) that are easier to interpret.

The elastic net combines ridge and lasso penalties: ${\lambda }_1\sum |{\beta }_j|+{\lambda }_2\sum {\beta }_j^2$. This combination retains the lasso's variable selection property while handling groups of correlated predictors more effectively (tending to select or exclude the group together rather than picking one member arbitrarily). The regularisation parameter $\lambda$ is typically chosen by cross-validation.

Regression with time series data

Time series regression introduces additional challenges. When the response and predictor are both trending over time, the regression may capture the shared trend rather than a genuine relationship. This is the problem of spurious regression, identified by Granger and Newbold (1974): regressing one random walk on another produces spuriously significant coefficients with high $R^2$, even when the series are independent.

The solution is to check for stationarity (constant mean and variance over time) and to difference non-stationary series before regression. The Augmented Dickey-Fuller test checks for unit roots (a form of non-stationarity). If both series are stationary, standard regression methods apply. If both are non-stationary but cointegrated (a linear combination is stationary), an error correction model captures both the short-run dynamics and the long-run equilibrium relationship.

Autocorrelation in the residuals violates the independence assumption of standard regression. The Durbin-Watson test checks for first-order autocorrelation. If present, the standard errors must be corrected (using Newey-West heteroscedasticity and autocorrelation consistent standard errors) or the model must be extended to include lagged variables (autoregressive distributed lag models). These methods are standard in econometrics and financial statistics.

Bibliography Master

• Galton, F., "Regression Towards Mediocrity in Hereditary Stature," Journal of the Anthropological Institute 15 (1886), 246-263. Discovery of regression toward the mean and the regression line.

• Pearson, K., "Regression, Heredity, and Panmixia," Philosophical Transactions of the Royal Society A 187 (1896), 253-318. Mathematical development of the correlation coefficient and regression.

• Gauss, C. F., Theoria Motus Corporum Coelestium (Perthes and Besser, 1809). The method of least squares and its connection to the normal distribution.

• Fisher, R. A., "The Goodness of Fit of Regression Formulae, and Its Distribution," Journal of the Royal Statistical Society 85(4) (1922), 597-612. Distributional theory for regression coefficients.

• Nelder, J. A. and Wedderburn, R. W. M., "Generalized Linear Models," Journal of the Royal Statistical Society A 135(3) (1972), 370-384. Unified framework for regression with non-normal responses.

• Hoerl, A. E. and Kennard, R. W., "Ridge Regression: Biased Estimation for Nonorthogonal Problems," Technometrics 12(1) (1970), 55-67. Introduction of ridge regression.

• Tibshirani, R., "Regression Shrinkage and Selection via the Lasso," Journal of the Royal Statistical Society B 58(1) (1996), 267-288. Introduction of the lasso.

• Stigler, S. M., The History of Statistics: The Measurement of Uncertainty before 1900 (Harvard University Press, 1986). Chapters 8-10 cover the development of regression and correlation.

• Hastie, T., Tibshirani, R., and Friedman, J., The Elements of Statistical Learning (2e, Springer, 2009). Modern treatment of regression and its extensions in the context of statistical learning.