Linear Regression in Python: Study Hours vs Score With SciPy
Linear regression in Python takes one function call: scipy.stats.linregress(x, y). It fits a straight line through your data and hands back the slope, the intercept, and a few numbers that tell you how good the line is.
In this post I use it on made-up data: hours studied against exam score. Then I open the box. I’ll show what R-squared measures, how least squares finds the line, why a good line can still give a silly prediction, what changes when you add a second input, and how the same fit looks in scikit-learn.
Linear regression shows up in my Big Data Analytics master’s at SDSU: I used it in one of my final projects, before the AI boom, when you built and tested your own models. This post pulls that together with notes from a machine learning course on Udemy.
The example is mine: study hours against score. Every code block runs top to bottom on NumPy 2.2.6 and SciPy 1.15.3, and the output shown is what it printed. If arrays are new to you, start with my NumPy tutorial first.
What is linear regression?
Linear regression predicts a number, called the dependent variable y, from one or more inputs, called independent variables x. It assumes the relationship is a straight line:
y = b0 + b1 * x + e
b0is the intercept: the predictedywhenxis zero.b1is the slope: how muchychanges whenxgoes up by one.eis the error: everything aboutythatxdoesn’t explain.
You already do this in your head. Say you earned $10,000 two years ago and $20,000 last year. If you guess $30,000 for this year, you just drew a line through two points and extended it.
The model answers two questions. Is there a linear relationship between x and y? And, with several inputs, which one matters most?
It has a big catch. It only works when the relationship really is roughly linear.
And it will never predict every point exactly. The error term is there for a reason.
How do you run linear regression in Python?
Install what you need:
pip install numpy scipy matplotlib
I generated 100 fake students so I’d know the true answer. Hours studied is normally distributed around 20. Each hour adds 1.1 points to the score, which is out of 100.
There’s also a second input, exercises completed, that adds 0.7 points per exercise, plus random noise.
For now I’ll ignore exercises and regress score on hours alone.
import numpy as np
from scipy import stats
rng = np.random.default_rng(1)
hours = rng.normal(20, 6, 100)
exercises = rng.normal(15, 4, 100)
score = 25 + 1.1 * hours + 0.7 * exercises + rng.normal(0, 4, 100)
print(f"hours: min {hours.min():.1f}, max {hours.max():.1f}")
print(f"score: min {score.min():.1f}, max {score.max():.1f}")
slope, intercept, r_value, p_value, std_err = stats.linregress(hours, score)
print(f"slope {slope:.3f} intercept {intercept:.3f} r {r_value:.3f} r2 {r_value**2:.3f} p {p_value:.2e} stderr {std_err:.3f}")
hours: min 3.7, max 32.7
score: min 39.2, max 75.4
slope 1.077 intercept 35.103 r 0.770 r2 0.593 p 7.28e-21 stderr 0.090
Here’s what each number means. The SciPy docs list them all.
- slope, 1.077: each extra hour of study goes with about one more point. That’s close to the 1.1 I built in. It says “goes with”, not “causes”.
- intercept, 35.103: the predicted score at zero hours. Nobody in this data studied less than 3.7 hours, so don’t read too much into it.
- r, 0.770: the correlation between hours and score.
- p, 7.28e-21: the p-value for a test where the starting assumption is that the slope is zero. A tiny value means that if the true slope were zero, you would almost never see a slope this large by chance.
- stderr, 0.090: the uncertainty on the slope.
What does R-squared tell you?
R-squared, also called the coefficient of determination, says how much of the variation in y your line captures. Here it’s 0.593, so the line explains about 59% of the spread in scores. The rest is everything else: exercises, sleep, a lucky guess.
SciPy’s docs confirm the shortcut: the square of rvalue is the coefficient of determination. But the formula is worth seeing once:
R-squared = 1 - (sum of squared errors) / (sum of squared distances from the mean)
The bottom part is the total variation in scores around their own average. The top part is what the line failed to explain. In code:
def predict(x):
return slope * x + intercept
sse = np.sum((score - predict(hours)) ** 2)
sst = np.sum((score - score.mean()) ** 2)
print("1 - sse/sst =", round(1 - sse / sst, 3))
1 - sse/sst = 0.593
Same number. Two ways to read it:
- 0 means the line explains nothing. 1 means it explains everything. A low value is a poor fit. A high value is a good one.
- It isn’t strictly stuck between 0 and 1. On new data, or for a model fitted without an intercept, R-squared can go negative. That means the model does worse than just guessing the average.
A high R-squared doesn’t prove your model is right. It just says the line follows the data you fitted it to.
How does least squares find the line?
The method is called ordinary least squares, or OLS. It picks the slope and intercept that make the total squared error as small as possible.
The error for one student is the actual score minus the predicted one. Then we square it, for two reasons:
- A prediction can miss high or low. Squaring makes every miss positive, so misses can’t cancel each other out.
- Squaring punishes big misses far more than small ones.
Squaring doesn’t guarantee a better model. It just defines what “best line” means.
You never have to do the math, because Python does it. But knowing what’s under the hood helps.
The answer has a closed form. For one input it looks like this:
slope = r * (standard deviation of y) / (standard deviation of x)
intercept = mean of y - slope * mean of x
I checked it against SciPy:
slope_hand = r_value * score.std() / hours.std()
intercept_hand = score.mean() - slope_hand * hours.mean()
print("by hand:", round(slope_hand, 3), round(intercept_hand, 3))
by hand: 1.077 35.103
Identical. So OLS gets its answer in one step.
Some tutorials show the line inching toward the data over many iterations. That’s gradient descent, which is how neural networks are trained, and it isn’t needed for a straight line.
I use it in my overfitting vs underfitting post if you want to see it running.
How do you predict with the fitted line?
Multiply the input by the slope and add the intercept. I already wrote predict() above. For a student who studies 30 hours:
print("30 hours ->", round(predict(30), 1))
print("100 hours ->", round(predict(100), 1))
30 hours -> 67.4
100 hours -> 142.8
The first answer is reasonable. 30 hours sits inside the range of the data, which ran up to 32.7.
The second answer is nonsense. A score of 142.8 isn’t possible on a test out of 100.
The math worked perfectly. The line just has no idea that scores stop at 100, and it has never seen anyone study for 100 hours.
Only trust predictions inside the range you fitted.
To draw the data and the fitted line with matplotlib:
import matplotlib.pyplot as plt
plt.scatter(hours, score)
plt.plot(hours, predict(hours), c="r")
plt.xlabel("Hours studied")
plt.ylabel("Score")
plt.show()
What is multiple linear regression?
Multiple linear regression uses two or more inputs. Add exercises completed to the model and it becomes:
score = b0 + b1 * hours + b2 * exercises
Now there are three numbers to estimate, and linregress can’t do it, since it only takes one x. NumPy’s least squares solver does the same job for any number of inputs. Add a column of ones so the model gets an intercept:
X = np.column_stack([np.ones_like(hours), hours, exercises])
coef, *_ = np.linalg.lstsq(X, score, rcond=None)
b0, b1, b2 = coef
print(f"b0 {b0:.3f} b1 {b1:.3f} b2 {b2:.3f}")
pred = X @ coef
r2 = 1 - np.sum((score - pred) ** 2) / sst
print("multiple r2", round(r2, 3))
print("30h, 15 ex ->", round(b0 + b1 * 30 + b2 * 15, 1))
b0 24.737 b1 1.091 b2 0.686
multiple r2 0.739
30h, 15 ex -> 67.8
The estimates land close to the truth I built in: 25, 1.1, and 0.7. R-squared jumps from 0.593 to 0.739, because exercises really do matter in this data.
Each hour now adds about 1.09 points with exercises held fixed. In this data, hours and exercises are almost uncorrelated (-0.03), which is why the one-input slope of 1.077 was already close to the truth.
In real data, inputs usually overlap, and leaving one out shifts the slope of the other.
One warning. Adding an input can never lower R-squared on the data you fitted, even if that input is pure noise.
So a rising R-squared isn’t proof that a new variable belongs. Check on data the model hasn’t seen, which is the whole point of the overfitting vs underfitting comparison.
Where have I used linear regression for real?
In the SDSU program, my part of the BDA 594 final group project was about NFL special teams. We were given league data and asked to build new metrics for the specialists: punters, kickers, and long snappers. Then we ranked them, hoping to find something that could help a team perform better.
I preprocessed the exploratory data and built predictions for two things: the average distance a special teams player travels during a game, and their kick success rate in a game. The model was a simple linear regression, built as an artificial neural network.
That’s less odd than it sounds. A neural network with no hidden layers and no activation function is linear regression: one weight per input, one bias, and a weighted sum.
What changes is how it’s fitted. The network nudges its weights a little on every pass, which is gradient descent, instead of using the closed-form formula from earlier.
Both are chasing the same line.
Ranking real people from predicted values is where the warnings above stop being academic. R-squared tells you how much of the spread your line explains.
It doesn’t tell you whether a ranking is fair. And a prediction is only as good as the range the model was fitted on, so a player far outside that range gets a guess, not a prediction.
My grad projects are in this GitHub repo, with a short write-up of each one.
How do you do linear regression with scikit-learn?
Most tutorials on this topic use scikit-learn’s LinearRegression instead of SciPy. It gives the same answer. Install it with pip install scikit-learn, then reuse the hours, exercises, and score arrays from above.
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(hours.reshape(-1, 1), score)
print(f"slope {model.coef_[0]:.3f} intercept {model.intercept_:.3f} r2 {model.score(hours.reshape(-1, 1), score):.3f}")
print("30 hours ->", round(model.predict([[30]])[0], 1))
slope 1.077 intercept 35.103 r2 0.593
30 hours -> 67.4
scikit-learn wants the inputs as a 2D array, with one row per student and one column per input. That’s why hours gets reshaped. Multiple regression is the same call with two columns:
X = np.column_stack([hours, exercises])
multi = LinearRegression().fit(X, score)
print(f"b0 {multi.intercept_:.3f} b1 {multi.coef_[0]:.3f} b2 {multi.coef_[1]:.3f} r2 {multi.score(X, score):.3f}")
b0 24.737 b1 1.091 b2 0.686 r2 0.739
The slope, intercept, and R-squared match SciPy and NumPy to three decimals. So which one should you use?
- SciPy’s
linregressis the quickest for one input, and it also gives you the p-value and standard error. - scikit-learn is the better fit once you have several inputs, or when the model is one step in a bigger workflow. It has the same
fitandpredictpattern as the classifier in my KNN tutorial, so train and test splits and cross-validation work the same way.
The LinearRegression docs list the options.
What should you remember about linear regression in Python?
scipy.stats.linregressgives you the slope, intercept, correlation, and p-value in one call.- R-squared is the share of variation your line explains, and
1 - SSE/SSTcomputes it. - OLS has a closed-form solution, so no iterations are needed.
- Predict inside your data’s range. Outside it, the line is guessing.
- For more than one input, use
np.linalg.lstsqor scikit-learn’sLinearRegression.
What’s the first thing you’d try to predict with a straight line?