Linear Regression

1. Regression problem：Simple linear regression

Slide 1：適合以下問題——

• 只有 $p = 1$ 個 predictor；
• X 同 Y 都係 numerical and continuous；
• X、Y 嘅 scatter plot 理想上似一條直線。

Y	X
Dependent variable / response / predicted target （statsmodels 叫 endogenous variable，endog）	Independent variable / predictor （statsmodels 叫 exogenous variable，exog）

Slides 2–5 例子：美國 48 個州嘅 fuel consumption vs number of driving licenses（dlic；Module 1 slide 8 叫 %Driving License）。Slide 3 嘅 scatter plot 顯示 dlic 越高、fuel 越高，大致線性。

• Problem：根據 X 同 Y 嘅數學關係，由 X 預測 Y。
• Linear model：$Y = a + bX + \text{Error}$（a = intercept，b = slope）。
• Regression：由 training data 估計 a、b 嘅統計方法 → $\hat a,\ \hat b$。新數據只知 X → 預測 $\hat Y = \hat a + \hat b X$。
• Slide 5 嘅 prediction problem：如果新政策將 license 數目減至 50%，執行後 fuel consumption 會係幾多？

2. Problem Solving with Python：OLS

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns
fueldata = pd.read_csv("fueldata.csv")
plt.scatter(fueldata['dlic'], fueldata['fuel'])    # 先睇 scatter plot 係咪線性
plt.show()

Slide 7：Simple linear regression 搵 $\hat a$、$\hat b$ 令 sum of squared errors 最細——OLS (Ordinary Least Squares)：

model = smf.ols('fuel~dlic', data=fueldata).fit()   # 'Y~X'：左邊 dependent，右邊 independent
print(model.summary())
print(model.params)                                  # 睇 a-hat 同 b-hat

$\hat a$ 讀做 “a-hat”，代表由數據估計出嚟嘅值。

2.1 讀 OLS Regression Results（slide 8）

項目	值	意思
Dep. Variable	fuel	Y
No. Observations	48	n
Df Residuals / Df Model	46 / 1	$n - 2$ / $p$
R-squared / Adj. R-squared	0.489 / 0.477	dlic 解釋到 fuel 約 48.9% 嘅 variation
F-statistic / Prob (F)	43.94 / 3.29e-08	整個 model 顯著（ANOVA F-test）
AIC / BIC	559.9 / 563.6	model selection 用，越低越好

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	−227.3091	121.862	−1.865	0.069	−472.604	17.986
dlic	14.0984	2.127	6.629	0.000	9.817	18.380

2.2 畫 fitted line 同 prediction（slides 9–10）

plt.scatter(fueldata['dlic'], fueldata['fuel'])
predicted_values = model.predict()                 # training data 嘅 fitted values
plt.plot(fueldata['dlic'], predicted_values)
plt.show()

newdata = pd.DataFrame({'dlic': [50]})             # test dataset：column 名要同 formula 一樣
prediction = model.predict(newdata)                # 477.61
newdata = pd.DataFrame({'dlic': [50, 55, 60]})
predictions = model.predict(newdata)               # 477.61, 548.10, 618.59

3. Step-by-step numerical example：Matrix OLS（slides 11–19）

3 個數據點：(1, 2)、(2, 3)、(3, 5)。Slide 12 嘅圖顯示 OLS line 穿過三點附近。

3.1 Derivation（slides 13–16）

Objective function（slide 14）：選 $\hat\beta$ 令 Sum of Squared Residuals 最細；column vector 嘅平方和 = inner product $e^Te$：

（因為 $y^TX\hat\beta$ 係 scalar，等於自己嘅 transpose $\hat\beta^TX^Ty$，兩個中間項可以合併。）

Matrix differentiation & first-order conditions（slide 15）：兩條 matrix calculus rules——Linear rule $\partial(A^T\beta)/\partial\beta = A$；Quadratic rule（A symmetric）$\partial(\beta^TA\beta)/\partial\beta = 2A\beta$：

Normal Equations（slide 16）：

前提：design matrix X 有 full column rank（冇 perfect multicollinearity），咁 $X^TX$ 先可逆。

3.2 手計（slides 17–18）

3.3 Python code（slide 19）

import numpy as np
import statsmodels.api as sm
x_data = np.array([1, 2, 3])
Y = np.array([2, 3, 5]).reshape(-1, 1)
X = np.column_stack((np.ones_like(x_data), x_data))   # design matrix：加一欄 1
XT = X.T
XTX = XT @ X
XTX_inv = np.linalg.inv(XTX)
XTY = XT @ Y
beta_hat = XTX_inv @ XTY
print(f"Intercept (Beta 0): {beta_hat[0][0]:.3f}")      # 0.333
print(f"Slope (Beta 1):     {beta_hat[1][0]:.3f}")      # 1.500
# 同 library 比較：endog = dependent，exog = independent（要自己加 constant 欄）
model = sm.OLS(endog=Y, exog=X).fit()
print(model.summary())

4. R squared（slide 20）

$R^2$ 係 0 至 1 之間嘅數，描述 Y 嘅 variation 有幾多比例被 model 解釋。

Explained / Regression SS (ESS)：$ESS = TSS - RSS$。$R^2 = 1$ 代表 perfect fit（$RSS = 0$）。

5. When does regression fail?（slides 21–22）

5.1 Visual diagnostics check（slides 23–25）

resid = model.resid
fitted = model.fittedvalues
n = len(resid)
fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].scatter(fitted, resid);            axs[0, 0].axhline(0, color='red', lw=2)  # e(i) vs fit(i)
axs[0, 1].scatter(np.arange(n), resid);      axs[0, 1].axhline(0, color='red', lw=2)  # e(i) vs i
axs[1, 0].scatter(resid[:-1], resid[1:]);    axs[1, 0].axhline(0, color='red', lw=2)  # e(i) vs e(i-1)
sm.qqplot(resid, ax=axs[1, 1]);              axs[1, 1].set_title('Normal Q-Q')
plt.show()
plt.hist(resid, bins=20); plt.show()         # residual histogram

Slide 24 係 fuel model 嘅四幅 diagnostic plots：residual vs predicted fuel、vs index、vs lagged error 大致都係圍繞 0 嘅散點，冇明顯 pattern；QQ-plot 大致沿直線，但右上尾有幾個大 residual（約 180、250）偏離 → 可能有 outliers／少少 right skew。（Summary 嘅 Skew = 0.747、Prob(Omnibus) = 0.049 亦提示 residual 唔算完全 normal；Durbin-Watson = 1.645 接近 2，autocorrelation 唔嚴重。）

6. Multiple linear regression（slides 26–29）

Applicable if：

• 所有 X 同 Y 都係 numerical and continuous；
• p 相對 n 細（$p > 10$ can be very problematic）。

Lecture 4 第二部分內容：basic ideas of regression for $p > 1$、doing prediction using Python、selecting important predictors。

例子（slide 29）：What affects the saving rate？savings.csv：50 個國家 1960–1970 數據（Belsley, Kuh & Welsch 1980, Regression Diagnostics）。

Variable	意思
$Y = \text{sr}$	personal saving ÷ disposable income
$X_1 = \text{pop15}$	% population under age 15
$X_2 = \text{pop75}$	% population over age 75
$X_3 = \text{dpi}$	per capita disposable income (US$)
$X_4 = \text{ddpi}$	% growth rate of dpi

7. Multiple regression with Python（slides 30–34）

要識：點寫 regression formula、計 coefficients、指定 training / test data。

model = smf.ols('sr~pop75+dpi', data=savings).fit()   # 多個 X 用 + 連接
print(model.summary())                                 # β-hat 喺 "coef" 欄

Slide 32 嘅 output：

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	7.0566	1.290	5.468	0.000	4.461	9.653
pop75	1.3050	0.778	1.678	0.100	−0.259	2.869
dpi	−0.0003	0.001	−0.337	0.738	−0.002	0.002

$R^2 = 0.102$、Adj. $R^2 = 0.064$、$F = 2.680$、Prob (F) = 0.0791、$n = 50$、Df Residuals = 47、Df Model = 2、AIC = 291.5、BIC = 297.2。

7.1 Matrix notation（slide 33）

Gradient = $-2X'Y + 2(X'X)\boldsymbol\beta$，設 = 0：

第一欄全 1 係「for intercept」。Simple regression 只係 p = 1 嘅特例。

7.2 Prediction for new data（slide 34）

newdata = pd.DataFrame({'pop15': [30], 'pop75': [3], 'dpi': [2000], 'ddpi': [5]})
newdata = pd.DataFrame({'pop15': [30, 25], 'pop75': [3, 4],
                        'dpi': [2000, 2000], 'ddpi': [5, 5]})
predictions = model.get_prediction(newdata)   # 有 summary_frame，可出 interval
predictions = model.predict(newdata)          # point prediction

例：用 slide 32 model，$\text{pop75} = 3$、$\text{dpi} = 2000$：$\widehat{sr} = 7.0566 + 1.305(3) - 0.0003(2000) = 7.0566 + 3.915 - 0.6 = 10.37$（用四捨五入後嘅 coef；多出嚟嘅 columns 如 pop15、ddpi 會被 model 忽略）。

8. Measuring prediction ability（slides 35–39）

savings_train = savings.iloc[:35]        # 頭 35 個國家做 training（70%）
savings_test = savings.iloc[35:50]       # 其餘 15 個做 test（30%）
model = smf.ols('sr~pop15+ddpi', data=savings_train).fit()
from sklearn.metrics import mean_squared_error
pre_val = model.predict(savings_test)
act_val = savings_test['sr']
bias = np.mean(pre_val - act_val)        # 正 = 系統性高估；負 = 系統性低估
rmse = np.sqrt(mean_squared_error(act_val, pre_val))
print("Bias:", bias)
print("RMSE:", rmse)                     # 亦可以用 MSE

8.1 Big model vs small model（slides 36–39）

# Small model vs Big model（slide 39）
small_model = smf.ols('sr~pop15+ddpi', data=savings_train).fit()
big_model = smf.ols('sr~pop15+pop75+dpi+ddpi', data=savings_train).fit()
for name, m in [("small", small_model), ("big", big_model)]:
    in_rmse = np.sqrt(mean_squared_error(savings_train['sr'], m.predict(savings_train)))
    out_rmse = np.sqrt(mean_squared_error(savings_test['sr'], m.predict(savings_test)))
    print(f"{name}: AIC={m.aic:.1f}, BIC={m.bic:.1f}, In-RMSE={in_rmse:.4f}, Out-RMSE={out_rmse:.4f}")

9. Optional topics（slides 40–46）

9.1 Interval prediction（slide 41）

newdata = pd.DataFrame({'dlic': [50]})          # simple regression 例子
predictions = model.get_prediction(newdata)
prediction_summary = predictions.summary_frame(alpha=0.05)
print(prediction_summary)
# mean_ci_lower/upper：population mean（平均 Y）嘅 CI
# obs_ci_lower/upper ：future prediction / new data point 嘅 prediction interval（更闊）

9.2 Is X helpful to predict Y?（slides 42–43）

	ANOVA (F-statistic)	T-tests
$H_0$	所有 regression coefficients（intercept 除外）= 0：$\beta_1=\beta_2=\cdots=\beta_p=0$	$\beta_j = 0$：第 j 個 predictor 對 Y 冇 effect（given other variables）
$H_1$	至少一個 $\ne 0$	$\beta_j \ne 0$：$X_j$ 同 Y 有顯著線性關係
測試咩	整個 model 係咪比 intercept-only mean model 解釋顯著更多 variance	個別 predictor 係咪顯著

Slide 44 係 F-table（$\alpha = 0.05$）：column = Numerator DF（1–10），row = Denominator DF（1–50）。例如 $F_{1,1} = 161.448$、$F_{1,3} = 10.128$、$F_{1,46} = 4.052$、$F_{2,47} = 3.195$。

9.3 Is the coefficient significantly different from zero?（slides 45–46）

統計估計同 prediction 通常會附上 standard error 量度準確度／不確定性。

Slide 46 係 t-table（同 Module 3a 一樣，row $\nu$、column $\alpha$）。

自我測驗 Self-test

Q1. 用 matrix 方法計 (0, 1), (1, 3), (2, 2) 嘅 OLS line。

$X^TX = \begin{bmatrix}3&3\\3&5\end{bmatrix}$，det = 15 − 9 = 6，$(X^TX)^{-1} = \frac16\begin{bmatrix}5&-3\\-3&3\end{bmatrix}$；$X^TY = \begin{bmatrix}6\\7\end{bmatrix}$；$\hat\beta = \frac16\begin{bmatrix}30-21\\-18+21\end{bmatrix} = \begin{bmatrix}1.5\\0.5\end{bmatrix}$ → $\hat Y = 1.5 + 0.5X$。

Q2. 繼續 Q1：計 TSS、RSS、$R^2$。

$\bar Y = 2$；$TSS = 1 + 1 + 0 = 2$。Fitted = 1.5, 2, 2.5；residuals = −0.5, 1, −0.5；$RSS = 0.25 + 1 + 0.25 = 1.5$。$R^2 = 1 - 1.5/2 = 0.25$。

Q3. 某 output：coef of X = −2.35，std err = 0.80，$n = 25$（simple regression）。做 two-sided t-test（$\alpha = 0.05$）並解讀 coefficient。

$t = -2.35/0.80 = -2.94$；$df = 23$，$t_{23,0.975} = 2.069$。|−2.94| > 2.069 → reject $H_0: \beta_1 = 0$。解讀：X 每增加 1 單位，Y 平均減少 2.35 單位。

Q4. Fuel model：$\text{dlic} = 65$ 時嘅預測？如果政策令 dlic 由 60 跌到 30，預測 fuel 變化？有咩問題？

$\text{dlic} = 65$：−227.3091 + 14.0984 × 65 = 689.09。60 → 30：變化 = 14.0984 × (−30) = −422.95。問題：$\text{dlic} = 30$ 可能遠超出 training data 範圍（extrapolation），linear 關係未必成立；而且 regression 只顯示 association，唔一定係 causation。

Q5. Residual vs fitted plot 呈漏斗形（fitted 越大，residual 越散），代表咩問題？違反咩 assumption？

Heteroskedasticity（non-constant variance），違反 identically distributed（constant Var(error)）assumption。可考慮 transformation（例如 log Y）或 GARCH 類 time-varying variance model。

Q6. Simple regression，$n = 22$，$SSR = 120$，$SSE = 200$。完成 ANOVA table 並做 F-test（$\alpha = 0.05$）。$R^2$ 係幾多？

df：Regression 1、Error 20、Total 21。$MSR = 120$；$MSE = 200/20 = 10$；$F = 12$。$F_{1,20,0.05} = 4.351$ → 12 > 4.351 → reject $H_0: \beta_1 = 0$。$SST = 320$，$R^2 = 120/320 = 0.375$。

Q7. 點解 big model 嘅 in-sample RMSE 一定 ≤ small model，但 out-of-sample 唔一定？點樣揀？

OLS 喺 training data 上最小化 SS；多咗 parameters 最差都可以將佢哋設 0 得返 small model，所以 in-sample error 唔會更大。但多餘 parameters 會 fit 埋 noise（overfitting），喺 unseen test data 表現可能更差。揀 model：比較 out-of-sample RMSE，或 AIC / BIC（越低越好）。

Q8. 寫 Python：用 savings 頭 35 行 fit sr ~ pop15 + pop75，計 test set 嘅 bias 同 RMSE。

m = smf.ols('sr~pop15+pop75', data=savings.iloc[:35]).fit(); te = savings.iloc[35:50]; pred = m.predict(te); bias = np.mean(pred - te['sr']); rmse = np.sqrt(mean_squared_error(te['sr'], pred))

Q9. Savings output 中 pop75 coef = 1.305、$p = 0.100$。用完整句子解讀，並講喺 $\alpha = 0.05 / 0.10$ 嘅結論。

喺 dpi 不變下，75 歲以上人口比例每增加 1 個百分點，saving rate 平均增加 1.305 個百分點。$\alpha = 0.05$：$p = 0.100 > 0.05$ → 唔顯著；$\alpha = 0.10$：$p = 0.100$ 啱啱喺邊界（p 唔細過 0.10），嚴格嚟講都唔 reject。

一頁總結 Cheat sheet