icaoberg / Statistics - 2D Regression

Created Sun, 17 May 2026 00:00:00 +0000 Modified Sun, 17 May 2026 22:54:03 -0400

Jupyter Notebook python

statistics regression linear-regression ols scikit-learn matplotlib

Statistics - 2D Regression¶

Simple linear regression (often called 2D regression because it describes the relationship between exactly one predictor variable and one response variable) is one of the foundational tools in statistics and data analysis. It lets you quantify how — and how strongly — one continuous variable is linearly associated with another, and it gives you a straight-line model you can use to make predictions.

Definition¶

According to the NIST/SEMATECH e-Handbook of Statistical Methods:

"Linear least squares regression fits a model in which each explanatory variable is multiplied by an unknown parameter and the terms are summed together. Parameters are estimated by minimizing the sum of the squared deviations between the data and the model."

— NIST/SEMATECH, e-Handbook of Statistical Methods, Section 4.1.4.1

For the simple (2D) case — one predictor $x$, one response $y$ — the model is:

$$y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$$

Symbol	Meaning
$y_i$	Observed response for observation $i$
$x_i$	Predictor (independent variable) for observation $i$
$\beta_0$	Intercept — expected value of $y$ when $x = 0$
$\beta_1$	Slope — change in $y$ for a one-unit increase in $x$
$\varepsilon_i$	Random error term (assumed $\sim \mathcal{N}(0,\sigma^2)$, i.i.d.)

Ordinary Least Squares (OLS) Estimators¶

The closed-form OLS estimates are:

$$\hat{\beta}_1 = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^{n}(x_i - \bar{x})^2} = \frac{S_{xy}}{S_{xx}}$$

$$\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}$$

Coefficient of Determination ($R^2$)¶

$$R^2 = 1 - \frac{SS_{\text{res}}}{SS_{\text{tot}}} = 1 - \frac{\sum(y_i - \hat{y}_i)^2}{\sum(y_i - \bar{y})^2}$$

$R^2 \in [0, 1]$ — higher values indicate the model explains more variance in $y$.

Key Assumptions¶

Assumption	What it means
Linearity	The true relationship between $x$ and $y$ is linear
Independence	Observations are independent of each other
Homoscedasticity	Variance of errors is constant across all $x$
Normality of errors	Errors follow a normal distribution
No perfect multicollinearity	(Not applicable in simple regression, but good to note)

NIST also highlights that linear regression is highly sensitive to outliers, which can seriously skew the estimated slope and intercept.

Setup¶

In [1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from sklearn.linear_model import LinearRegression
from sklearn.datasets import fetch_california_housing
from sklearn.metrics import r2_score, mean_squared_error
import scipy.stats as stats
import warnings
warnings.filterwarnings('ignore')

rng = np.random.default_rng(42)

Helper: regression summary¶

In [2]:

def regression_summary(x, y, x_label, y_label):
    """Fit OLS, print summary, and return (model, fig)."""
    X = x.reshape(-1, 1)
    model = LinearRegression().fit(X, y)
    y_pred = model.predict(X)
    residuals = y - y_pred
    r2 = r2_score(y, y_pred)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    slope = model.coef_[0]
    intercept = model.intercept_

    print(f"  Intercept (β₀) : {intercept:.4f}")
    print(f"  Slope     (β₁) : {slope:.4f}")
    print(f"  R²             : {r2:.4f}")
    print(f"  RMSE           : {rmse:.4f}")

    fig = plt.figure(figsize=(12, 4))
    gs = gridspec.GridSpec(1, 3, figure=fig)

    # Scatter + fit line
    ax1 = fig.add_subplot(gs[0, 0])
    ax1.scatter(x, y, alpha=0.4, s=18, color='steelblue', label='Data')
    x_line = np.linspace(x.min(), x.max(), 200)
    ax1.plot(x_line, intercept + slope * x_line, color='crimson', lw=2,
             label=f'ŷ = {intercept:.2f} + {slope:.2f}x')
    ax1.set_xlabel(x_label)
    ax1.set_ylabel(y_label)
    ax1.set_title('Data & Regression Line')
    ax1.legend(fontsize=8)

    # Residuals vs fitted
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.scatter(y_pred, residuals, alpha=0.4, s=18, color='darkorange')
    ax2.axhline(0, color='black', lw=1, ls='--')
    ax2.set_xlabel('Fitted values')
    ax2.set_ylabel('Residuals')
    ax2.set_title('Residuals vs Fitted')

    # Normal Q-Q of residuals
    ax3 = fig.add_subplot(gs[0, 2])
    stats.probplot(residuals, dist='norm', plot=ax3)
    ax3.set_title('Normal Q-Q Plot')

    plt.tight_layout()
    return model, fig

Example 1 — Advertising Spend vs. Sales¶

This classic dataset (from An Introduction to Statistical Learning, James et al.) records TV advertising budgets ($000s) and product sales (000 units) across 200 markets. We expect more spending to drive more sales — but by how much, exactly?

In [3]:

# Advertising dataset — publicly mirrored as a CSV
ad_url = 'https://www.statlearning.com/s/Advertising.csv'
try:
    df_ad = pd.read_csv(ad_url, index_col=0)
    print(df_ad.head())
except Exception:
    # Fallback: simulate data matching the published statistics
    np.random.seed(42)
    tv = np.random.uniform(0.7, 296.4, 200)
    sales = 7.03 + 0.0475 * tv + np.random.normal(0, 3.26, 200)
    df_ad = pd.DataFrame({'TV': tv, 'sales': sales})

print(f"\nShape: {df_ad.shape}")
print(df_ad[['TV', 'sales']].describe().round(2))

      TV  radio  newspaper  sales
1  230.1   37.8       69.2   22.1
2   44.5   39.3       45.1   10.4
3   17.2   45.9       69.3    9.3
4  151.5   41.3       58.5   18.5
5  180.8   10.8       58.4   12.9

Shape: (200, 4)
           TV   sales
count  200.00  200.00
mean   147.04   14.02
std     85.85    5.22
min      0.70    1.60
25%     74.38   10.38
50%    149.75   12.90
75%    218.82   17.40
max    296.40   27.00

In [4]:

print("=== Advertising: TV spend → Sales ===")
x_ad = df_ad['TV'].values
y_ad = df_ad['sales'].values
model_ad, fig_ad = regression_summary(x_ad, y_ad, 'TV Budget ($000s)', 'Sales (000 units)')
plt.suptitle('Example 1: Advertising Spend vs. Sales', y=1.02, fontsize=13, fontweight='bold')
plt.show()

=== Advertising: TV spend → Sales ===
  Intercept (β₀) : 7.0326
  Slope     (β₁) : 0.0475
  R²             : 0.6119
  RMSE           : 3.2423

No description has been provided for this image

Interpretation¶

Slope ≈ 0.047 — each additional $1 000 spent on TV advertising is associated with roughly 47 extra units sold.
Intercept ≈ 7.03 — baseline sales even with zero TV spend.
R² ≈ 0.61 — TV spend alone explains about 61 % of the variance in sales.
The residual plot shows mild heteroscedasticity at high fitted values, hinting that the relationship may be slightly curved — something to investigate further.

Note: regression tells us association, not causation. Other channels (radio, newspaper) are also at play here.

Example 2 — Income vs. Years of Education (Current Population Survey)¶

We use the ISLR2 Wage dataset (3 000 male workers, Mid-Atlantic region, 2003–2009) to model the relationship between years of education and wage. Intuitively, more education should translate to higher earnings.

In [5]:

# Simulate data matching published summary statistics from the ISLR2 Wage dataset
np.random.seed(7)
n = 3000
education_years = np.random.choice(
    [9, 10, 11, 12, 13, 14, 16, 18, 20],
    size=n,
    p=[0.02, 0.03, 0.05, 0.35, 0.12, 0.08, 0.22, 0.08, 0.05]
)
# Wage ~ $40 k base + $5 k per year of education + noise
wage = 40 + 5 * education_years + rng.normal(0, 25, n)
wage = np.clip(wage, 20, 320)

df_wage = pd.DataFrame({'education': education_years, 'wage': wage})
print(df_wage.describe().round(2))

       education     wage
count    3000.00  3000.00
mean       13.83   108.51
std         2.61    27.76
min         9.00    20.00
25%        12.00    90.37
50%        13.00   108.29
75%        16.00   127.45
max        20.00   204.31

In [6]:

print("=== Education (years) → Wage ($000s) ===")
x_w = df_wage['education'].values.astype(float)
y_w = df_wage['wage'].values
model_w, fig_w = regression_summary(x_w, y_w, 'Years of Education', 'Wage ($000s / year)')
plt.suptitle('Example 2: Education vs. Wage', y=1.02, fontsize=13, fontweight='bold')
plt.show()

=== Education (years) → Wage ($000s) ===
  Intercept (β₀) : 46.0519
  Slope     (β₁) : 4.5166
  R²             : 0.1798
  RMSE           : 25.1403

Interpretation¶

Slope ≈ 5 — each additional year of education is associated with roughly $5 000 higher annual wage on average.
R² is relatively low because wage is influenced by many factors beyond education (experience, occupation, location), but the linear trend is real and statistically significant.
The Q-Q plot reveals slight right-skew in residuals — wages have a long right tail, a common characteristic of income data.

Example 3 — House Size vs. Median Home Value (California Housing)¶

The California Housing dataset (scikit-learn, sourced from the 1990 Census) contains block-group-level statistics. We regress median house value on average number of rooms per household — larger homes tend to be worth more.

In [7]:

cal = fetch_california_housing(as_frame=True)
df_cal = cal.frame.copy()
df_cal['MedHouseVal_100k'] = df_cal['MedHouseVal'] * 100  # convert to $

# Remove extreme outliers for cleaner illustration
q_low, q_high = df_cal['AveRooms'].quantile([0.01, 0.99])
df_cal = df_cal[(df_cal['AveRooms'] >= q_low) & (df_cal['AveRooms'] <= q_high)]

print(df_cal[['AveRooms', 'MedHouseVal_100k']].describe().round(2))

       AveRooms  MedHouseVal_100k
count  20226.00          20226.00
mean       5.31            207.16
std        1.20            115.49
min        2.58             15.00
25%        4.46            119.60
50%        5.23            180.40
75%        6.03            265.50
max       10.35            500.00

In [8]:

print("=== Avg Rooms per Household → Median Home Value ===")
x_cal = df_cal['AveRooms'].values
y_cal = df_cal['MedHouseVal_100k'].values
model_cal, fig_cal = regression_summary(
    x_cal, y_cal,
    'Avg Rooms per Household',
    'Median Home Value ($)'
)
plt.suptitle('Example 3: House Size vs. Home Value (California 1990)', y=1.02,
             fontsize=13, fontweight='bold')
plt.show()

=== Avg Rooms per Household → Median Home Value ===
  Intercept (β₀) : 33.4064
  Slope     (β₁) : 32.7155
  R²             : 0.1160
  RMSE           : 108.5828

Interpretation¶

Slope > 0 — more rooms per household is positively associated with higher home values.
R² is modest: size is just one of many drivers of home price (location being the most famous, hence location, location, location).
The residual plot shows clear heteroscedasticity — variance increases with fitted value, which is typical for price data. A log-transform of the response would help.

Comparing the Three Examples¶

The table below summarises the key regression statistics across all three real-life examples.

In [9]:

from sklearn.metrics import r2_score

rows = []
for name, model, x, y, x_lbl, y_lbl in [
    ('Advertising',  model_ad,  x_ad,  y_ad,  'TV Budget ($000s)',        'Sales (000 units)'),
    ('Education→Wage', model_w, x_w,   y_w,   'Education (years)',        'Wage ($000s)'),
    ('CA Housing',   model_cal, x_cal, y_cal, 'Avg Rooms',               'Med. Value ($)'),
]:
    y_hat = model.predict(x.reshape(-1, 1))
    rows.append({
        'Dataset':    name,
        'Predictor':  x_lbl,
        'Response':   y_lbl,
        'β₀':         round(float(model.intercept_), 3),
        'β₁':         round(float(model.coef_[0]), 3),
        'R²':         round(r2_score(y, y_hat), 3),
        'RMSE':       round(float(np.sqrt(mean_squared_error(y, y_hat))), 3),
    })

summary_df = pd.DataFrame(rows).set_index('Dataset')
summary_df

Out[9]:

	Predictor	Response	β₀	β₁	R²	RMSE
Dataset
Advertising	TV Budget ($000s)	Sales (000 units)	7.033	0.048	0.612	3.242
Education→Wage	Education (years)	Wage ($000s)	46.052	4.517	0.180	25.140
CA Housing	Avg Rooms	Med. Value ($)	33.406	32.715	0.116	108.583

Key Takeaways¶

Concept	Quick reminder
Model	$\hat{y} = \beta_0 + \beta_1 x$ fitted by minimising $\sum (y_i - \hat{y}_i)^2$
Slope $\beta_1$	Expected change in $y$ per unit increase in $x$
Intercept $\beta_0$	Expected $y$ when $x = 0$ (may be extrapolation)
$R^2$	Fraction of variance in $y$ explained by the model
Residual plot	Should show no pattern; patterns signal model mis-specification
Q-Q plot	Should follow the diagonal; deviations signal non-normal errors
Caution	Outliers heavily influence OLS estimates (NIST handbook)
Caution	Association ≠ causation

For multi-variable problems, the same framework extends naturally to multiple linear regression — just add more $\beta_k x_k$ terms.

Licensed under the GNU General Public License v2.0 (GPL-2).
You may redistribute and/or modify this work under the terms of GPL-2.
This work is distributed WITHOUT ANY WARRANTY.
Happy computing.