Linear Regression

Here are the first 8 lines of the advertising data set. Each line shows the sales of a product for given advertising budgets on TV, radio and newspapers.

There are some questions we could ask based on the data.

1. Is there a relationship between advertising budget and sales?
2. How strong is the relationship between advertising budget and sales?
3. Which media are associated with sales?
4. How large is the association between each medium and sales?
5. How accurately can we predict future sales?
6. Is the relationship linear?
7. Is there synergy among the advertising media?

Let's read in the advertising data.

from os import getenv
from pathlib import Path

import polars as pl
from dotenv import load_dotenv

load_dotenv(".envrc")

ISLP_DATADIR = Path(
    getenv("ISLP_DATADIR")
    or (
        "/home/bon/projects/auc/courses/ml/website/"
        "islp/data/ALL CSV FILES - 2nd Edition"
    )
)

ads_csv = ISLP_DATADIR / "Advertising.csv"

ads = pl.read_csv(ads_csv)
print(ads.shape)

(200, 5)

Are sales linear in TV budget? The Seaborn library, built on top of Matplotlib has powerful functionality for plotting regressions. This plot shows not only the regression line but also the confidence level.

import polars_seaborn
import matplotlib.pyplot as plt

plt.style.use("seaborn-v0_8-whitegrid")

ads.sns.lmplot(x="TV", y="sales")
plt.savefig("images/sales-v-TV.png")

How about radio?

ads.sns.lmplot(x="radio", y="sales")
plt.savefig("images/sales-v-radio.png")

And finally newspaper?

ads.sns.lmplot(x="newspaper", y="sales")
plt.savefig("images/sales-v-newspaper.png")

[<Expr ['col("newspaper")'] at 0x7FB581188D50>, <Expr ['col("sales")'] at 0x7FB581188C50>]

We might conclude that sales are quite linear in TV, a little less so in radio and not very much at all in newspaper.

We assume a model of the form

\[ Y = \beta_0 + \beta_1 X + \epsilon \]

where

• \(Y\) is the response or dependent variable
• \(X\) is the predictor or input variable
• \(\beta_0\) and \(\beta_1\) are two numbers representing the intercept and slope
• \(\epsilon\) is an error term

If we have good estimates \(\hat\beta_0\) and \(\hat\beta_1\) we can predict future sales using

\[ \hat{y} = \hat\beta_0 + \hat\beta_1 x \]

When \(\hat{y_i} = \hat\beta_0 + \hat\beta_1 x_i\) is the prediction for \(Y\) for the ith value of \(X\) then we say \(e_i = y_i - \hat{y_i}\) is the ith residual.

The residual sum of squares (RSS) is

We want to choose \(\beta_0\) and \(\beta_1\) to minimize the RSS.

With some calculus we can show the minimizing values are

where \(\bar{y} = \frac{1}{n} \sum_{i=1}^n y_i\) and \(\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\) are the sample means.

We now move to a linear regression that depends on more that one predictor.

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p + \epsilon \]

We interpret each \(\beta_j\) as the average effect on \(Y\) of an increase in \(X_j\) of one unit while holding all other predictor values fixed.

For our advertising data set we have

\[ \mathtt{sales} = \beta_0 + \beta_1 \mathtt{TV}+ \beta_2 \mathtt{radio}+ \beta_3 \mathtt{newspaper} + \epsilon \]

Some caveats:

• When the predictors are uncorrelated, each coefficient can be estimated and tested separately.
• Correlations among predictors cause problems.
• We should not claim causality.

All models are wrong but some models are useful. – George Box

Now we need to estimate \(\hat\beta_0\), \(\hat\beta_1\), \(\hat\beta_2\), …, \(\hat\beta_p\) so that we can predict using

\[ \hat{y} = \hat\beta_0 + \hat\beta_1 x_1 + \hat\beta_2 x_2 + \dots + + \hat\beta_p x_p \]

This works exactly as in the case where we had only one predictor, except that we are now working in \(p+1\) dimensions.

The scikit-learn library provides us with a flexible Linear Regression library. It takes numpy arrays.

import pandas as pd
from sklearn.linear_model import LinearRegression

X = ads[["TV", "radio", "newspaper"]].to_numpy()
y = ads["sales"].to_numpy()

reg = LinearRegression().fit(X, y)

np.set_printoptions(precision=4)

print(
    f"""R^2: {reg.score(X, y):.4f}
Intercept: {reg.intercept_:.4f}
Coefficients: {reg.coef_}
Prediction: {reg.predict([[3, 4, 5]])}"""
)

R^2: 0.8972
Intercept: 2.9389
Coefficients: [ 0.0458  0.1885 -0.001 ]
Prediction: [3.8251]

These are the correlations between the predictors.

print(ads[["TV", "radio", "newspaper", "sales"]].corr())

shape: (4, 4)
┌──────────┬──────────┬───────────┬──────────┐
│ TV       ┆ radio    ┆ newspaper ┆ sales    │
│ ---      ┆ ---      ┆ ---       ┆ ---      │
│ f64      ┆ f64      ┆ f64       ┆ f64      │
╞══════════╪══════════╪═══════════╪══════════╡
│ 1.0      ┆ 0.054809 ┆ 0.056648  ┆ 0.782224 │
│ 0.054809 ┆ 1.0      ┆ 0.354104  ┆ 0.576223 │
│ 0.056648 ┆ 0.354104 ┆ 1.0       ┆ 0.228299 │
│ 0.782224 ┆ 0.576223 ┆ 0.228299  ┆ 1.0      │
└──────────┴──────────┴───────────┴──────────┘

Here are full statistical results for the advertising data as given by the statsmodels library.

import statsmodels.api as sm

estimate = sm.OLS(y, sm.add_constant(X)).fit()

print(estimate.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.897
Model:                            OLS   Adj. R-squared:                  0.896
Method:                 Least Squares   F-statistic:                     570.3
Date:                Thu, 25 Sep 2025   Prob (F-statistic):           1.58e-96
Time:                        07:57:13   Log-Likelihood:                -386.18
No. Observations:                 200   AIC:                             780.4
Df Residuals:                     196   BIC:                             793.6
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.9389      0.312      9.422      0.000       2.324       3.554
x1             0.0458      0.001     32.809      0.000       0.043       0.049
x2             0.1885      0.009     21.893      0.000       0.172       0.206
x3            -0.0010      0.006     -0.177      0.860      -0.013       0.011
==============================================================================
Omnibus:                       60.414   Durbin-Watson:                   2.084
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              151.241
Skew:                          -1.327   Prob(JB):                     1.44e-33
Kurtosis:                       6.332   Cond. No.                         454.
==============================================================================

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

Note:

• The standard errors are low except for the intercept
• The t-statistics are high except for newspaper
• The p-values are low except for newspaper

Linear regression is a parametric approach. The \(\beta\) are the parameters.

Parametric methods are

• easy to fit
• easy to interpret
• easy to test for statistical significance

However they also make strong assumptions about the functional form of \(f\) which may be far from the true relationship.

Non-parametric methods do not make such assumptions and can be more flexible in practice. An example of a non-parametric method we will explore in more detail later is \(K\)-Nearest Neighbours (KNN) regression. We choose a whole-number value for \(K\), say 3. Now to make a prediction for a point \(x_0\) we find the 3 points in our training set which are closest to \(x_0\). We call this set of points the neighbourhood of \(x_0\), written \(\mathcal{N_0}\). We let \(\hat{f}(x_0)\) be the average of the responses corresponding to \(\mathcal{N_0}\).

\[ \hat{f}(x_0) = {1 \over K} \sum_{x_i \in \mathcal{N_0}} y_i \]

The choice of value for \(K\) is a trade-off. Smaller values give a "blocky" function while larger values give a smoother function which can be "over-smoothed". We will study this bias-variance trade-off in more detail later in the course.

So which to choose, a parametric or non-parametric method? Quoting ISLP,

The answer is simple: the parametric approach will outperform the non-parametric approach if the parametric form that has been selected is close to the true form of \(f\).

Let's look at KNN on some almost-linear data.

The black line is the true relationship. The red dots are samples from the true relationship by adding some noise. The blue line is the KNN regression with \(K=1\) on the left and \(K=9\) on the right.

Here are the same data with the least-squares linear regression line in dashed blue. On the right we see the mean squared error for the linear regression in dashed black and for KNN in green for various values of \(K\). Note that the plot is in terms of the inverse \(1/K\).

Because the parametric form we chose is close to the true form of our data, the parametric approach beats the non-parametric approach for all values of \(K\).

Now let's consider data which are not so linear and data that are very non-linear. We see that the non-parametric method begins to shine and indeed outdo the parametric method.

Summing up, we might conclude that KNN is more flexible that linear regression but strange things happen in higher numbers of dimensions. In particular, our notions of distance that we are familiar with in three dimensions or less will begin to let us down. We call this the Curse of Dimensionality and will revisit it later in the course. For now we look at how the MSE behaves for the above example in higher numbers of dimensions \(p\).