Linear Regression with Python

Linear regression is one of the oldest and most widely used predictive techniques in statistics and machine learning. At its core, it answers a simple question: given a set of measurable inputs, what continuous output value does the data suggest? Unlike classification, which predicts categories, regression predicts a number — for example, the selling price of a house.

In this tutorial you will build a regression pipeline on the Boston housing dataset. You will explore the data, select features by correlation, train a LinearRegression model with scikit-learn, and evaluate it using three standard metrics. You will also plot learning curves to understand how the model generalises as more training data is added.

Prerequisites: Python 3.x, scikit-learn, Pandas, NumPy, Matplotlib, Seaborn.

What is Regression?

Regression is a statistical method that attempts to determine the strength and character of the relationship between one dependent variable (the value you want to predict, usually written as $y$) and one or more independent variables (the inputs, written as $x$).

The diagram below shows the core idea: raw data goes into a regression model, which learns the $x\to y$ relationship and produces a continuous prediction.

Where Regression Applies

Linear regression is used across many domains because it is fast to train, easy to interpret, and often gives a strong baseline. A few practical examples follow.

Stock price prediction uses historical prices and news events as input features to forecast future stock values.

Tweet popularity prediction estimates how many people will retweet a post based on the number of followers, hashtag popularity, and text features.

House price prediction uses property attributes such as number of rooms, distance to employment centres, and neighbourhood crime rate to estimate sale price — the exact task you will work through in this tutorial.

Simple and Multiple Linear Regression

Linear regression comes in two forms depending on how many input features you use.

Simple Linear Regression

Simple linear regression uses a single independent variable to predict the dependent variable. The relationship is described by a straight line:

Multiple Linear Regression

Multiple linear regression extends the formula to more than one independent variable, allowing the model to capture more complex patterns in the data. Instead of fitting a line in 2D space, you fit a hyperplane in higher-dimensional space. When you add too many features, you can also fit more complex polynomial relationships.

The diagram below illustrates how adding a second input dimension (e.g. number of bathrooms alongside house size) lifts the problem from a 2D line to a 3D plane fit.

Theory: Best-Fit Line, Cost Function, and Gradient Descent

The Best-Fit Line and Residuals

The best-fit line (also called the regression line) is the line for which the error between predicted values and observed values is minimised. These errors are called residuals — the vertical distances between each observed data point and the line.

The Cost Function

To find the best-fit line algorithmically, you minimise a cost function — a measure of how wrong your current parameter estimates are. The standard choice is Mean Squared Error (MSE), which penalises large errors more heavily than small ones.

The diagram below shows the hypothesis function, the MSE cost function $J({\theta }_0,{\theta }_1)$, and the optimisation goal:

Gradient Descent

Gradient descent is the optimisation algorithm used to minimise the cost function. Starting from random parameter values for ${\theta }_0$ and ${\theta }_1$, gradient descent iteratively takes steps in the direction of steepest descent on the cost surface until it reaches the minimum. The step size is controlled by the learning rate — too large and the algorithm overshoots; too small and training is very slow.

The diagram below shows both the goodness-of-fit metric being defined for each candidate line and the gradient descent surface being navigated to reach optimal parameters:

Bias-Variance Tradeoff

Every supervised model must balance two sources of error. Bias refers to simplifying assumptions a model makes to learn the target function — high bias produces underfitting. Variance refers to how sensitive the model is to fluctuations in the training data — high variance produces overfitting. The goal is to reach the sweet spot where both are low.

The diagram below shows the classic bias-variance tradeoff curve where MSE = bias² + variance, and the optimal model complexity sits at the "sweet spot":

Loading the Boston Dataset

Start by importing all libraries you will need for this tutorial. Keeping all imports in one block at the top makes dependencies immediately visible.

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error

Load the dataset and inspect its structure:

boston = load_boston()
type(boston)
sklearn.utils.Bunch
boston.keys()

dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])

Print the full dataset description to understand what each feature represents:

print(boston.DESCR)

.. _boston_dataset:

Boston house prices dataset
---------------------------

**Data Set Characteristics:**

    :Number of Instances: 506

    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

    :Attribute Information (in order):
        - CRIM     per capita crime rate by town
        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
        - INDUS    proportion of non-retail business acres per town
        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
        - NOX      nitric oxides concentration (parts per 10 million)
        - RM       average number of rooms per dwelling
        - AGE      proportion of owner-occupied units built prior to 1940
        - DIS      weighted distances to five Boston employment centres
        - RAD      index of accessibility to radial highways
        - TAX      full-value property-tax rate per $10,000
        - PTRATIO  pupil-teacher ratio by town
        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
        - LSTAT    % lower status of the population
        - MEDV     Median value of owner-occupied homes in $1000's

    :Missing Attribute Values: None

    :Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/

This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic
prices and the demand for clean air', J. Environ. Economics & Management,
vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics
...', Wiley, 1980.   N.B. Various transformations are used in the table on
pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression
problems.

Inspect the feature names and the first five target values (median house prices in thousands of dollars):

boston.feature_names

array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
       'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')

boston.target[: 5]

array([24. , 21.6, 34.7, 33.4, 36.2])

Inspect the raw array shape to confirm 506 samples and 13 features:

data = boston.data
type(data)
numpy.ndarray
data.shape

(506, 13)

Building the DataFrame

A DataFrame is a two-dimensional tabular data structure provided by Pandas. Converting the NumPy array into a DataFrame gives you labelled columns and makes exploration much easier.

data = pd.DataFrame(data = data, columns= boston.feature_names)
data.head()

Add the target (house price) as a new column so all data lives in one place:

data['Price'] = boston.target
data.head()

Use describe() to get a statistical summary of every feature including count, mean, standard deviation, and percentiles:

data.describe()

data.info() gives a concise structural summary — column types, non-null counts, and memory usage — useful for spotting missing values quickly:

data.info()

RangeIndex: 506 entries, 0 to 505
Data columns (total 14 columns):
 #   Column   Non-Null Count  Dtype
---  ------   --------------  -----
 0   CRIM     506 non-null    float64
 1   ZN       506 non-null    float64
 2   INDUS    506 non-null    float64
 3   CHAS     506 non-null    float64
 4   NOX      506 non-null    float64
 5   RM       506 non-null    float64
 6   AGE      506 non-null    float64
 7   DIS      506 non-null    float64
 8   RAD      506 non-null    float64
 9   TAX      506 non-null    float64
 10  PTRATIO  506 non-null    float64
 11  B        506 non-null    float64
 12  LSTAT    506 non-null    float64
 13  Price    506 non-null    float64
dtypes: float64(14)
memory usage: 55.5 KB

The dataset has no missing values, which you can confirm explicitly:

data.isnull().sum()

CRIM       0
ZN         0
INDUS      0
CHAS       0
NOX        0
RM         0
AGE        0
DIS        0
RAD        0
TAX        0
PTRATIO    0
B          0
LSTAT      0
Price      0
dtype: int64

Exploratory Data Analysis

Before selecting features, you should understand how they relate to each other and to the target. A scatterplot matrix (pairplot) lets you see all pairwise relationships and identify obvious linear correlations at a glance.

sns.pairplot(data)
plt.show()

The pairplot below shows every feature plotted against every other feature. Diagonal cells show the marginal distribution of each variable:

You can also inspect the marginal distribution of each individual feature with distribution plots to spot skew and outliers:

rows = 2
cols = 7
fig, ax = plt.subplots(nrows= rows, ncols= cols, figsize = (16,4))
col = data.columns
index = 0
for i in range(rows):
    for j in range(cols):
        sns.distplot(data[col[index]], ax = ax[i][j])
        index = index + 1
plt.tight_layout()
plt.show()

The grid of distribution plots below shows the density of each feature — several are heavily right-skewed (e.g. CRIM, ZN), while RM and Price are roughly bell-shaped:

Correlation Matrix

A correlation matrix quantifies the pairwise linear relationships between all variables using Pearson's r — a coefficient that ranges from −1 (perfect negative correlation) to +1 (perfect positive correlation). You will use this matrix to select which features to feed into the model.

corrmat = data.corr()
corrmat

A heatmap makes it much easier to spot high-correlation pairs visually. Warm colours indicate positive correlation and cool colours indicate negative correlation. Run the block below to render the annotated heatmap:

fig, ax = plt.subplots(figsize = (18, 10))
sns.heatmap(corrmat, annot = True, annot_kws={'size': 12})
plt.show()

The output heatmap shows the full Pearson correlation matrix for all 14 features. Notice that RM (average number of rooms) has the strongest positive correlation with Price (0.70) and LSTAT (lower-status population percentage) has the strongest negative correlation (−0.74).

Feature Selection by Correlation Threshold

Rather than using all 13 features, you will select only those with an absolute Pearson correlation above a chosen threshold against the target price. The helper function below returns a DataFrame of features that pass the threshold:

corrmat.index.values

array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
       'TAX', 'PTRATIO', 'B', 'LSTAT', 'Price'], dtype=object)

def getCorrelatedFeature(corrdata, threshold):
    feature = []
    value = []

    for i, index in enumerate(corrdata.index):
        if abs(corrdata[index])> threshold:
            feature.append(index)
            value.append(corrdata[index])

    df = pd.DataFrame(data = value, index = feature, columns=['Corr Value'])
    return df

Apply the function at a threshold of 0.50 to find features with a moderate-to-strong relationship with Price:

threshold = 0.50
corr_value = getCorrelatedFeature(corrmat['Price'], threshold)
corr_value

Three features pass the 0.50 threshold: RM (positive), PTRATIO (negative), and LSTAT (negative). Retrieve their names as an array:

corr_value.index.values

array(['RM', 'PTRATIO', 'LSTAT', 'Price'], dtype=object)

Subset the DataFrame to these correlated features and inspect the first few rows:

correlated_data = data[corr_value.index]
correlated_data.head()

Plot the pairwise relationships within this reduced feature set to confirm the linear trends are visible:

sns.pairplot(correlated_data)
plt.tight_layout()

Training and Evaluating the Model

With the features selected, split the data into training (80 %) and test (20 %) sets. Shuffling the data when splitting removes any ordering bias in the dataset:

X = correlated_data.drop(labels=['Price'], axis = 1)
y = correlated_data['Price']
X.head()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state = 0)
X_train.shape, X_test.shape

((404, 3), (102, 3))

You have 404 samples for training and 102 for testing. Now fit the model and generate predictions on the test set:

model = LinearRegression()
model.fit(X_train, y_train)
LinearRegression()
y_predict = model.predict(X_test)
df = pd.DataFrame(data = [y_predict, y_test])
df.T

102 rows × 2 columns

Regression Evaluation Metrics

Good regression models require a quantitative measure of prediction quality. Three metrics are widely used, each with a different sensitivity to large errors.

Mean Absolute Error (MAE) is the average of the absolute differences between predictions and actuals — easy to interpret because it is in the same units as the target:

Mean Squared Error (MSE) squares each error before averaging, so large errors are penalised much more heavily than small ones:

Root Mean Squared Error (RMSE) takes the square root of MSE so the metric is back in the same units as the target, making it more interpretable:

In practice, MAE is easiest to interpret as the average error, MSE punishes large mistakes more, and RMSE gives you a single number in the units of your target variable. All three are loss functions — you want to minimise them.

The R² Score

The coefficient of determination, $R^2$, measures how much of the variance in the target the model explains. It is defined as:

The diagram below illustrates the difference between a bad model (R² near 0) and a good model (R² near 1):

Now compute all three metrics for the first model:

from sklearn.metrics import r2_score

correlated_data.columns

Index(['RM', 'PTRATIO', 'LSTAT', 'Price'], dtype='object')

score = r2_score(y_test, y_predict)
mae = mean_absolute_error(y_test, y_predict)
mse = mean_squared_error(y_test, y_predict)
print('r2_score: ', score)
print('mae: ', mae)
print('mse: ', mse)

r2_score:  0.48816420156925056
mae:  4.404434993909258
mse:  41.67799012221684

An R² of 0.49 means the three-feature model explains about half the variance in house prices. You will now experiment with different correlation thresholds to see whether adding or removing features improves this.

Feature Selection Experiments

The helper function below records metrics for each experiment so you can compare all threshold configurations in one table:

total_features = []
total_features_name = []
selected_correlation_value = []
r2_scores = []
mae_value = []
mse_value = []

def performance_metrics(features, th, y_true, y_pred):
    score = r2_score(y_true, y_pred)
    mae = mean_absolute_error(y_true, y_pred)
    mse = mean_squared_error(y_true, y_pred)
    total_features.append(len(features)-1)
    total_features_name.append(str(features))
    selected_correlation_value.append(th)
    r2_scores.append(score)
    mae_value.append(mae)
    mse_value.append(mse)

    metrics_dataframe = pd.DataFrame(data= [total_features_name, total_features, selected_correlation_value, r2_scores, mae_value, mse_value],
                                    index = ['features name', '#feature', 'corr_value', 'r2_score', 'MAE', 'MSE'])
    return metrics_dataframe.T

Record the baseline result (threshold = 0.50, three features):

performance_metrics(correlated_data.columns.values, threshold, y_test, y_predict)

The following scatter plots show the linear relationship between each selected feature and house price, with regression lines fitted:

rows = 2
cols = 2
fig, ax = plt.subplots(nrows=rows, ncols=cols, figsize = (16, 4))
ax[0, 0].set_title("House price with respect  to RM")
ax[0, 1].set_title("House price with respect to PTRATIO")
ax[1, 0].set_title("House price with respect to LSTAT")
ax[1, 1].set_title("House price with respect to PRICE")
col = correlated_data.columns
index = 0

for i in range(rows):
    for j in range(cols):
        sns.regplot(x = correlated_data[col[index]], y = correlated_data['Price'], ax = ax[i][j])
        index = index + 1
fig.tight_layout()

The four-panel plot below shows regression lines for each of the three selected features against Price, confirming positive correlation for RM and negative correlations for PTRATIO and LSTAT:

The helper function below retrains the model from scratch for any given correlated feature subset, making it easy to compare different threshold configurations:

corrmat['Price']

CRIM      -0.388305
ZN         0.360445
INDUS     -0.483725
CHAS       0.175260
NOX       -0.427321
RM         0.695360
AGE       -0.376955
DIS        0.249929
RAD       -0.381626
TAX       -0.468536
PTRATIO   -0.507787
B          0.333461
LSTAT     -0.737663
Price      1.000000
Name: Price, dtype: float64

def get_y_predict(corr_data):
    X = corr_data.drop(labels = ['Price'], axis = 1)
    y = corr_data['Price']

    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state = 0)
    model = LinearRegression()
    model.fit(X_train, y_train)
    y_predict = model.predict(X_test)
    return y_predict

Threshold 0.60 — RM and LSTAT only: Tightening the threshold drops PTRATIO and keeps only the two most correlated features:

threshold = 0.60
corr_value = getCorrelatedFeature(corrmat['Price'], threshold)
corr_value

correlated_data = data[corr_value.index]
correlated_data.head()

y_predict = get_y_predict(correlated_data)
performance_metrics(correlated_data.columns.values, threshold, y_test, y_predict)

Dropping PTRATIO and keeping only RM and LSTAT improves R² from 0.49 to 0.54 — fewer but more informative features can outperform a larger set.

Threshold 0.70 — LSTAT only: Raising the threshold further leaves just LSTAT:

corrmat['Price']

CRIM      -0.388305
ZN         0.360445
INDUS     -0.483725
CHAS       0.175260
NOX       -0.427321
RM         0.695360
AGE       -0.376955
DIS        0.249929
RAD       -0.381626
TAX       -0.468536
PTRATIO   -0.507787
B          0.333461
LSTAT     -0.737663
Price      1.000000
Name: Price, dtype: float64

threshold = 0.70
corr_value = getCorrelatedFeature(corrmat['Price'], threshold)
corr_value

correlated_data = data[corr_value.index]
correlated_data.head()

y_predict = get_y_predict(correlated_data)
performance_metrics(correlated_data.columns.values, threshold, y_test, y_predict)

Using LSTAT alone drops R² back to 0.43, confirming that the combination of RM and LSTAT at threshold 0.60 was the sweet spot. Try RM alone as a further comparison:

correlated_data = data[['RM', 'Price']]
correlated_data.head()

y_predict = get_y_predict(correlated_data)
performance_metrics(correlated_data.columns.values, threshold, y_test, y_predict)

Threshold 0.40 — six features: Lowering the threshold to 0.40 brings in more features (INDUS, NOX, RM, TAX, PTRATIO, LSTAT):

threshold = 0.40
corr_value = getCorrelatedFeature(corrmat['Price'], threshold)
corr_value

correlated_data = data[corr_value.index]
correlated_data.head()

y_predict = get_y_predict(correlated_data)
performance_metrics(correlated_data.columns.values, threshold, y_test, y_predict)

Adding weakly correlated features at threshold 0.40 does not improve over the two-feature model, which confirms that feature quality matters more than feature quantity for linear regression.

Normalisation and Standardisation

When features have very different scales (e.g. TAX ranges into the hundreds while NOX is between 0 and 1), models that rely on distances or gradient-based updates can be affected. Standardisation rescales each feature to have zero mean and unit variance (Gaussian distribution):

Normalisation rescales samples to have unit norm and is particularly useful when using kernel-based similarity measures. Scikit-learn provides several ready-made scalers:

The LinearRegression estimator in scikit-learn accepts a normalize parameter that applies feature-wise standardisation before fitting. Try it and observe whether the R² changes:

model = LinearRegression(normalize=True)
model.fit(X_train, y_train)
LinearRegression(normalize=True)
y_predict = model.predict(X_test)
r2_score(y_test, y_predict)

0.48816420156925067

The R² is essentially identical, which is expected for ordinary least squares — linear regression coefficients adjust to scale, so normalisation does not change the fit quality here.

Plotting Learning Curves

Learning curves show how model performance changes as you add more training data. They help you diagnose two common problems: if training and validation scores are both low you have underfitting (high bias); if the training score is high but the validation score is low you have overfitting (high variance).

Import the cross-validation utilities needed:

from sklearn.model_selection import learning_curve, ShuffleSplit

Define a reusable function that plots both the training score and the cross-validation score as a function of training set size:

def plot_learning_curve(estimator, title, X, y, ylim=None, cv=None, n_jobs=None, train_sizes=np.linspace(.1, 1.0, 10)):

    plt.figure()
    plt.title(title)
    plt.xlabel("Training examples")
    plt.ylabel("Score")

    train_sizes, train_scores, test_scores = learning_curve(
        estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes)

    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)

    plt.grid()

    plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std, alpha=0.1,
                     color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
                     test_scores_mean + test_scores_std, alpha=0.1, color="g")

    plt.plot(train_sizes, train_scores_mean, 'o-', color="r",
             label="Training score")
    plt.plot(train_sizes, test_scores_mean, 'o-', color="g",
             label="Cross-validation score")

    plt.legend(loc="best")
    return plt

X = correlated_data.drop(labels = ['Price'], axis = 1)
y = correlated_data['Price']

title = "Learning Curves (Linear Regression) " + str(X.columns.values)

cv = ShuffleSplit(n_splits=100, test_size=0.2, random_state=0)

estimator = LinearRegression()
plot_learning_curve(estimator, title, X, y, ylim=(0.7, 1.01), cv=cv, n_jobs=-1)
plt.show()

Conclusion

In this tutorial you trained a linear regression model on the Boston housing dataset. You explored the data with pairplots and distribution plots, selected features using Pearson correlation thresholds, and compared five different feature configurations. The two-feature model using RM and LSTAT (threshold 0.60) achieved the best R² of 0.54 and MAE of 4.14, outperforming both more restrictive and more relaxed feature sets. You also verified that feature normalisation has no effect on ordinary least squares fit quality, and built reusable learning-curve tooling to diagnose bias and variance.

Key takeaways:

• Linear regression fits a line (or hyperplane) by minimising the mean squared error between predictions and actual values using gradient descent.
• Feature selection by Pearson correlation is a fast and interpretable way to identify which inputs are worth including — but always validate empirically.
• R² measures how much target variance your model explains; MAE and RMSE measure average prediction error in interpretable units.
• More features do not always improve regression — weakly correlated features can add noise rather than signal.
• Learning curves reveal whether your model is underfitting or overfitting and guide decisions about adding more data or regularisation.

Next steps:

• Read Logistic Regression with Python to see how the linear model is adapted for binary classification problems.
• Explore Random Forest Regressor to see how an ensemble of trees can capture non-linear relationships that linear regression misses.
• Experiment with all 13 Boston features (threshold = 0.0) and add Ridge or Lasso regularisation to compare against the correlation-filtered models built here.