Lesson 12 · Statistics for data analysts

Multiple Regression & Model Diagnostics | Statistics #12

Video twelve of the 15-part series: multiple real predictors at once, and the diagnostics that catch when they're fighting each other. Real mtcars data,…

⬇ Download notebook ⬇ Notebook + data (.zip)Open in Colab ↗

CourseStatistics for data analysts
Lesson12 of 15
Video13 min
FormatJupyter notebook · 10 code cells
Data1 dataset

statsmodels pandas NumPy Matplotlib

What you'll learn

A Three-Predictor Model
Diagnosing Multicollinearity with VIF
Fixing It - Dropping the Redundant Predictor
Checking Residual Assumptions
Interpreting the Final Model
Wrap-Up: What You Learned

Datasets used in this lesson

Save these next to the notebook. In Google Colab, upload them with the 📁 icon on the left first.

mtcars.csv1.7 KB

📓 Full notebook

Download .ipynb

Statistics for Data Analysts, Video 12: Multiple Regression and Diagnostics#

Video twelve of the 15-part series: multiple real predictors at once, and the diagnostics that catch when they're fighting each other.
Real mtcars data, extended with statsmodels for a full regression summary and multicollinearity checks.
Let's get into it.

Before You Start#

Open a new Jupyter Notebook in VS Code and select your Python interpreter as the kernel.
Install statsmodels if you don't have it yet: pip install statsmodels.
Place mtcars.csv in the same folder as this notebook.

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan

mtcars = pd.read_csv('mtcars.csv')
print(mtcars[['model', 'mpg', 'wt', 'hp', 'disp']].head())

               model   mpg     wt   hp   disp
0          Mazda RX4  21.0  2.620  110  160.0
1      Mazda RX4 Wag  21.0  2.875  110  160.0
2         Datsun 710  22.8  2.320   93  108.0
3     Hornet 4 Drive  21.4  3.215  110  258.0
4  Hornet Sportabout  18.7  3.440  175  360.0

Part 1: A Three-Predictor Model#

model_full = smf.ols('mpg ~ wt + hp + disp', data=mtcars).fit()
print(model_full.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                    mpg   R-squared:                       0.827
Model:                            OLS   Adj. R-squared:                  0.808
Method:                 Least Squares   F-statistic:                     44.57
Date:                Sun, 16 Aug 2026   Prob (F-statistic):           8.65e-11
Time:                        19:30:09   Log-Likelihood:                -74.321
No. Observations:                  32   AIC:                             156.6
Df Residuals:                      28   BIC:                             162.5
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     37.1055      2.111     17.579      0.000      32.782      41.429
wt            -3.8009      1.066     -3.565      0.001      -5.985      -1.617
hp            -0.0312      0.011     -2.724      0.011      -0.055      -0.008
disp          -0.0009      0.010     -0.091      0.929      -0.022       0.020
==============================================================================
Omnibus:                        5.269   Durbin-Watson:                   1.367
Prob(Omnibus):                  0.072   Jarque-Bera (JB):                4.038
Skew:                           0.856   Prob(JB):                        0.133
Kurtosis:                       3.310   Cond. No.                     1.50e+03
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.5e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Real R-squared jumps to about 0.827, up from 0.753 with weight alone, so the fuller model explains more real variance. But look closer: the real coefficient on disp has a p-value of about 0.93, nowhere near significant, even though disp is plausibly related to mpg on its own. That's the first hint something is off, and the summary's own footnote agrees: 'condition number is large... this might indicate strong multicollinearity.'

Part 2: Diagnosing Multicollinearity with VIF#

X = mtcars[['wt', 'hp', 'disp']]
X_const = sm.add_constant(X)
vif_data = pd.DataFrame()
vif_data['feature'] = X_const.columns
vif_data['VIF'] = [variance_inflation_factor(X_const.values, i) for i in range(X_const.shape[1])]
print(vif_data)

  feature        VIF
0   const  20.473619
1      wt   4.844618
2      hp   2.736633
3    disp   7.324517

print(X.corr().round(3))

         wt     hp   disp
wt    1.000  0.659  0.888
hp    0.659  1.000  0.791
disp  0.888  0.791  1.000

Part 3: Fixing It - Dropping the Redundant Predictor#

model_reduced = smf.ols('mpg ~ wt + hp', data=mtcars).fit()
print(model_reduced.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                    mpg   R-squared:                       0.827
Model:                            OLS   Adj. R-squared:                  0.815
Method:                 Least Squares   F-statistic:                     69.21
Date:                Sun, 16 Aug 2026   Prob (F-statistic):           9.11e-12
Time:                        19:33:43   Log-Likelihood:                -74.326
No. Observations:                  32   AIC:                             154.7
Df Residuals:                      29   BIC:                             159.0
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     37.2273      1.599     23.285      0.000      33.957      40.497
wt            -3.8778      0.633     -6.129      0.000      -5.172      -2.584
hp            -0.0318      0.009     -3.519      0.001      -0.050      -0.013
==============================================================================
Omnibus:                        5.303   Durbin-Watson:                   1.362
Prob(Omnibus):                  0.071   Jarque-Bera (JB):                4.046
Skew:                           0.855   Prob(JB):                        0.132
Kurtosis:                       3.332   Cond. No.                         588.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

print(f'Full model:    R2={model_full.rsquared:.4f}, Adj R2={model_full.rsquared_adj:.4f}, AIC={model_full.aic:.2f}, BIC={model_full.bic:.2f}')
print(f'Reduced model: R2={model_reduced.rsquared:.4f}, Adj R2={model_reduced.rsquared_adj:.4f}, AIC={model_reduced.aic:.2f}, BIC={model_reduced.bic:.2f}')

Full model:    R2=0.8268, Adj R2=0.8083, AIC=156.64, BIC=162.51
Reduced model: R2=0.8268, Adj R2=0.8148, AIC=154.65, BIC=159.05

X2 = mtcars[['wt', 'hp']]
X2_const = sm.add_constant(X2)
vif2 = pd.DataFrame()
vif2['feature'] = X2_const.columns
vif2['VIF'] = [variance_inflation_factor(X2_const.values, i) for i in range(X2_const.shape[1])]
print(vif2)

  feature        VIF
0   const  12.161539
1      wt   1.766625
2      hp   1.766625

Part 4: Checking Residual Assumptions#

bp_stat, bp_pvalue, bp_fstat, bp_fpvalue = het_breuschpagan(model_reduced.resid, model_reduced.model.exog)
print(f'Real Breusch-Pagan statistic: {bp_stat:.4f}')
print(f'Real Breusch-Pagan p-value: {bp_pvalue:.4f}')

Real Breusch-Pagan statistic: 0.8807
Real Breusch-Pagan p-value: 0.6438

import matplotlib.pyplot as plt
plt.figure(figsize=(8, 4))
plt.scatter(model_reduced.fittedvalues, model_reduced.resid, color='steelblue')
plt.axhline(0, color='darkred', linestyle='--')
plt.title('Real Residuals vs. Fitted Values (Reduced Model)')
plt.xlabel('Fitted mpg')
plt.ylabel('Residual')
plt.show()

No description has been provided for this image

Part 5: Interpreting the Final Model#

wt_coef = model_reduced.params['wt']
hp_coef = model_reduced.params['hp']
print(f'Holding horsepower constant, each additional 1,000 lbs of real weight predicts {abs(wt_coef):.2f} fewer real mpg')
print(f'Holding weight constant, each additional real horsepower predicts {abs(hp_coef):.4f} fewer real mpg')
print(f'Real 100-horsepower difference, weight held constant: {abs(hp_coef) * 100:.2f} fewer real mpg')

Holding horsepower constant, each additional 1,000 lbs of real weight predicts 3.88 fewer real mpg
Holding weight constant, each additional real horsepower predicts 0.0318 fewer real mpg
Real 100-horsepower difference, weight held constant: 3.18 fewer real mpg

Wrap-Up: What You Learned#

Fitting a multiple regression with statsmodels and reading the full real summary table.
Diagnosing multicollinearity with VIF and a real correlation matrix, catching a predictor whose coefficient looked meaningless because it was redundant, not because it was irrelevant.
Comparing real models with adjusted R-squared, AIC, and BIC, which properly reward simplicity, not just raw fit.
The Breusch-Pagan test for constant residual variance, alongside the residual plot.
Interpreting multiple regression coefficients correctly, each one holding the others constant.
Video thirteen turns to categorical data specifically: a deeper chi-square treatment, odds ratios, and a first look at logistic regression. Subscribe so it lands automatically see you there.

Found this useful?

All lessons, notebooks and datasets here are free. If they helped you, a coffee keeps new lessons coming.

☕ Buy me a coffee Other ways to help