Python correlation and regression measure how two variables move together: Pearson and Spearman correlation, correlation matrices, simple linear regression from scratch and with SciPy, R-squared, and residual analysis, every number run on real code.

“Correlation does not imply causation, but it does waggle its eyebrows suggestively.”

Randall Munroe, xkcd

Correlation answers one simple question: when one number goes up, does the other one tend to go up too, go down, or just do its own thing? Regression goes one step further and draws the best straight line through the cloud of points, so you can actually predict one value from the other. That is the whole idea. Everything below is just turning that idea into Python you can run.

Here is a quick way to feel it. Picture a kid’s height and their shoe size. Taller kids usually have bigger feet, so the two move together: that is positive correlation. Now picture the temperature outside and how many layers of clothing people wear. Hotter day, fewer layers: that is negative correlation. Correlation runs from -1 (a perfect downhill line) to +1 (a perfect uphill line), with 0 meaning the two numbers basically ignore each other.

Last Updated: July 2026 | Tested on: Python 3.14.6, NumPy 2.4.6, pandas 2.3.3, SciPy 1.18.0 | Difficulty: Advanced | Reading Time: 19 minutes

A short story shows why both tools matter. Vinay, a 27-year-old analyst, found a strong correlation (r = 0.85) between advertising spend and revenue. His manager wanted to double the ad budget that same afternoon. Vinay ran a proper regression first and saw the relationship flatten out after a certain spend level, so doubling the budget would lift revenue by about 15%, not 100%. Correlation told him the link was real. Regression told him the shape of that link, and saved the company a pile of money.

One warning to keep taped to your monitor: a strong correlation never proves that one thing causes the other. Ice cream sales and sunburns rise together, but ice cream does not cause sunburn. The sun does both. Hold that thought, because it is the single most common mistake people make with this whole topic.

The diagram maps the whole Python correlation and regression flow you will follow. You start on the left by measuring relationship strength with a correlation coefficient: Pearson for straight-line links, Spearman for ranked or monotonic links, Kendall for small ordinal samples. Then you move right to regression, where you actually model the relationship and predict one variable from another. Correlation asks “how strongly are these two related?” Regression asks “given x, what is my best guess for y?” Keep one thing in mind as you read across: a strong correlation never means one variable causes the other.

Prerequisites

You do not need heavy math for this post. School-level algebra is enough: if you can read y = mx + b and remember that it draws a straight line, you are ready. It also helps to have finished the hypothesis testing tutorial (so p-values are not new) and the NumPy linear algebra tutorial (so arrays and the mean of a column feel familiar). On the tools side, install the data-science stack: pip install numpy pandas scipy. The code here was run on NumPy 2.4.6, pandas 2.3.3, and SciPy 1.18.0.

Python Correlation: Measuring Relationships

Let us start with the most common Python correlation measure: Pearson r. It measures how close your points sit to a straight line. Its cousin Spearman rho does the same job on the ranks of the values instead of the raw values, so it still catches a relationship even when the line bends, and it shrugs off a few wild outliers. Think of Pearson as a strict ruler that only likes straight lines, and Spearman as a more relaxed friend who just asks “do they generally go up together?”

The example below makes up 100 students. More study hours should push exam scores up, and more study hours should leave a little less time for sleep. We then ask SciPy to measure both relationships and print a full correlation matrix with pandas.

📄 correlation.py: Pearson and Spearman correlation

import numpy as np
import pandas as pd
from scipy import stats

rng = np.random.default_rng(42)

# Generate correlated data
n = 100
study_hours = rng.uniform(1, 10, n)
exam_score = 50 + 4 * study_hours + rng.normal(0, 5, n)

# Pearson correlation (linear relationships)
r, p = stats.pearsonr(study_hours, exam_score)
print(f"Pearson r: {r:.4f}, p-value: {p:.4e}")
print(f"Interpretation: {'Strong' if abs(r) > 0.7 else 'Moderate' if abs(r) > 0.4 else 'Weak'} {'positive' if r > 0 else 'negative'} correlation")

# Spearman correlation (monotonic relationships, robust to outliers)
rho, p_s = stats.spearmanr(study_hours, exam_score)
print(f"\nSpearman rho: {rho:.4f}, p-value: {p_s:.4e}")

# Correlation matrix with Pandas
df = pd.DataFrame({
    "hours": study_hours,
    "score": exam_score,
    "sleep": 8 - 0.3 * study_hours + rng.normal(0, 1, n)
})

print(f"\nCorrelation Matrix:\n{df.corr().round(3)}")

▶ Output

Pearson r: 0.8965, p-value: 2.0429e-36
Interpretation: Strong positive correlation

Spearman rho: 0.8950, p-value: 3.9817e-36

Correlation Matrix:
       hours  score  sleep
hours  1.000  0.896 -0.674
score  0.896  1.000 -0.598
sleep -0.674 -0.598  1.000

What happened here: Pearson r came out at 0.896, very close to +1, so study hours and exam score move up together in an almost straight line. The tiny p-value (2e-36, basically zero) means there is essentially no chance this strength is a fluke of random data. Spearman rho landed at 0.895, right next to Pearson, which tells you the relationship is not just strong but also nicely straight. Read the correlation matrix like a grid of every pairing: each variable scores a perfect 1.000 against itself down the diagonal, hours and score share that 0.896, and sleep shows a negative number against both because more studying eats into sleep.

One note: this data uses a fixed seed (default_rng(42)), so you will see these exact numbers. Drop the seed and your figures will shift a little, but the story stays the same.

The Math, Worked by Hand

Before you trust a library, it pays to compute correlation and a regression line once by hand, on numbers small enough to check in your head. It is like learning to add on paper before you lean on a calculator: once you have done it yourself, you know exactly what the machine is doing when you hit the button. This is the part most tutorials skip, and it is exactly the part that makes .fit() stop feeling mysterious. Every Python correlation function is doing this same arithmetic underneath.

Take three points: (1, 2), (2, 4), (3, 5). We want two things: how strongly they line up (Pearson r), and the straight line y = mx + b that fits them best.

• Mean of x: (1 + 2 + 3) / 3 = 2.0
• Mean of y: (2 + 4 + 5) / 3 = 3.6667
• For each point, take (x minus mean x) and (y minus mean y), multiply the pair, then add them up: that sum is 3.0
• Add up (x minus mean x) squared: that is 2.0. Add up (y minus mean y) squared: that is 4.6667
• Pearson r = 3.0 divided by the square root of (2.0 times 4.6667) = 0.982
• Slope m = 3.0 / 2.0 = 1.5. Intercept b = mean y minus (m times mean x) = 3.6667 minus 3.0 = 0.6667

So our best-fit line is y = 1.5x + 0.6667, and the points are 98.2% lined up. Now let us prove those exact numbers in code, using only NumPy, no statistics library doing the work for us.

📄 from_scratch.py: correlation and a regression line with plain NumPy

import numpy as np

# The same 3 points we worked by hand: (1, 2), (2, 4), (3, 5)
x = np.array([1, 2, 3], dtype=float)
y = np.array([2, 4, 5], dtype=float)

mean_x = x.mean()
mean_y = y.mean()
print(f"mean of x = {mean_x}")
print(f"mean of y = {mean_y:.4f}")

# Pearson r: covariance divided by the product of standard deviations
dx = x - mean_x
dy = y - mean_y
r = np.sum(dx * dy) / np.sqrt(np.sum(dx**2) * np.sum(dy**2))
print(f"Pearson r  = {r:.4f}")

# Least-squares line: slope first, then intercept
slope = np.sum(dx * dy) / np.sum(dx**2)
intercept = mean_y - slope * mean_x
print(f"slope (m)     = {slope}")
print(f"intercept (b) = {intercept:.4f}")

# Predictions and residuals (how far each real point sits from the line)
predicted = slope * x + intercept
residuals = y - predicted
print(f"predicted = {np.round(predicted, 4)}")
print(f"residuals = {np.round(residuals, 4)}")

▶ Output

mean of x = 2.0
mean of y = 3.6667
Pearson r  = 0.9820
slope (m)     = 1.5
intercept (b) = 0.6667
predicted = [2.1667 3.6667 5.1667]
residuals = [-0.1667  0.3333 -0.1667]

What happened here: Every number matches the hand calculation exactly: r = 0.982, slope 1.5, intercept 0.6667. The residuals line is the new piece worth staring at. A residual is just the gap between a real point and the line’s guess for it. Our three residuals are tiny (-0.17, 0.33, -0.17) and they add up to roughly zero, which is the signature of a good fit. When the library spits out a slope and intercept later, this is the arithmetic it ran for you.

Simple Linear Regression

Doing it by hand is great for understanding, painful for real data. With 50 or 50,000 points you reach for scipy.stats.linregress, which fits the same least-squares line in one call and hands back the slope, intercept, R-squared, p-value, and standard error. Think of regression as the line a tailor draws through a row of buttons: not touching every button, but sitting as close to all of them as possible.

📄 linear_regression.py: SciPy plus a from-scratch NumPy check

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
X = rng.uniform(1, 10, 50)
y = 20 + 3.5 * X + rng.normal(0, 4, 50)

# scipy.stats.linregress: the quick way
result = stats.linregress(X, y)
print(f"Slope:     {result.slope:.4f}")
print(f"Intercept: {result.intercept:.4f}")
print(f"R-squared: {result.rvalue**2:.4f}")
print(f"P-value:   {result.pvalue:.4e}")
print(f"Std error: {result.stderr:.4f}")

# Prediction: feed the line a new x, read off the y
new_x = 7.5
predicted_y = result.slope * new_x + result.intercept
print(f"\nPrediction: x={new_x} -> y={predicted_y:.1f}")

# R-squared interpretation
r2 = result.rvalue**2
print(f"\nR-squared = {r2:.4f}")
print(f"This means {r2:.1%} of the variance in y is explained by x")

# From scratch with NumPy least squares (the same line, no stats library)
A = np.column_stack([X, np.ones(len(X))])
coeffs = np.linalg.lstsq(A, y, rcond=None)[0]
print(f"\nNumPy lstsq: slope={coeffs[0]:.4f}, intercept={coeffs[1]:.4f}")

▶ Output

Slope:     3.5537
Intercept: 19.0210
R-squared: 0.8959
P-value:   3.1713e-25
Std error: 0.1748

Prediction: x=7.5 -> y=45.7

R-squared = 0.8959
This means 89.6% of the variance in y is explained by x

NumPy lstsq: slope=3.5537, intercept=19.0210

What happened here: We built the data with a true slope of 3.5 and an intercept of 20, then asked SciPy to find them back from the noisy points. It guessed a slope of 3.5537 and an intercept of 19.0210, close to the truth, with the gap coming from the random noise we sprinkled in. The R-squared of 0.8959 reads as plain English: about 90% of the up-and-down in y is explained by x, and the leftover 10% is noise and whatever we did not measure.

The standout line is the last one. The hand-rolled NumPy lstsq returns slope 3.5537 and intercept 19.0210, character for character identical to SciPy. That is the payoff of the from-scratch section: the library is not doing anything mysterious, it is running the same least-squares arithmetic you just saw.

Residual Analysis

A residual is the leftover: real value minus the line’s prediction. Residuals are how you tell a good fit from a lucky-looking one. The quick mental picture: imagine hanging laundry on a clothesline. If the line is at the right height, some shirts hang a bit above it and some a bit below, scattered evenly. If every shirt on the left droops below and every shirt on the right pokes above, the line is at the wrong angle. Residuals catch that.

📄 residuals.py: are the leftovers just random noise?

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
X = rng.uniform(1, 10, 50)
y = 20 + 3.5 * X + rng.normal(0, 4, 50)

result = stats.linregress(X, y)
predicted = result.slope * X + result.intercept
residuals = y - predicted

print(f"Mean of residuals: {residuals.mean():.6f}")   # should be about 0
print(f"Std of residuals:  {residuals.std():.4f}")
print(f"Largest residual:  {residuals.max():.4f}")
print(f"Smallest residual: {residuals.min():.4f}")

# A clean fit keeps almost every point within 2 std of the line
within_2sd = np.sum(np.abs(residuals) < 2 * residuals.std())
print(f"Points within 2 std of the line: {within_2sd} of {len(X)}")

▶ Output

Mean of residuals: 0.000000
Std of residuals:  3.0291
Largest residual:  6.6876
Smallest residual: -6.1171
Points within 2 std of the line: 48 of 50

What happened here: The mean of the residuals is 0.000000, which is no accident: the least-squares line is built so the misses above and below cancel out. The spread (standard deviation) is about 3.03, close to the noise level of 4 we baked into the data, so the line soaked up the real signal and left only the random part behind. And 48 of the 50 points sit within two standard deviations of the line, which is exactly the kind of even, boring scatter you want. If instead you saw a clear curve in the residuals, or a fan that grows wider to the right, that would be the data telling you a straight line is the wrong tool.

When Correlation Lies to You

Here is the trap that catches people running their first Python correlation analysis. Pearson r only measures straight-line strength. A relationship can be perfect and obvious to the eye, yet Pearson r reports zero, because the shape is not a straight line. The classic case is a U-shape.

📄 zero_correlation.py: a perfect relationship with r = 0

import numpy as np
from scipy import stats

# A perfect U-shape: y = x squared, x running from -5 to 5
x = np.linspace(-5, 5, 11)
y = x**2

r, p = stats.pearsonr(x, y)
print(f"x: {x}")
print(f"y: {y}")
print(f"Pearson r: {r:.4f}")
print("The relationship is perfect, yet Pearson r is 0.")

▶ Output

x: [-5. -4. -3. -2. -1.  0.  1.  2.  3.  4.  5.]
y: [25. 16.  9.  4.  1.  0.  1.  4.  9. 16. 25.]
Pearson r: 0.0000
The relationship is perfect, yet Pearson r is 0.

What happened here: Every y is exactly x squared, so the link is as strong as a link can be. But the left half slopes down and the right half slopes up, and those two halves cancel out when Pearson does its straight-line math, leaving r = 0.0000. The lesson is short and important: never judge a relationship by the correlation number alone. Plot the points first. Your eyes catch a U-shape in one second; a single correlation coefficient never will. For curved data like this, reach for polynomial regression or a non-linear model instead.

One caution: Spearman will not rescue a symmetric U-shape either, because the ranks fall and then rise, so rho also lands near zero; Spearman only helps when the curve keeps heading in one direction.

Common Mistakes

Mistake 1: Reading correlation as causation

Ice cream sales and drowning deaths rise together every summer. Ice cream does not push anyone into the water. A third thing, hot weather, drives both. That hidden third variable is called a confounder, and it is everywhere. Before you tell your boss that x drives y, ask what else moves both of them at the same time.

🚫 Wrong takeaway

# r = 0.95 between ice cream sales and drownings
# Conclusion: "ban ice cream to save lives"  (nonsense)
# The real driver is summer heat, which lifts both numbers.

✅ Right approach

# To claim x CAUSES y, you need more than correlation:
# 1. A controlled experiment or A/B test
# 2. Randomization (so other factors even out)
# 3. The cause must come before the effect in time

Mistake 2: Trusting Pearson r on curved data

Why: as the U-shape demo proved, Pearson r can be 0 on a flawless relationship just because the shape bends. Always plot the data, and switch to a polynomial or non-linear model when the cloud of points is not roughly a straight band. Spearman only helps when the curve is monotonic, rising or falling the whole way.

Mistake 3: Chasing a high R-squared and ignoring the residuals

Why: a high R-squared feels reassuring, but it can hide a curved pattern the straight line is missing. The residual plot is your honest second opinion. If the leftovers show a clear shape instead of even noise, the model is wrong no matter how high R-squared climbs.

Practice Exercises

1. Exercise 1: Generate 100 points where y rises with x plus some noise. Run both Python correlation coefficients, Pearson r and Spearman rho, and explain in one sentence why they are close.
2. Exercise 2: Fit a line with stats.linregress, then compute the residuals by hand and confirm their mean is essentially zero.
3. Exercise 3: Build a U-shaped dataset (y = x**2 with x from -5 to 5), show that Pearson r is near zero, then check Spearman rho and confirm it is near zero too, and write a note on why both fail here. Then repeat with a monotonic curve like y = np.exp(x) and watch Spearman stay at 1.0 while Pearson drops.

Conclusion

You now have the full toolkit for measuring how two variables move together and modelling that link in Python. You learned Pearson r and Spearman rho and when to reach for each, built a Python correlation matrix with pandas, worked the math by hand on three points, and then watched SciPy and a from-scratch NumPy lstsq return the exact same slope and intercept. You also saw R-squared for what it is, read residuals as an honest second opinion on the fit, and caught Pearson r reporting zero on a perfect U-shape. Above all, keep the one rule that saves careers: a strong correlation waggles its eyebrows, but it never proves causation.

Next up is turning these numbers into pictures. Regression lines, residual clouds, and correlation heatmaps all land harder when you can see them, which is exactly what the Matplotlib basics tutorial covers. For the full path from first script to machine learning, browse the Python + AI/ML tutorial series home.

Frequently Asked Questions

Interview Questions on Correlation and Regression

How interviewers actually probe this topic: real scenarios, with answers you can say out loud.

Q: What is the core difference between correlation and regression?

Correlation gives you a single number for how strongly two variables move together and in which direction, but it treats them symmetrically and predicts nothing. Regression fits an actual equation (like y = mx + b) so you can predict one variable from another, and it treats x and y differently: x is the input, y is the output. In short, correlation measures strength, regression models the relationship and lets you forecast.

Q: How is the least-squares slope actually computed?

You take the deviations of x and y from their means, multiply the paired deviations and sum them (that is the covariance term), then divide by the sum of the squared x deviations. In the post that was slope = sum(dx * dy) / sum(dx**2). The intercept then falls out of mean_y - slope * mean_x, since the least-squares line always passes through the point of means.

Q: Why does the mean of the residuals come out to essentially zero in ordinary least squares?

Least squares chooses the slope and intercept that minimise the sum of squared residuals, and the math of that minimisation forces the residuals to sum to zero whenever the model includes an intercept term. So the positive misses above the line and the negative misses below it exactly cancel. A near-zero residual mean is therefore not evidence of a good fit on its own; you still need to check the residuals for structure.

Q: Your regression reports R-squared of 0.94 on training data, but predictions on fresh data are badly off. What do you check first?

A high R-squared on the data you fit only says the line hugs those points; it says nothing about new points. First plot the residuals: a curve or a widening fan means a straight line is the wrong shape and the high R-squared is misleading. Then evaluate on a held-out test set, because you may be overfitting, and watch for a few high-leverage outliers inflating the score. R-squared measures fit, not out-of-sample accuracy.

Q: A teammate says Pearson r is 0.02 so two variables are unrelated, but the scatter plot shows an obvious arch. What is going on and what do you recommend?

Pearson r only measures straight-line strength, so a symmetric curve like a U-shape or arch can hand back r near zero even when the relationship is perfect, because the down-slope and up-slope cancel. The variables are clearly related, just not linearly. Recommend switching to Spearman rho for a monotonic trend, or fitting a polynomial or non-linear model for the curve, and always trusting the plot over a single coefficient.

Q: Your correlation matrix shows two predictors correlated at 0.95 before a regression. Why is that a problem?

That is multicollinearity: the two predictors carry almost the same information, so the model cannot cleanly separate their individual effects. The coefficients become unstable and their standard errors balloon, meaning small changes in the data can flip a coefficient’s sign or size. The usual fixes are to drop one of the pair, combine them into a single feature, or use a regularised model such as ridge regression.

Q: When would you prefer Spearman or Kendall over Pearson for a real dataset?

Reach for Spearman when the relationship is monotonic but not straight, when the data is ranked or ordinal (such as satisfaction ratings), or when outliers would distort Pearson, since Spearman works on ranks and is far more robust. Kendall tau is a good choice for small samples or when you want a coefficient that is easier to interpret as the probability of concordant versus discordant pairs. Pearson stays the right pick only when you genuinely expect a linear link between continuous, roughly normal variables.

Go deeper: when you outgrow this post, the official Python documentation is the next stop.

Related Posts

Previous: Statistics: Hypothesis Testing with t-test and chi-square

Next: Matplotlib Basics: Line, Bar, Scatter, Histogram

Series Home: Python + AI/ML Tutorial Series