Linear Regression: From the Maths to scikit-learn

Linear regression fits a straight line through your data and uses it to predict a number. It’s the oldest algorithm in the toolkit and still the right first thing to try: it trains instantly, it’s fully interpretable, and it gives you the baseline every fancier model has to beat.

This page works one example all the way through — the maths by hand first, then the same thing in scikit-learn, then the error metrics on the results.

The model

$$y = mx + c$$

• y — the value you’re predicting (the dependent variable)
• x — the input (the independent variable / feature)
• m — the slope (called the coefficient): how much y changes per unit of x
• c — the intercept: the value of y when x is 0

Training a linear regression means finding the m and c that make the line fit best. “Best” needs a definition, which is where least squares comes in.

What “best fit” means: least squares

For any candidate line, each point has a residual — the vertical gap between the actual value and the line:

$$\text{residual}_i = y_i - \hat{y}_i$$

We can’t just add residuals up (positives and negatives cancel), so we square them first and add those. The best line is the one that makes this sum as small as possible:

$$\text{minimise} \quad \sum_{i=1}^{n}(y_i - \hat{y}_i)^2$$

Hence Ordinary Least Squares (OLS). Squaring does two jobs: it kills the sign, and it punishes big misses disproportionately — being 4 off counts 16, while being 1 off counts 1. The line is pulled hard toward reducing large errors.

Solving that minimisation gives closed-form answers:

$$m = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2} \qquad c = \bar{y} - m\bar{x}$$

Worked example — experience vs salary

x = years of experience, y = salary in thousands.

x = [1, 2, 3, 4, 5]
y = [2, 4, 5, 4, 5]

By hand

Means: x̄ = 3, ȳ = 4.

x	y	x − x̄	y − ȳ	(x−x̄)(y−ȳ)	(x−x̄)²
1	2	−2	−2	4	4
2	4	−1	0	0	1
3	5	0	1	0	0
4	4	1	0	0	1
5	5	2	1	2	4
			Σ	6	10

$$m = \frac{6}{10} = 0.6 \qquad c = 4 - 0.6 \times 3 = 2.2$$

The fitted line is y = 0.6x + 2.2.

Read it in plain words: each extra year of experience is worth about 0.6k in salary, and someone with zero experience starts around 2.2k. That interpretability is linear regression’s biggest selling point — the coefficient means something.

In scikit-learn

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

x = [1, 2, 3, 4, 5]
y = [2, 4, 5, 4, 5]

# Reshape x into a 2-D column — scikit-learn requires 2-D features
x = np.array(x).reshape(-1, 1)
y = np.array(y)

model = LinearRegression()
model.fit(x, y)

print("Slope (Coefficient):", model.coef_)      # => [0.6]
print("Intercept:", model.intercept_)           # => 2.2

Identical to the hand calculation.

Why .reshape(-1, 1)? Scikit-learn always expects X as a 2-D array of shape (n_samples, n_features). A plain list of 5 numbers has shape (5,); reshaping makes it (5, 1) — five samples, one feature each. The -1 means “work out this dimension yourself”.

print(x)
# [[1]
#  [2]
#  [3]
#  [4]
#  [5]]

Forgetting this is the single most common scikit-learn error for beginners: “Expected 2D array, got 1D array instead.” Note y stays 1-D — only X gets reshaped.

Predicting

y_pred = model.predict(x)
print(y_pred)
# [2.8 3.4 4.  4.6 5.2]

Check one by hand: at x = 1, 0.6(1) + 2.2 = 2.8 ✓

Plotting it

plt.scatter(x, y, color='red', label='Actual data')
plt.plot(x, y_pred, color='blue', label='Fitted line')
plt.scatter(x, y_pred, color='green', marker='o', label='Predicted data')

# Annotate each predicted point with its coordinates
for i in range(len(x)):
    plt.annotate(f'({x[i][0]},{y_pred[i]})', (x[i][0], y_pred[i]),
                 textcoords="offset points", xytext=(0, 10), ha='center')

plt.title('x vs y (Linear Regression)')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

Measuring the error

The residuals from our fit:

x	Actual y	Predicted ŷ	Residual (y − ŷ)	|Residual|	Residual²
1	2	2.8	−0.8	0.8	0.64
2	4	3.4	+0.6	0.6	0.36
3	5	4.0	+1.0	1.0	1.00
4	4	4.6	−0.6	0.6	0.36
5	5	5.2	−0.2	0.2	0.04
				Σ 3.2	Σ 2.40

Mean Absolute Error (MAE)

$$\text{MAE} = \frac{1}{n}\sum |y_i - \hat{y}_i| = \frac{3.2}{5} = 0.64$$

The average miss, in the original units. Our predictions are off by about 0.64k on average. MAE is the most directly interpretable metric — and because it doesn’t square, it’s robust to outliers.

Mean Squared Error (MSE)

Three steps: subtract, square, average.

$$\text{MSE} = \frac{1}{n}\sum (y_i - \hat{y}_i)^2 = \frac{2.40}{5} = 0.48$$

MSE punishes large errors harder, which is what you want when one huge miss is worse than several small ones. The catch is the units: this is 0.48 squared thousands, which means nothing physically.

Root Mean Squared Error (RMSE)

$$\text{RMSE} = \sqrt{\text{MSE}} = \sqrt{0.48} = 0.693$$

The square root puts it back into the original units — so RMSE is the metric to report. It keeps MSE’s sensitivity to big errors while staying readable.

import numpy as np
from sklearn.metrics import mean_absolute_error, mean_squared_error

mae = mean_absolute_error(y, y_pred)
mse = mean_squared_error(y, y_pred)
rmse = np.sqrt(mse)

print(f"Mean Absolute Error (MAE): {mae}")   # 0.6399999999999999
print(f"Mean Squared Error (MSE): {mse}")    # 0.47999999999999987
print(f"RMSE: {rmse}")                       # 0.6928203230275508

# MSE by hand gives exactly the same number:
mse_manual = np.mean((y - y_pred) ** 2)      # 0.47999999999999987

Is 0.693 good? Only relative to the scale of the target. Our y values run up to 5, so being off by ~0.69 is roughly a 14% error — reasonable for five data points. The same RMSE on a target scaled in millions would be superb; on a target running 0–1 it would be catastrophic. RMSE is never good or bad in isolation.

R² — the score

$$R^2 = 1 - \frac{\sum(y_i - \hat{y}_i)^2}{\sum(y_i - \bar{y})^2} = 1 - \frac{2.4}{6.0} = 0.60$$

R² is the fraction of the variance in y that the model explains, on a 0–1 scale:

• 1.0 — perfect fit
• 0.60 — our model explains 60% of the variation; the rest is unexplained
• 0.0 — no better than always predicting the mean
• negative — worse than predicting the mean

print(model.score(x, y))    # => 0.6

R² is the one metric that’s comparable across problems, because it’s unitless.

Which metric to report

Metric	Units	Outlier-sensitive	Use when
MAE	original	Low	You want the plain average error
MSE	squared	High	Optimising (it’s the training loss)
RMSE	original	High	Reporting — the usual default
R²	none	Medium	Comparing across different datasets

A complete workflow on a real dataset

The example above fit the model on all five points and scored it on those same points. That’s fine for learning the maths, but it isn’t how you build a model — you’d have no idea whether it generalises.

Here’s the same algorithm done properly on a real 170-row dataset, with a train/test split, scaling, and evaluation on data the model never saw.

Step 1 — Load and look before modelling

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

dataset = pd.read_csv('weight-Height.csv')
dataset.head()

   Weight  Height
0      74     242
1      69     162
2      74     213
3      72     220
4      70     206

Always plot the relationship first. Linear regression assumes a straight line — a scatter plot tells you in two seconds whether that’s plausible:

plt.scatter(dataset['Weight'], dataset['Height'])
plt.xlabel('Weight')
plt.ylabel('Height')

Sanity-check your columns while you’re here. In this file the values look swapped relative to their names — 74 with 242 reads like 74 inches tall, 242 lbs, not the reverse. The maths works either way, but the interpretation of your coefficient depends on knowing which is which. Verifying that column names match their contents is part of exploratory analysis, and it costs nothing to check.

Step 2 — Build X and y with the right shapes

Scikit-learn wants X as 2-D and y as 1-D. Pandas gives you both, depending on how you index:

# Double brackets -> DataFrame -> 2-D. This is what X needs.
x = dataset[['Weight']]

# Single brackets -> Series -> 1-D. This is what y needs.
y = dataset['Height']
y = np.array(y)          # convert the Series to a plain 1-D array

x.shape    # => (170, 1)   rows AND columns
y.shape    # => (170,)     rows only

df[['col']] vs df['col'] is the whole trick — the extra pair of brackets is the difference between a 2-D DataFrame and a 1-D Series. If you already have a Series, .reshape(-1, 1) gets you to the same place:

x_series = dataset['Weight']                    # 1-D Series
x_series = np.array(x_series).reshape(-1, 1)    # -> 2-D array
x_series.shape    # => (170, 1)

Both routes produce identical input. Checking .shape at this point is the cheapest way to avoid the “Expected 2D array, got 1D array instead” error.

Step 3 — Split before you touch anything else

from sklearn.model_selection import train_test_split

x_train, x_test, y_train, y_test = train_test_split(
    x, y, test_size=0.2, random_state=1
)

x_train.shape    # => (136, 1)    80% of 170

170 rows → 136 train, 34 test.

Step 4 — Standardise (z-score normalisation)

from sklearn.preprocessing import StandardScaler

sc = StandardScaler()
x_train = sc.fit_transform(x_train)   # LEARN mean/SD from train, then apply
x_test = sc.transform(x_test)         # apply the SAME mean/SD to test

Each value becomes a z-score — how many standard deviations it sits from the mean:

$$z = \frac{x - \mu}{\sigma}$$

x_train
# array([[ 0.52709199],
#        [-0.2437092 ],
#        [-0.62910979],
#        [ 1.29789317],
#        ...

Values now centre on 0 and mostly fall in −3…+3.

fit_transform on train, transform on test — never fit_transform on both. The scaler must be calibrated on training data only; scaling the test set with its own statistics leaks information about it into your evaluation and quietly inflates your score. This is the data leakage trap, and it’s the single most common mistake in a beginner pipeline.

Step 5 — Fit

from sklearn.linear_model import LinearRegression

regressor = LinearRegression()
regressor.fit(x_train, y_train)

print('coefficient or slope: ', regressor.coef_)      # => [15.05079326]
print('intercept: ', regressor.intercept_)            # => 184.88235294117644

Reading coefficients after scaling

This is where scaling changes the meaning of your numbers, and it trips people up constantly.

Unscaled, a coefficient reads as “per one unit of x”. Scaled, it reads as “per one standard deviation of x”:

15.05 — a one-standard-deviation increase in Weight is associated with a 15.05 unit increase in Height.

That’s the trade-off: you lose the plain-units interpretation, but you gain comparability — with several features, the largest scaled coefficient is the most influential one, which is not true of unscaled coefficients.

There’s also a neat consequence hiding in the intercept. Standardised data has mean 0, and an OLS line always passes through the point of means — so at x = 0 the prediction is the training mean of y:

intercept 184.88 = the mean Height of the training set.

Whenever you standardise your features, the intercept stops being “the value at x = 0” in any physical sense and becomes “the average outcome”. A useful check that your scaling did what you thought.

Step 6 — Predict and plot the fit

y_predict_train = regressor.predict(x_train)

plt.scatter(x_train, y_train, color='red')
plt.plot(x_train, y_predict_train, color='blue')
plt.title('Training data')
plt.xlabel('Weight')
plt.ylabel('Height')
plt.show()

Red points are the actual training data; the blue line is the model. Note the x-axis is now in standard deviations, not the original units — a side effect of plotting scaled data.

Step 7 — Evaluate on the test set

This is the number that actually matters:

y_predict_test = regressor.predict(x_test)

from sklearn.metrics import mean_squared_error, mean_absolute_error

mse = mean_squared_error(y_test, y_predict_test)
mae = mean_absolute_error(y_test, y_predict_test)
rmse = np.sqrt(mse)

print('Mean Squared Error: ', mse)        # 169.1522176510848
print('Mean Absolute Error: ', mae)       # 10.666581636991307
print('Root Mean square Error:', rmse)    # 13.005853207348022

Metric	Value	Reading
MSE	169.15	Squared units — for optimising, not reporting
MAE	10.67	Typical miss of about 10.7 units
RMSE	13.01	About 13 units, penalising big misses more

Is that good? Judge it against the scale of the target. The Height values run roughly 160–240, averaging around 200 — so an RMSE of 13 is about a 6–7% error. Reasonable for a single predictor.

Note also that RMSE (13.0) exceeds MAE (10.7). That gap is always present when errors vary in size, and the wider it is, the more a few large misses are driving your error — a quick diagnostic for outliers, as covered in skewness and kurtosis.

plt.scatter(x_test, y_test, color='red')
plt.plot(x_test, y_predict_test, color='blue')
plt.title('Test data')
plt.xlabel('Weight')
plt.ylabel('Height')

The test plot is the honest picture: those red points were never seen during training, so how close the blue line runs to them is real predictive performance.

The whole thing, start to finish

import pandas as pd, numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

# 1. Load
dataset = pd.read_csv('weight-Height.csv')
x = dataset[['Weight']]          # 2-D
y = np.array(dataset['Height'])  # 1-D

# 2. Split FIRST
x_train, x_test, y_train, y_test = train_test_split(
    x, y, test_size=0.2, random_state=1
)

# 3. Scale using train statistics only
sc = StandardScaler()
x_train = sc.fit_transform(x_train)
x_test = sc.transform(x_test)

# 4. Fit
regressor = LinearRegression()
regressor.fit(x_train, y_train)

# 5. Evaluate on unseen data
y_pred = regressor.predict(x_test)
print('RMSE:', np.sqrt(mean_squared_error(y_test, y_pred)))
print('MAE :', mean_absolute_error(y_test, y_pred))
print('R2  :', r2_score(y_test, y_pred))

Multiple linear regression

With more than one feature the line becomes a plane (or hyperplane) — but nothing about the code changes:

$$y = b_0 + b_1x_1 + b_2x_2 + \dots + b_nx_n$$

# X now has several columns; y is unchanged
model = LinearRegression()
model.fit(X_train, y_train)

print(model.coef_)       # one coefficient per feature
print(model.intercept_)

Each coefficient reads as: “holding everything else constant, a one-unit increase in this feature changes y by this much.” That’s why scaling matters for interpretation — unscaled coefficients aren’t comparable to each other.

The assumptions

Linear regression is only valid when these roughly hold:

Assumption	Meaning	If violated
Linearity	The relationship really is a straight line	Use polynomial or a non-linear model
Independence	Observations don’t influence each other	Time series needs different methods
Homoscedasticity	Residual spread is constant across x	Transform y (often a log)
Normality of residuals	Errors are roughly normally distributed	Check with a residual histogram
No multicollinearity	Features aren’t strongly correlated	Drop one, or use Ridge

The quickest diagnostic is a residual plot — predictions on x, residuals on y. You want a formless cloud around zero. Any curve, funnel or pattern means an assumption is broken.

residuals = y - y_pred
plt.scatter(y_pred, residuals)
plt.axhline(0, color='red', linestyle='--')
plt.xlabel('Predicted'); plt.ylabel('Residual')
plt.show()

Takeaways

• Linear regression fits y = mx + c by minimising the sum of squared residuals (least squares).
• Our worked fit: y = 0.6x + 2.2 — every extra year of experience is worth ~0.6k, and the coefficient is directly interpretable.
• X must be 2-D, y 1-D. Use df[['col']] (double brackets) or .reshape(-1, 1), and check .shape before you fit.
• Split before you scale, and fit_transform the scaler on training data only — otherwise you leak the test set into training.
• After standardising, a coefficient means “per one standard deviation”, and the intercept becomes the mean of y.
• Evaluate on the test set. On the real dataset: RMSE 13.0 against a target averaging ~200, so roughly a 6–7% error.
• MAE = average miss (robust); MSE = punishes big errors (training loss); RMSE = MSE back in real units (report this); R² = fraction of variance explained (compare across problems). A gap between RMSE and MAE points to a few large errors.
• A metric is only meaningful relative to the scale of your target.
• Plot the scatter before fitting and the residuals after.

Next: regression algorithms beyond the straight line.