Lesson 24 · Probability and Statistics in python

Understanding Multiple Linear Regression with Python’s Statsmodels Library

Welcome! In this lesson, you will learn the basics of probability, descriptive statistics, and multiple linear regression using Python. By the end, you will…

⬇ Download notebook Open in Colab ↗

CourseProbability and Statistics in python
Lesson24 of 35
Video11 min
FormatJupyter notebook · 14 code cells

statsmodels NumPy seaborn pandas Matplotlib

What you'll learn

Probability: Flipping a Coin
Descriptive Statistics: Summarize Your Data
Probability Distributions
Exploring Relationships: Correlation
Linear Regression: Predicting MPG From Just Horsepower
Multiple Linear Regression: Predicting MPG from Several Features
Best Practices and Troubleshooting
Challenge Exercise

Data

No separate download needed — the notebook creates or downloads everything it uses.

📓 Full notebook

Download .ipynb

Multiple Linear Regression in Python: Absolute Beginner Lesson#

Welcome! In this lesson, you will learn the basics of probability, descriptive statistics, and multiple linear regression using Python.

By the end, you will understand how to apply regression to real-world data and explore relationships between variables.

Let us get started!

import warnings; warnings.filterwarnings("ignore")
# All warnings will be suppressed for a clean output

Probability: Flipping a Coin#

Let's start with a simple probability example using a coin flip.

Probability tells us how likely something is to happen.

import numpy as np

# Simulate 10 coin flips (0=heads, 1=tails)
coin_flips = np.random.choice([0,1], size=10)
print("Coin flips:", coin_flips)

# Calculate the probability of getting heads
prob_heads = np.mean(coin_flips == 0)
print("Probability of heads:", prob_heads)

Coin flips: [1 0 0 1 1 0 0 1 1 1]
Probability of heads: 0.4

Descriptive Statistics: Summarize Your Data#

Descriptive statistics help us quickly describe data using numbers like mean, median, and standard deviation.

Let us work with some fun real data next!

# Data setup: Load the MPG (Miles per Gallon) dataset
import seaborn as sns
import pandas as pd

mpg = sns.load_dataset("mpg")
print("MPG data shape:", mpg.shape)
mpg.head()

MPG data shape: (398, 9)

	mpg	cylinders	displacement	horsepower	weight	acceleration	model_year	origin	name
0	18.0	8	307.0	130.0	3504	12.0	70	usa	chevrolet chevelle malibu
1	15.0	8	350.0	165.0	3693	11.5	70	usa	buick skylark 320
2	18.0	8	318.0	150.0	3436	11.0	70	usa	plymouth satellite
3	16.0	8	304.0	150.0	3433	12.0	70	usa	amc rebel sst
4	17.0	8	302.0	140.0	3449	10.5	70	usa	ford torino

# Quickly look at summary statistics for mpg
mpg_describe = mpg['mpg'].describe()
print(mpg_describe)

# Find the median explicitly
median_mpg = mpg['mpg'].median()
print("Median MPG:", median_mpg)

count    398.000000
mean      23.514573
std        7.815984
min        9.000000
25%       17.500000
50%       23.000000
75%       29.000000
max       46.600000
Name: mpg, dtype: float64
Median MPG: 23.0

# Check for missing data in all columns
missing = mpg.isnull().sum()
print("Missing data by column:\n", missing)

# Remove any rows with missing mpg or horsepower values
mpg_clean = mpg.dropna(subset=['mpg', 'horsepower'])
print("Data shape after dropna:", mpg_clean.shape)

Missing data by column:
 mpg             0
cylinders       0
displacement    0
horsepower      6
weight          0
acceleration    0
model_year      0
origin          0
name            0
dtype: int64
Data shape after dropna: (392, 9)

Probability Distributions#

Probability distributions help us describe how values are spread out.

Let us simulate rolling a die using a discrete distribution!

# Simulate 1000 dice rolls (1 through 6)
rolls = np.random.choice([1,2,3,4,5,6], size=1000)

# Probability of getting a 6
prob_6 = np.mean(rolls == 6)
print("Probability of 6:", prob_6)

Probability of 6: 0.177

# Visualize the distribution of MPG values
import matplotlib.pyplot as plt

plt.figure(figsize=(6, 3))
plt.hist(mpg_clean['mpg'], bins=20, alpha=0.7, color='skyblue')
plt.title('Distribution of Miles Per Gallon')
plt.xlabel('Miles Per Gallon')
plt.ylabel('Count')
plt.show()

No description has been provided for this image

Exploring Relationships: Correlation#

Correlation measures how two variables move together.

Values close to 1 or -1 mean strong relationships. Values close to 0 mean weak or no relationship.

Let us check how mpg relates to horsepower.

# Calculate correlation between mpg and horsepower
corr = mpg_clean['mpg'].corr(mpg_clean['horsepower'])
print("Correlation (mpg vs horsepower):", corr)

Correlation (mpg vs horsepower): -0.7784267838977762

Linear Regression: Predicting MPG From Just Horsepower#

Regression predicts the value of one variable from one or more other variables.

Let us start by predicting mpg from horsepower using only one variable.

# Fit a simple linear regression (mpg ~ horsepower) with statsmodels
import statsmodels.api as sm

X_simple = mpg_clean[['horsepower']]
X_simple = sm.add_constant(X_simple)
y = mpg_clean['mpg']

model_simple = sm.OLS(y, X_simple).fit()
print(model_simple.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                    mpg   R-squared:                       0.606
Model:                            OLS   Adj. R-squared:                  0.605
Method:                 Least Squares   F-statistic:                     599.7
Date:                Wed, 03 Sep 2025   Prob (F-statistic):           7.03e-81
Time:                        08:48:20   Log-Likelihood:                -1178.7
No. Observations:                 392   AIC:                             2361.
Df Residuals:                     390   BIC:                             2369.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         39.9359      0.717     55.660      0.000      38.525      41.347
horsepower    -0.1578      0.006    -24.489      0.000      -0.171      -0.145
==============================================================================
Omnibus:                       16.432   Durbin-Watson:                   0.920
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               17.305
Skew:                           0.492   Prob(JB):                     0.000175
Kurtosis:                       3.299   Cond. No.                         322.
==============================================================================

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

Multiple Linear Regression: Predicting MPG from Several Features#

Multiple linear regression lets us use several predictors at once.

We can use horsepower, weight, and cylinders to better predict mpg.

Let us do it!

# Build the model with more predictors
X_multi = mpg_clean[['horsepower', 'weight', 'cylinders']]
X_multi = sm.add_constant(X_multi)
model_multi = sm.OLS(y, X_multi).fit()
print(model_multi.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                    mpg   R-squared:                       0.708
Model:                            OLS   Adj. R-squared:                  0.705
Method:                 Least Squares   F-statistic:                     313.1
Date:                Wed, 03 Sep 2025   Prob (F-statistic):          3.22e-103
Time:                        08:49:14   Log-Likelihood:                -1120.1
No. Observations:                 392   AIC:                             2248.
Df Residuals:                     388   BIC:                             2264.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         45.7368      0.796     57.461      0.000      44.172      47.302
horsepower    -0.0427      0.012     -3.677      0.000      -0.066      -0.020
weight        -0.0053      0.001     -8.208      0.000      -0.007      -0.004
cylinders     -0.3890      0.299     -1.302      0.194      -0.977       0.199
==============================================================================
Omnibus:                       37.624   Durbin-Watson:                   0.862
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               50.959
Skew:                           0.697   Prob(JB):                     8.60e-12
Kurtosis:                       4.085   Cond. No.                     1.15e+04
==============================================================================

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

# Visualize predictions vs. actual MPG
pred_mpg = model_multi.predict(X_multi)
plt.figure(figsize=(5, 5))
plt.scatter(y, pred_mpg, alpha=0.7, color='teal')
plt.xlabel('Actual MPG')
plt.ylabel('Predicted MPG')
plt.title('Actual vs Predicted MPG (Multiple Regression)')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--')
plt.show()

# Check for multicollinearity (predictors that are too related to each other)
from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame()
vif_data['feature'] = X_multi.columns
vif_data['VIF'] = [variance_inflation_factor(X_multi.values, i) for i in range(X_multi.shape[1])]
print(vif_data)

      feature        VIF
0       const  13.837987
1  horsepower   4.358007
2      weight   6.485732
3   cylinders   5.660847

# Use the model to predict mpg for a new car
print("Enter new horsepower:")
hp = float(input())
print("Enter new weight:")
wt = float(input())
print("Enter number of cylinders:")
cy = float(input())

new_X = pd.DataFrame({'const': [1], 'horsepower': [hp], 'weight': [wt], 'cylinders': [cy]})
predicted_mpg = model_multi.predict(new_X)[0]
print(f"Predicted MPG for your car: {predicted_mpg:.2f}")

Enter new horsepower:
Enter new weight:
Enter number of cylinders:
Predicted MPG for your car: 26.73

Best Practices and Troubleshooting#

Always check for missing data and outliers.
Keep an eye on multicollinearity using VIF.
Try different models and be sure to visualize your results.

If your predictions do not make sense, double-check your data and calculations.

# Try making your model with only two predictors
X_two = mpg_clean[['weight', 'cylinders']]
X_two = sm.add_constant(X_two)
model_two = sm.OLS(y, X_two).fit()
print(model_two.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                    mpg   R-squared:                       0.697
Model:                            OLS   Adj. R-squared:                  0.696
Method:                 Least Squares   F-statistic:                     448.4
Date:                Wed, 03 Sep 2025   Prob (F-statistic):          1.03e-101
Time:                        08:52:03   Log-Likelihood:                -1126.9
No. Observations:                 392   AIC:                             2260.
Df Residuals:                     389   BIC:                             2272.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         46.2923      0.794     58.305      0.000      44.731      47.853
weight        -0.0063      0.001    -10.922      0.000      -0.007      -0.005
cylinders     -0.7214      0.289     -2.493      0.013      -1.290      -0.152
==============================================================================
Omnibus:                       45.369   Durbin-Watson:                   0.831
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               69.921
Skew:                           0.748   Prob(JB):                     6.56e-16
Kurtosis:                       4.429   Cond. No.                     1.13e+04
==============================================================================

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

Challenge Exercise#

Build a regression model to predict mpg using 'horsepower' and 'acceleration'.

How does the R-squared value compare to your model with three predictors?

Try visualizing the predictions too!

Recap: What You Learned#

How to load and clean real datasets in Python
How to calculate descriptive statistics, visualize distributions, and check relationships
How to build and interpret multiple linear regression models
How to predict new values using regression

Keep practicing and you will master these skills!

Thanks for Learning Multiple Regression in Python!#

If you found this helpful, like and subscribe to the channel for more beginner data science lessons!

See you next time.

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