Diagnosing Regression Models with Residual Plots in Python

ByJordan Brown

Linear regression is one of the most widely used modeling techniques in data science and statistics. It‘s a powerful tool for understanding relationships between variables and making predictions. However, blindly trusting regression output without proper validation is risky. Even a model with high R-squared and significant p-values can be misleading if its assumptions are violated.

One of the key assumptions of ordinary least squares (OLS) regression is that the model residuals — the differences between the predicted and actual values — are well-behaved. Residual plots provide a visual way to check this assumption and diagnose potential issues with model fit. In this post, we‘ll dive deep into residual analysis for regression, using Python code examples throughout.

What are Residuals and Residual Plots?

In linear regression, we fit a line to the data that minimizes the sum of squared residuals. The residuals are the vertical distances between the observed data points and the predicted values on the line. Mathematically, the residual e_i for the i-th observation is:

e_i = y_i – ŷ_i

where y_i is the observed value and ŷ_i is the fitted value from the regression line.

If the linear model is a good fit for the data, we expect the residuals to be normally distributed with mean zero and constant variance at all fitted values. We also assume the residuals are independent from each other.

Residual plots let us visualize the residuals and check these assumptions. The basic residual plot is a scatter plot of residuals on the y-axis against the fitted values on the x-axis. We look for random scatter around the horizontal line at 0. Other types of residual plots test for normality, constant variance, outliers, and influential points. We‘ll see examples of these plots shortly.

Regression Assumptions and Residual Plots

Before we get to the Python code, let‘s review the key OLS assumptions and how residual plots help diagnose violations.

Linearity: The relationship between the predictors and the outcome variable is linear. If this is violated, the residuals will show a curved pattern when plotted against fitted values. Nonlinearity can be remedied with transformations or polynomial terms.

Homoscedasticity: The variance of the residuals is constant at all levels of the predictors. Heteroscedasticity (non-constant variance) shows up as a fan or cone shape in the residual plots. Possible fixes include weighted least squares or variance-stabilizing transformations.

Normality: The residuals are normally distributed at each level of the predictors. We check this with a normal Q-Q plot of the residuals. Non-normality can appear as an S-shaped curve or heavy tails. Transformations can sometimes correct this, or we can use robust regression methods.

Independence: The residuals are independent of each other, without correlation or autocorrelation. This is harder to see in residual plots, but patterns or clustering can indicate dependence. Remedies include time series methods, mixed models, or generalized least squares.

Residual plots can also identify outliers (isolated points far from the regression line) and high leverage points (extreme x-values that pull the line towards them). These unusual points aren‘t necessarily bad, but they warrant scrutiny to make sure they aren‘t data errors or overly influential.

Creating Residual Plots in Python

Now let‘s see how to create and interpret diagnostic plots for regression in Python. We‘ll use a simple example with one predictor and one outcome variable, but these techniques generalize to multiple regression. We‘ll use the seaborn library for plotting, along with numpy and statsmodels.

First, let‘s generate some data and fit a regression line:

import numpy as np
import statsmodels.api as sm
import seaborn as sns

# Generate random data
x = np.random.normal(10, 1, 50)  
y = 2 + 0.5*x + np.random.normal(0, 1, 50)

# Fit regression model
X = sm.add_constant(x)
model = sm.OLS(y, X).fit()

We can get the fitted values and residuals from the model object:

fitted = model.fittedvalues
residuals = model.resid

Residuals vs Fitted Plot

The basic diagnostic plot is a scatter plot of residuals vs fitted values. We look for random scatter around the horizontal line at 0:

sns.residplot(x=fitted, y=residuals, lowess=True, line_kws={‘color‘: ‘red‘}) 
plt.title(‘Residuals vs Fitted‘)
plt.xlabel(‘Fitted values‘)
plt.ylabel(‘Residuals‘)

Residuals vs Fitted Plot

The red line is a LOWESS smooth that helps see patterns. Here, the residuals look randomly scattered, which is good. But what if we simulated some data with a nonlinear relationship?

x = np.random.normal(10, 1, 50)
y_nonlinear = 10 + 0.2*x + 0.05*x**2 + np.random.normal(0, 1, 50)

X_nonlinear = sm.add_constant(x) 
model_nonlinear = sm.OLS(y_nonlinear, X_nonlinear).fit()

sns.residplot(x=model_nonlinear.fittedvalues, y=model_nonlinear.resid, lowess=True, line_kws={‘color‘: ‘red‘})

Residuals vs Fitted Nonlinear

The U-shaped pattern indicates the linear model doesn‘t capture the quadratic relationship. We could fix this by adding a squared term to the model.

Normal Q-Q Plot

To check normality of residuals, we use a Q-Q (quantile-quantile) plot. This plots the quantiles of the residual distribution against the quantiles of a normal distribution. If the residuals are normal, the points will fall along the diagonal line. Here‘s how to create a Q-Q plot in seaborn:

sns.set_style(‘whitegrid‘)
fig, ax = plt.subplots()
fig = sm.qqplot(model.resid, fit=True, line=‘45‘, ax=ax)

Normal Q-Q Plot

The residuals fall along the line, indicating they are approximately normally distributed. Let‘s see what happens if we generate non-normal residuals:

x = np.random.normal(10, 1, 50)
residuals_nonnormal = np.random.exponential(1, 50)
y_nonnormal = 5 + 0.5*x + residuals_nonnormal

model_nonnormal = sm.OLS(y_nonnormal, sm.add_constant(x)).fit()
fig = sm.qqplot(model_nonnormal.resid, fit=True, line=‘45‘)

Non-normal Q-Q Plot

The curved pattern indicates the residuals have a right-skewed, non-normal distribution. We might consider a log transformation to make the residuals more normal.

Scale-Location Plot

Also known as the spread-location plot, this shows if the residual variance is constant at all levels of the predicted values (homoscedasticity). It‘s a scatter plot of the square root of the absolute value of the standardized residuals against the fitted values. We look for a horizontal line with no trend:

fitted = model.fittedvalues
resid_standardized = model.get_influence().resid_studentized_internal

sns.regplot(x=fitted, y=np.sqrt(np.abs(resid_standardized)), 
            ci=None, lowess=True, line_kws={‘color‘: ‘red‘})
plt.title(‘Scale-Location‘)
plt.xlabel(‘Fitted values‘)
plt.ylabel(r‘$\sqrt{|Standardized Residuals|}$‘);

Scale-Location Plot

The flat red line shows the residual variance is roughly constant at all fitted values. Now let‘s generate heteroscedastic data and see how this plot looks:

x_het = np.random.normal(10, 1, 50)
y_het = 5 + 0.5*x_het +  np.random.normal(0, 0.5 + 0.1*x_het, 50)

model_het = sm.OLS(y_het, sm.add_constant(x_het)).fit()
resid_standardized_het = model_het.get_influence().resid_studentized_internal

sns.regplot(x=model_het.fittedvalues, y=np.sqrt(np.abs(resid_standardized_het)), 
            ci=None, lowess=True, line_kws={‘color‘: ‘red‘})

Heteroscedastic Scale-Location Plot

The upward trend in this plot signals heteroscedasticity – the residual variance increases with the level of the fitted values. Weighted least squares or robust standard errors might be needed here.

Residuals vs Leverage Plot

Finally, let‘s check for influential outliers with residuals vs leverage. Leverage measures how far away the predictor values are from their mean. High leverage points can have a strong effect on the regression line. We identify high leverage points and large residuals (potential outliers) with reference lines:

from statsmodels.stats.outliers_influence import OLSInfluence

fig, ax = plt.subplots(figsize=(8,6))
fig = sm.graphics.influence_plot(model, alpha=0.05, ax=ax, criterion="cooks")

Influence Plot

This plot shows the standardized residuals against leverage. Dashed lines indicate values that are above a threshold for high leverage or large standardized residuals. The index labels on the points show which observations exceed the thresholds. No points are flagged here.

Let‘s add an influential outlier to the data and see how this plot changes:

x_out = np.append(x, 20)
y_out = np.append(y, 0)

model_out = sm.OLS(y_out, sm.add_constant(x_out)).fit()
fig = sm.graphics.influence_plot(model_out, alpha=0.05, criterion="cooks")

Influence Plot with Outlier

Observation 50 is flagged as both high leverage (far from the mean on the x-axis) and a large residual. If we remove this point and refit the model, the regression line would change noticeably.

Dealing with Problematic Residual Plots

If our residual plots show problems, what can we do? It depends on the issue, but some common tactics include:

  • Nonlinearity: Try polynomial terms, splines, or other nonlinear transformations of the predictors. Plot the residuals against these transformed predictors to check for improvement.

  • Heteroscedasticity: Use weighted least squares, which gives less weight to observations with large residual variance. Or use heteroscedasticity-robust standard errors for inference.

  • Non-normal residuals: Apply variance-stabilizing transformations like log or square root to the outcome variable. For heavy-tailed distributions, use robust regression methods like Huber or bisquare weights.

  • Correlated residuals: For time series data, use ARIMA models or Newey-West standard errors. For clustered data, use mixed models with random effects.

  • High leverage points: Examine them to see if they are data errors or outliers that don‘t fit the model. Consider removing them if they have undue influence on the results.

It‘s important to remember that residual plots alone don‘t prove a model is "correct". They only show if the model is a reasonable fit for the data and its assumptions are not grossly violated. We should also think about the practical and scientific context.

Conclusion

Residual plots are a key tool for diagnosing linear regression models. They help us check assumptions, spot problems, and identify ways to improve model fit. Using Python libraries like seaborn and statsmodels, it‘s easy to create and customize these diagnostic plots.

The examples we‘ve covered show some common patterns to look for: random scatter of residuals around 0, normally distributed residuals, constant variance, and no overly influential points. But we‘ve also seen how real data often violates these ideal conditions. The residual plots let us see these problems and consider solutions.

Still, residual diagnostics have limitations. They don‘t work well for small sample sizes, and they can miss subtle patterns. We should combine them with numerical tests and cross-validation to get a full picture of model performance.

I hope this deep dive into residual plots has been helpful! For more on regression diagnostics and Python data science, check out the references below. As always, the key is to think critically, visualize your data, and iterate to find the best model for your problem.

References

How useful was this post?

Click on a star to rate it!

Average rating 0 / 5. Vote count: 0

No votes so far! Be the first to rate this post.

You May Like to Read,

Similar Posts