Adjusted R-Squared 小整理

Adjusted R2 是 R2 的進階版, 它考慮了模型中變數 (預測因子) 的數量, 懲罰那些濫竽充數的變數, 在多變量回歸 (multiple regression ) 中達到去蕪存菁的效果.

聽起來相當美好.  但我們怎麼知道何時需要它? 要是我們的 model 已經是最好了, 再調整它是否浪費時間? 我們可以先看 R2 的破綻在哪裡?

R2 = 1 – SSR/SST

若 SSR/SST < 1, 表示這是個真分數. 當我新增一個因子, 表現得和平均 (猜的期望值) 一樣好, 真分數 (分子分母加同樣數字) 愈加愈大, R 會愈來愈小. 所以加到沒用的因子, R 會變差 (Adjusted R 變差). 反過來, 如果我們是要減少變數, 也可以預測減掉一個變數是否有感?

Adjusted R² = 1 – [(1 – R²) * (n – 1) / (n – k – 1)]

  • 如果加入的變數有貢獻 (提高解釋力), Adjusted R² 會上升.
  • 如果加入的變數無貢獻 (p-value 高, 不顯著), Adjusted R² 會下降或持平.

基本上, k 如果增加一個變數, 但它對 R2 又無額外貢獻, 分母的 (n – k – 1) 變小, … 以小學程度的數學可知 Adjusted R² 變小.

順便偷渡一下 n. n -> ∞, (n – 1) / (n – k – 1) -> 1, Adjusted R² 越接近 R² (增加變數的懲罰很小). n 越小, 分母越小, 懲罰項越大, 讓 Adjusted R² 更容易小於 R², 這能更敏感地偵測冗餘變數 (overfitting). 但因為它效力強到難以分辨 feature 好壞, 所以一般不調整 n, 只調整 k.

舉例而言, 現在要估計房價模型. 變數有坪數、屋齡、噪音 (亂數) 三個, k = 3. 統計樣本 n = 100 個. 為了說明方便, 我們先在 Python 固定亂數種子. 做出一堆可以收斂的假資料.

import statsmodels.api as sm
import pandas as pd
import numpy as np

# 步驟 1: 固定亂數種子與資料(同之前)
np.random.seed(0)
n = 100
X = pd.DataFrame({
    '坪數': np.random.uniform(20, 200, n),
    '屋齡': np.random.uniform(1, 50, n),
    '噪音': np.random.normal(0, 1, n)
})
y = 5000 * X['坪數'] - 2000 * X['屋齡'] + np.random.normal(0, 10000, n)  # real y

安裝 library ( pip install statsmodels) 後, 很多東西都不用算. 包括 OLS (Ordinary Least Squares) model 都會建好. 因為太方便了, 所以定義一個會印出過程的 function.

# 函數:計算 R²/adj_r2、係數和 p-value
def compute_and_compare(model_vars, model_name):
    X_sub = sm.add_constant(X[model_vars])  # 加常數項
    model = sm.OLS(y, X_sub).fit()          # 擬合 OLS 模型
    
    # 計算 R² 和 adj_r2
    r2 = model.rsquared
    adj_r2 = model.rsquared_adj
    print(f"\n{model_name} (變數: {model_vars}):")
    print(f"  R² = {r2:.4f}")
    print(f"  Adjusted R² = {adj_r2:.4f}")
    
    # 得到估計係數
    params = model.params
    print("  估計係數(params):")
    for var in params.index:
        print(f"    {var}: {params[var]:.2f}")
    
    # 得到每個變數的 p-value
    pvalues = model.pvalues
    print("  每個變數的 p-value(顯著性檢定):")
    for var in pvalues.index:
        print(f"    {var}: {pvalues[var]:.4f} (顯著: {pvalues[var] < 0.05})")
    
    print()  # 空行分隔

讓 3 個變數的 model A 跑一下, 然後讓除掉噪音 (random noise, 不是房子旁邊有噪音), 用只有 2 個變數的 model B 跑一下.

# 步驟 2: 計算模型 A (有噪音)
compute_and_compare(['坪數', '屋齡', '噪音'], "模型 A (有噪音)")

# 步驟 3: 計算模型 B (無噪音)
compute_and_compare(['坪數', '屋齡'], "模型 B (無噪音)")

執行上面程式, 會印出:

模型 A (有噪音) (變數: ['坪數', '屋齡', '噪音']):
  R² = 0.9745
  Adjusted R² = 0.9736
  估計係數(params):
    const: -949.47
    坪數: 4995.64
    屋齡: -2003.71
    噪音: 42.60
  每個變數的 p-value(顯著性檢定):
    const: 0.9644 (顯著: False)
    坪數: 0.0000 (顯著: True)
    屋齡: 0.0000 (顯著: True)
    噪音: 0.9655 (顯著: False)


模型 B (無噪音) (變數: ['坪數', '屋齡']):
  R² = 0.9744
  Adjusted R² = 0.9737
  估計係數(params):
    const: -947.35
    坪數: 4995.64
    屋齡: -2003.72
  每個變數的 p-value(顯著性檢定):
    const: 0.9623 (顯著: False)
    坪數: 0.0000 (顯著: True)
    屋齡: 0.0000 (顯著: True)

我們先關心 Model A 的 R² = 0.9745 > Adjusted R² = 0.9736, 或曰 Adjusted R² < , 這個訊號是提示我們的 model 可能有冗餘的變數存在! 調一調會更好! 至於誰是老鼠屎? 這時還看不出來.

但我們也印出了 p-value, 噪音的 p-value = 0.9655 很不顯著, 所以就抓它了! 在 Model B, 我們只用兩個變數.

此時 R² = 0.9744 > Adjusted R² = 0.9737 仍然成立, 或曰 Adjusted R² < 就是成立, 那這個時候戰犯要抓誰呢?

const 這個 p-value 值雖然大了點, 但是通常不會是標的. 而且坪數和屋齡的 p-value 都為 0 了, 表示極為顯著! 基本上就沒得挑了. 這也告訴我們指標並非絕對的, 要綜合起來使用.

Model B 的 R² 雖然輸給 Model A (0.9744 < 0.9745), 但我們本來就是抓 overfitting (模型在訓練資料上表現好, 但在未見過的資料上表現差), 所以不會太在意這裡掉了一點. 反而是看 Adjusted R² 從 0.9736 升到 0.9737, 還拔掉一個 p-value 很大的變數, 故改用 model B 是值得肯定的方向.

R-Squared 小註解

原本寫完兩篇 R, 就要來寫 R2, 我覺得這個哏不錯! 不過呢, 因為寫著寫著歪樓了, 所以把那邊的主題讓給 P-value, R2 放到後面來.

前面說到 P-value 用來衡量變數的顯著性. 由於我舉的是單一變數的 t-test, 好像 p-value 就足以解釋整個 model. 但真正能解釋整個 model 的是 R2 , 說明如下. 假設:

全部的誤差 SST = SSR + SSE .

SSR = residual variation = model 無法解釋的 (實際觀測值 – 模型預測值)

SSE = regression model variation = model 解釋掉的部分 (模型預測值 – y 的平均值)

R2 代表模型解釋變異的能力, 所以是 model 可以解釋的變異的比重 = SSE / SST = 1 – SSR / SST. 所以怎麼看 R 都是介於 0~1 之間. R 越接近 1, 表示 model 解釋得更好.

換個生活化的例子來說明.

  • 假設抽樣調查了三對父子, 記下父親與兒子的身高 (單位:公分):
父親身高 (x)兒子身高 (y)
170172
175177
180179
  • 新手根據網路謠言, 建構了下面這個回歸模型:
  • 我們來檢驗看看. 首先將每個父親的身高 x 帶入計算, 預測兒子身高 y^.
父親身高 (x)兒子實際身高 (y)預測兒子身高 y^
170172177.4
175177178
180179178.6
  • 計算 R2
  1. 計算 yˉ​= (172+177+179) /3 ​=528/3​ ≈ 176 — 兒子的平均身高
  2. 計算 SST (Sum of squared Total) = (172−176)2 + (177−176)2 + (179−176)2 = (−4)2+(1)2+(3)2=16+1+9=26 — 兒子的身高和平均身高的差的平方和
  3. 計算 SSR (Sum of Squared Residual) = (172-177.4)2 + (177-178)2 + (179-178.6)2 = 29.16+1+0.16=30.32 — 兒子的身高和預測身高的差的平方和

R2 = 1 – SST / SSR = 1 – 30.32 /26 = 1 – 1.166 ≈ -0.166

因為 R2 < 0, 或 SST > SSR, 等於是說: 全部猜平均身高的話, 還比用 model 估來得準. 所以此 model 很不好!

如果我們一項一項分開算. SSE = (177.4−176)2 + (178−176)2 + (178.6−176)2 = 1+4+6.76 = 12.72

SST ≠ SSR + SSE 或 26 ≠30.32 + 12.72

為何會有這樣的結果呢? 主要是這個 model 不正確, 所以什麼都解釋不了. 發生這種問題可能的原因包括:

  • 模型套錯:使用的模型和資料的型態根本不吻合(例如線性模型去解釋非線性資料)
  • 資料點太少:樣本數極少,容易受到極端值影響
  • 預測變數選錯:所用的 x、公式、或回歸參數明顯偏離資料特性
  • 資料雜訊過高:資料本身的變異無法用簡單模型捕捉

讓我們只用這三筆資料來生成合理的 linear regression model, 並重新計算 R2. 這個部分就請 AI 來執行了.

1. 原始資料
父親身高 (xx)兒子身高 (yy)
170172
175177
180179
2. 計算步驟
步驟一:計算平均值
步驟二:計算斜率 (b)


b=0.7

步驟三:計算截距 (a)

a=53.5

3. 生成線性回歸方程式

​解釋為: 父親每增加 1 公分身高, 兒子的預測身高增加 0.7 公分.

4 . 重新計算 R2

這個新模型不就超準了嗎? 但我們這裡只是建立 model, 並且評估其效果. 還不可以把它視為真正能實用. 下圖可以看到網路謠言版和依照數學原理建立的 model 有何差異? 網路謠言 model 不能解釋我們的資料庫, 但實測會不會更好? 這件事無法反映得出來.

T-Test 和 P-Value

雖然今天有望創下投資市值新高 (舊紀錄是 2025/2/20) – 台股平平、日股創高、美股小漲還沒收市, 結果還不好說 . 但這篇的重點是整理上班會用到的一門統計技術.

假設我有大把的 features, 要怎麼知道每個因子對 model 有沒有影響? 以及因子之間是否過於相關 (例如華氏度和攝氏度) 而顯得多餘呢? 讓我們從基本問題開始理解.

首先我們會遇到每個因子的 scale 差異很大. 例如預測房價時, 我們可能有屋齡、坪數、坪單價、附近實價登錄總價等因子. 前兩個算它們從 0~1,000 應該就夠用. 甚至普通人應該不會買千坪豪宅、千年老宅 (古剎?). 至於後兩個就有趣了. 假設我們用萬元為單位, 那麼 後兩者的 order 可能差個 10 倍到上百倍 (10~100 坪), 更不用說和前兩個因子擺在一起比較, 總會有些因子顯得參差不齊~~~

甚至那些 one hot encoding 來的, 根本就是單位向量, 這時要怎麼辦呢? 我們可以對 skewed 的變數取對數 log, 或用標準化 (z-score) 將所有因子拉到類似範圍。這類轉換能處理加法和乘法效應的混合 [1],讓模型更公平. 這部分先帶過, 減少數學、訴求直覺就好.

接下來要篩選出有用的因子, 把雜質打掉, 才能真正抓到關聯性. 我們要先檢查因子之間的相關性 (避免如華氏/攝氏度的多餘), 並確保數據的有效性. 如果數據本身沒意義, 其他動作都是多餘的.

以下就介紹 t-test 和 p-value:它們是用來測試每個因子是否真的對目標變數 (如房價) 有顯著影響. 簡單說, t-test 測量因子的 “效果強度", p-value 則告訴你這種效果在隨機情況下發生的機率有多高— p 小表示影響顯著!

顯然, 我們要做一些測試才會產生測量因子. 假設我們用了單樣本測試 (One-sample t-test), t 公式表示如下 (t 公式還有很多型態).

例如測試新藥是否使血壓平均降至 120 mmHg, 多筆測量出的血壓為 x, 假設值 μ0 為吃了藥但無效果的血壓 (虛無假設), 此時可以看出 x 的平均值 xˉ 並不是 μ0. 假設此藥對某人藥效不佳, 吃了高血壓藥, 收縮壓從 140 變成 150, 差異是 +10, 或是對某人藥效太強, 收縮壓從 135 變成 90, 差異就是 -45.

  • 假設資料(血壓變化值,單位 mmHg):
    • 受試者 1: -5(血壓下降 5)
    • 受試者 2: -3
    • 受試者 3: -7
    • 受試者 4: -4
    • 受試者 5: -6
    • 樣本平均 xˉ=−5xˉ=−5(整體平均下降 5 mmHg)
    • 樣本標準差 s ≈ 1.58
    • 樣本大小 n = 5
    • 虛無假設 μ₀ = 0(藥物無效果)

此時 t = (-5-0)/(1.58 / √5)=-5/​0.707≈−7.07,

t 可以表示出藥效過高或試過低. 由於 t 就是一個正規化的過程. 所以 x 不需要預先做正規化處理.

另外還有個 p-value, 它的計算不是公式, 而是機率的累積函數 CDF (cumulative distribution function). 再以血壓藥為例, p 是紅線以下的面積 + 右邊對稱位置的面積 (所以 x2).

p=2×(1−CDFt​(∣−7.07∣,df=4)) = 0.0021

也就是在自由度 [2] df = n – 1 = 5 – 1 = 4 的情況下, 落在分布長尾部分的機率 (tail probability).

用實際的數字翻譯上述公式就是: t ≈ -7.07, df=4, p ≈ 0.0021 解釋為: 如果藥物無效(H0:μ=0),則發生五個受試者 t = -7.07 或更極端 (更負或更正) 的機率只有 0.2%. 隱含這種事不太可能發生, 藥物相當可能是有效的.

換個角度, 如果大家試完藥, t -> 0, 測起來像是這個藥吃了沒效. 但 p -> 1 只是表示這個測試的證據力很低, 表示無法證明這個藥有效, 而不是證明這個藥非常無效.

[REF]

  1. http://ch.whu.edu.cn/cn/article/doi/10.13203/j.whugis20210659?viewType=SUP
  2. 自由度 – 在計算樣本標準差(s)時,我們先用樣本平均值 (xˉ) 來「約束」資料: 一旦知道平均值, 前 n-1 個值可以自由變化, 但最後一個值必須調整以符合平均值. 這 “犧牲" 了 1 個自由度. 當自由度愈低, 表示愈不確定, 長尾就比較肥. 對同一個 t, p-value 相對大.

R 小筆記 – 2

R 真是不一般, 首先它不像別的語言, index 是從 0 開始. 而是從 1 開始. 再來呢, 很多地方的 1 看起來是 1 但不是 1; 看起來是 0 但也不是 0; 2 也不是平常的 2. 底下舉例說明.

上次說到用 X predict Y, 像是身高預測. 因為最低身高不是 0, 這條 minimum error 的線是有截距的. 此時一般線性迴歸的斜率公式是:

  • 模型:y = α + βx
  • 斜率:β = r(x,y) × sd(y)/sd(x)

假設強制要求通過原點, 他的迴歸如下.

  • 模型:y = βx(無截距項)
  • 斜率:β = Σ(xy)/Σ(x²)

最後, 假設原本是有截距的迴歸, 但我把它的數據先做正規化, 那它會跟誰一樣嗎?

x <- c(1.2, 2.8, 3.5, 1.9, 4.1, 2.3, 3.8, 1.6, 2.9, 3.2)
y <- c(2.4, 5.1, 6.8, 3.2, 7.9, 4.1, 7.2, 2.8, 5.5, 6.1)

# 標準化(z-score normalization)
x_norm <- (x - mean(x)) / sd(x)
y_norm <- (y - mean(y)) / sd(y)

# 方法0:一般線性迴歸斜率公式
slope0 <- cor(x, y) * sd(y) / sd(x)

# 方法1:通過原點迴歸
slope1 <- sum(x * y) / sum(x^2)

# 方法2:先正規化後當作一般線性迴歸
slope2 <- cor(x_norm, y_norm) * sd(y_norm) / sd(x_norm) 

cat("有截距 - 不通過原點迴歸:", slope0, "\n")
cat("通過原點迴歸:", slope1, "\n")
cat("標準化後做一般回歸:", slope2, "\n")

結果我們會得到 slope0 ≈ 1.85, slope1 ≈ 1.94, slope2 ≈ 0.93, 三個都不一樣.

方法 0 和方法 2 當中, 因為按照定義 :

sd(x_norm) = sd(y_norm) = 1, mean = 0

而 correlation 方面: cor(x,y) 又等於 cor(x_norm, y_norm),

所以 sd(y)/sd(x) ≠ 1 的情況下, 會導致 slope0 ≠ slope2.

方法 1 是強制通過原點, 而原本最佳解可能是有截距的, 所以 slope0 也不太容易巧合和 slope1 、slope2 一樣.

上面講的還沒有發揮到 R 的精神, 我們用 lm() 重做一次.

x <- c(1.2, 2.8, 3.5, 1.9, 4.1, 2.3, 3.8, 1.6, 2.9, 3.2)
y <- c(2.4, 5.1, 6.8, 3.2, 7.9, 4.1, 7.2, 2.8, 5.5, 6.1)

# 標準化
x_norm <- (x - mean(x)) / sd(x)
y_norm <- (y - mean(y)) / sd(y)

# === 手動計算(原始方法) ===
slope0_manual <- cor(x, y) * sd(y) / sd(x)
slope1_manual <- sum(x * y) / sum(x^2)
slope2_manual <- cor(x_norm, y_norm) * sd(y_norm) / sd(x_norm)

# === 使用 lm() 實現 ===
# 方法0:一般線性迴歸(有截距)
model0 <- lm(y ~ x)
slope0_lm <- coef(model0)[2]  # 取斜率係數

# 方法1:通過原點迴歸(無截距)
model1 <- lm(y ~ x - 1)  # -1 表示移除截距項
slope1_lm <- coef(model1)[1]  # 只有一個係數

# 方法2:標準化後的一般迴歸
model2 <- lm(y_norm ~ x_norm)
slope2_lm <- coef(model2)[2]  # 取斜率係數

# === 比較結果 ===
cat("=== 手動計算 vs lm() 比較 ===\n")
cat("方法0 - 有截距迴歸:\n")
cat("  手動計算:", slope0_manual, "\n")
cat("  lm() 結果:", slope0_lm, "\n")
cat("  差異:", abs(slope0_manual - slope0_lm), "\n\n")

cat("方法1 - 通過原點迴歸:\n")
cat("  手動計算:", slope1_manual, "\n")
cat("  lm() 結果:", slope1_lm, "\n")
cat("  差異:", abs(slope1_manual - slope1_lm), "\n\n")

cat("方法2 - 標準化後迴歸:\n")
cat("  手動計算:", slope2_manual, "\n")
cat("  lm() 結果:", slope2_lm, "\n")
cat("  差異:", abs(slope2_manual - slope2_lm), "\n\n")

# === 顯示完整模型資訊 ===
cat("=== 完整模型摘要 ===\n")
cat("\n--- 方法0:一般線性迴歸 ---\n")
print(summary(model0))

cat("\n--- 方法1:通過原點迴歸 ---\n")
print(summary(model1))

cat("\n--- 方法2:標準化後迴歸 ---\n")
print(summary(model2))

我們可以看到 model1 <- lm(y ~ x – 1) , 它不是真的減 1, -1 表示移除截距項. 還可以表示為 y ~ 0 + x. 兩者等效.

再來 slope2_lm <- coef(model2)[2] 的 [2] 不是第三個係數. 其中, 第一個係數 [1] 就固定是截距, 第二個係數 [2] 就固定是斜率.

Rrrrrrrrrrrrrrrr 啊啊啊啊啊啊啊啊啊啊….. (聲音漸尖)

R 小筆記

Realtek 公司也簡稱 R, 不過應該沒有多少人用過 R 語言. 最近適逢要用到 regression model, 但網路課程幾乎都是用 R (歷史包袱), 而不是 Python, 所以我也試著瞭解用 R 來理解 regression model.

首先 library(UsingR) 載入基本套件, 然後用 x 來 predict y.

y = β0 + β1 * X

在程式中, 以 beta0 取代 β0 , 而 beta1 取代β 1 , 因為我們只是 predict, 並不是 assign 值給 y, 所以要求的是 minimum square error.

i=1n​(Yi−(β0​+β1Xi))2

為了增加趣味性, 課程舉例是用父母的身高預測小孩的身高. 我們可以理解到身高不會為 0, 所以有一個 beta0 撐著, 然後父母身高 X 乘上的 weighting beta1, 就可以預測小孩的身高. 當然, 我們又可以理解到這是在統計上才有意義, 對個案沒有意義.

lm(y~x) 就可以得到用 x 預測 y 的 linear model 的結果.

coef(lm(y~x)) 回報 coefficients beta0 和 beta1.

c(beta0, beta1) 用來產生向量.

舉例來說, x 和 y 用 c() 產生向量, 得出 linear model 後, model 再用 coef() 取出係數, 然後印出來.


x <- c(1, 2, 3, 4, 5)
y <- c(2, 3, 5, 7, 11)

# 建立線性回歸模型
model <- lm(y ~ x)

# 獲取模型係數
coefficients <- coef(model)  # This gives c(beta0, beta1)

intercept <- coefficients[1] # beta0
slope <- coefficients[2] # beta1

# 印模型係數
print(intercept)
print(slope)

# OR 一行印模型係數
print(coefficients)

其中 <- 表示 assign, 正式名稱為賦值運算子 (assignment operator), 和 C 的 pointer 反方向.

用 Python 其實也不會變囉嗦. 以這個例子來說, Python 可以輕鬆應對. 雖然 R 確實更簡潔一點點.

import numpy as np
from sklearn.linear_model import LinearRegression

x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 3, 5, 7, 11])

# 建立線性回歸模型
model = LinearRegression()
model.fit(x.reshape(-1, 1), y)

# 獲取模型係數
intercept = model.intercept_
slope = model.coef_[0]

# 印模型係數
print(f"Intercept (beta0): {intercept}")
print(f"Slope (beta1): {slope}")

# OR 一行印模型係數
print(model)