A Comprehensive Guide to OLS Regression: Part 1
Introduction
Ordinary least squares (OLS) regression is the workhorse of statistics and machine learning. It‘s a simple yet powerful technique for modeling linear relationships between variables. OLS has an illustrious history dating back over 200 years, and remains widely used to this day.
In this guide, we‘ll take a deep dive into the fundamentals of OLS regression. I‘ll explain what OLS is, the math behind how it works, and how to implement it in Python. I‘ll also discuss diagnostics for checking the assumptions of OLS and what to do when they‘re violated.
By the end of this guide, you‘ll have a solid understanding of OLS regression and be well-equipped to apply it to your own data. Let‘s get started!
What is OLS Regression?
OLS regression is a method for estimating the unknown parameters in a linear regression model. In a linear regression, we have a dependent variable y that we want to predict based on one or more independent variables (also called predictors, features, or covariates) X. The goal is to find the linear function of X that best fits the observed data.
Mathematically, a linear regression model takes the form:
y = β0 + β1×1 + β2×2 + … + βkxk + ε
where:
- y is the dependent variable
- x1, x2, …, xk are the independent variables
- β0, β1, …, βk are the unknown parameters to be estimated
- ε is the error term representing the difference between the observed and predicted y values
The parameters β0, β1, …, βk define a hyperplane in the feature space that best fits the data. β0 is the intercept (the predicted value of y when all x = 0) and β1, …, βk are the slopes (the expected change in y for a 1 unit increase in the corresponding x).
A Brief History of OLS
The method of least squares was independently developed by Carl Friedrich Gauss and Adrien-Marie Legendre in the early 1800s. Gauss used it to fit curves to astronomical data, while Legendre applied it to geodesy (the measurement of the Earth‘s shape).
In 1821, Gauss published a detailed account of the method in his monumental work Theoria Combinationis Observationum Erroribus Minimis Obnoxiae (Theory of the Combination of Observations Least Subject to Error). He showed that least squares provides the "best" linear unbiased estimates of the unknown parameters, where best means minimum variance.
Least squares became widely used in the early 1900s with the rise of modern statistics. Many of the foundational results in regression analysis and analysis of variance (ANOVA) are based on least squares.
Today, OLS remains the most commonly used type of regression, both in academia and industry. Nearly all statistical software packages include functions for fitting OLS models.
The Math Behind OLS
Now let‘s dive into the mathematical details of how OLS works. The key idea is to choose the parameter values that minimize the sum of squared residuals (SSR). The residuals are the differences between the observed y values and the predicted y values from the model:
ei = yi – ŷi
where yi is the ith observed value and ŷi is the ith predicted value. The sum of squared residuals is then:
SSR = Σ(yi – ŷi)^2 = e1^2 + e2^2 + … + en^2
Minimizing the SSR leads to the least squares parameter estimates. This can be done using calculus by taking the partial derivatives of the SSR with respect to each parameter and setting them equal to zero:
∂SSR/∂β0 = -2 Σ(yi – β0 – β1xi1 – … – βkxik) = 0
∂SSR/∂β1 = -2 Σ(yi – β0 – β1xi1 – … – βkxik) xi1 = 0
…
∂SSR/∂βk = -2 Σ(yi – β0 – β1xi1 – … – βkxik) * xik = 0
This results in a system of k+1 linear equations called the normal equations. In matrix notation, the normal equations can be written as:
(X^T X) β = X^T y
where X is the n x (k+1) design matrix (a matrix with a column of 1‘s for the intercept plus the independent variable values), X^T is its transpose, and y is the vector of observed values. This system can be solved for the least squares estimates:
β = (X^T X)^(-1) X^T y
With modern computing power, the normal equations can be solved very efficiently even for large datasets with many variables. Most statistical software uses the QR decomposition rather than the normal equations for enhanced numerical stability.
OLS in Python
Let‘s see how to fit OLS regression models in Python. We‘ll use the Boston Housing dataset which contains median house values and various predictor variables for 506 neighborhoods around Boston.
First we load the data from scikit-learn and create the design matrix X and response vector y:
from sklearn.datasets import load_boston
boston = load_boston()
X = boston.data
y = boston.targetThe boston dataset contains the following variables:
| Variable | Description |
|---|---|
| CRIM | per capita crime rate by town |
| ZN | proportion of residential land zoned for lots over 25,000 sq.ft. |
| INDUS | proportion of non-retail business acres per town |
| CHAS | Charles River dummy variable (= 1 if tract bounds river; 0 otherwise) |
| NOX | nitric oxides concentration (parts per 10 million) |
| RM | average number of rooms per dwelling |
| AGE | proportion of owner-occupied units built prior to 1940 |
| DIS | weighted distances to five Boston employment centres |
| RAD | index of accessibility to radial highways |
| TAX | full-value property-tax rate per $10,000 |
| PTRATIO | pupil-teacher ratio by town |
| B | 1000(Bk – 0.63)^2 where Bk is the proportion of blacks by town |
| LSTAT | % lower status of the population |
| MEDV | Median value of owner-occupied homes in $1000‘s |
Next we fit the model using statsmodels:
import statsmodels.api as sm
model = sm.OLS(y, sm.add_constant(X))
results = model.fit()
print(results.summary())This prints a summary of the fitted model:
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.741
Model: OLS Adj. R-squared: 0.734
Method: Least Squares F-statistic: 108.1
Date: Fri, 09 Jun 2023 Prob (F-statistic): 2.16e-135
Time: 10:41:53 Log-Likelihood: -1498.8
No. Observations: 506 AIC: 3026.
Df Residuals: 492 BIC: 3085.
Df Model: 13
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 36.4595 5.103 7.144 0.000 26.432 46.487
x1 -0.1080 0.033 -3.287 0.001 -0.173 -0.043
x2 0.0464 0.014 3.380 0.001 0.019 0.073
x3 2.7187 0.126 21.587 0.000 2.471 2.966
x4 -17.7958 3.820 -4.659 0.000 -25.302 -10.289
x5 3.8048 0.417 9.116 0.000 2.984 4.625
x6 0.0007 0.002 0.354 0.723 -0.003 0.004
x7 -1.4756 0.199 -7.398 0.000 -1.868 -1.083
x8 0.3057 0.061 5.019 0.000 0.186 0.425
x9 -0.0123 0.003 -3.661 0.000 -0.019 -0.006
x10 -0.9527 0.130 -7.308 0.000 -1.210 -0.696
x11 0.0094 0.004 2.203 0.028 0.001 0.018
x12 -0.5248 0.051 -10.267 0.000 -0.625 -0.424
x13 4.6334 0.151 30.641 0.000 4.336 4.931
==============================================================================
Omnibus: 178.029 Durbin-Watson: 1.078
Prob(Omnibus): 0.000 Jarque-Bera (JB): 783.126
Skew: 1.521 Prob(JB): 8.84e-171
Kurtosis: 8.281 Cond. No. 1.51e+04
==============================================================================The R-squared of 0.741 indicates that the model explains about 74% of the variance in median home values. The coefficients give the estimated effect of each predictor on the median home value, holding the other predictors constant. For example, an increase of one room (RM) is associated with an increase of $3,805 in median value, on average.
Model Diagnostics
After fitting the model, it‘s important to check the assumptions of OLS using diagnostic plots.
The residuals vs fitted plot checks the assumptions of linearity and homoscedasticity. The residuals should be randomly scattered around 0 with no systematic patterns.
import seaborn as sns
sns.residplot(x=results.fittedvalues, y=results.resid, lowess=True,
line_kws={‘color‘: ‘red‘, ‘lw‘: 1, ‘alpha‘: 1})
Here we see some curvature and evidence of heteroscedasticity (non-constant variance). This suggests that a linear model may not be appropriate and a nonlinear transformation (like log) may help.
The Q-Q plot checks the assumption of normality of the residuals. If the residuals are normally distributed, the points should fall along a straight diagonal line.
import scipy.stats as stats
stats.probplot(results.resid, dist="norm", plot=sns.mpl.pyplot)
The tails falling off the line indicate that the residuals are not quite normally distributed. Outliers or skewness in the predictors or response can cause this.
If the diagnostics reveal problems with the assumptions, some remedies include:
- Transforming the response or predictors
- Adding interaction terms or polynomials
- Using robust standard errors
- Fitting a generalized linear model
- Using a nonlinear regression method
Comparing OLS to Other Regression Methods
OLS is just one type of regression. It‘s the most basic and widely used, but there are many other techniques that can be used depending on the structure of the data and the goals of the analysis.
Some common linear regression variants include:
- Ridge Regression: Adds L2 regularization to handle multicollinearity
- Lasso Regression: Uses L1 regularization for feature selection
- Elastic Net: Combines L1 and L2 penalties
- Stepwise Regression: Iteratively adds or removes predictors
- Partial Least Squares: Projects predictors onto a smaller set of components
- Generalized Least Squares: Accounts for heteroscedasticity and correlation in errors
For nonlinear relationships, some options are:
- Polynomial Regression: Fits a polynomial function of the predictors
- Spline Regression: Uses piecewise polynomials with constraints for smoothness
- Local Regression (LOESS): Fits a low-degree polynomial at each point
- Generalized Additive Models: Build a model by summing smooth functions of each predictor
There are also regression methods for non-continuous response variables:
- Logistic Regression: For binary outcomes
- Multinomial Regression: For unordered categorical outcomes
- Ordinal Regression: For ordered categorical outcomes
- Poisson Regression: For count outcomes
- Negative Binomial Regression: For overdispersed count data
In machine learning, there are even more types of regression methods like:
- Decision Trees and Random Forests
- Support Vector Machines
- Neural Networks
- Gradient Boosting Machines
- Gaussian Process Regression
Each of these methods has its own strengths and weaknesses. The choice of which to use depends on factors like the size and structure of the data, the presence of nonlinearity and interaction effects, the need for interpretability vs predictive accuracy, and computational constraints.
Conclusion
In this guide, we took a comprehensive look at ordinary least squares regression. We covered what it is, the math behind how it works, and how to implement it in Python. We also discussed model diagnostics for checking the assumptions of OLS and some potential remedies when they‘re violated.
OLS is a powerful and widely used tool for understanding relationships between variables. However, it‘s important to remember that it‘s just one type of regression and there are many other techniques that may be more appropriate depending on the data and the goals of the analysis.
Proper use of OLS requires carefully checking the assumptions and trying different transformations and models to ensure reliable results. Blindly applying OLS without verification is a common pitfall that can lead to misleading conclusions.
In future guides, we‘ll dive into some of the more advanced regression techniques and discuss how to select and tune models for optimal performance. Stay tuned!
