跳轉到

9.2 係數的估計與標準誤#

本節目標#

  • 說明「估計值」與「標準誤」為什麼必須一起報導
  • 會算殘差變異數 \(s^2\) 與係數的標準誤
  • 逐行讀懂 summary(lm) 的輸出：Estimate、Std. Error、t value、Pr(>|t|)
  • 理解 \(R^2\) 的意義，以及單元迴歸中 \(R^2 = r^2\) 這個好用的關係
  • 知道信賴區間是怎麼從「估計值 ± 標準誤」來的（完整版在 9.5）

問題場景：一個數字不夠#

上一節我們用最小平方法得到 Galton 資料的擬合線：

\(\hat{y} = 23.9415 + 0.6463 x\)

如果只把這兩個數字寫在報告上，讀者會以為「斜率就是 0.6463，這麼精確」。但我們都知道這是從 928 個家庭抽出來的估計值——如果重新抽 928 個家庭，斜率會不會很不一樣？

這就是統計上「估計值（estimate）」與「估計值的精度（standard error）」的分工：

  • 估計值：點估計，回答「大概是多少」。
  • 標準誤：估計值本身的變異程度，回答「這個數字有多可靠」。

沒有標準誤的估計值，就像只報身高不報體重誤差的體檢報告。

歷史：這套工具是怎麼來的#

  • 1908 年，William Sealy Gosset（以筆名「Student」發表，因為他任職於愛爾蘭的都柏林 Guinness 啤酒廠，不允許員工發表論文）提出了 t 分配，用來處理小樣本下平均數與標準差的問題。他面對的正是「釀造原料含水量」這種只有十幾筆資料的困境。
  • 1920 年代，R. A. Fisher 把「最小平方法」和「過度決定（overdetermined system）」的想法結合起來：一組方程、未知數比方程少時沒有解，多於方程時用最小平方找「最接近的解」。Fisher 同時引入了變異數分析（把總變異拆解成來源）與 F 分布。
  • 1930 年代，Fisher 與 Yates、以及後來 Harold Hotelling、Jerzy Neyman 與 Egon Pearson，把這些工具整理成今天的推論統計學架構：估計量、標準誤、檢定、信賴區間。

換句話說，「係數 + 標準誤 + t 檢定」這三件套是同一個時代的產物，它們回答的是同一個問題的不同面向：這個係數可能是多少？它有多不確定？它跟 0 有差別嗎？

理論：從殘差到標準誤#

殘差變異數 \(s^2\)#

擬合之後，每個資料點都留下一個殘差 \(e_i = y_i - \hat{y}_i\) 。如果模型對了，這些殘差應該是「隨機的、不含系統性的」，它們的大小代表資料本身的噪音程度。我們用殘差平方和估計這個噪音：

\(s^2 = \hat{\sigma}^2 = \frac{\sum e_i^2}{n - 2}\)

為什麼除以 \(n-2\) 而不是 \(n\) ？因為自由度：估計 \(b_0\) 與 \(b_1\) 兩個參數之後，資料只剩下 \(n-2\) 個獨立的自由度。分母用 \(n\) 會系統性低估變異數（這與第 6 章「用樣本變異數估計母體變異數」是同一件事）。

\(s = \sqrt{s^2}\) 稱為殘差標準差（residual standard error），它的單位和 \(y\) 一樣，所以最好解釋。

斜率的標準誤#

簡單迴歸中，斜率估計值 \(b_1\) 的標準誤是

\(\mathrm{SE}(b_1) = s \sqrt{\frac{1}{S_{xx}}} = \frac{s}{\sqrt{\sum (x_i - \bar{x})^2}}\)

這個式子告訴我們兩件事：

  1. ** \(x\) 的離散程度越大，斜率估得越準。** 因為 \(S_{xx}\) 在分母， \(x\) 攤得越開，斜率就越容易被「看清楚」。
  2. 樣本數越多， \(S_{xx}\) 越大，斜率越準。 這是「資料越多越好」的第一個量化理由。

截距的標準誤#

\(\mathrm{SE}(b_0) = s \sqrt{\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}}\)

注意截距的標準誤比斜率大得多（當 \(\bar{x} \ne 0\) 時）。原因很直覺：截距是直線在 \(x=0\) 處的高度，必須靠「整條線往回推」來確定，而這正是最不穩定的部分。上一節 Galton 的 \(\mathrm{SE}(b_0) = 2.81\) 吋 vs \(\mathrm{SE}(b_1) = 0.041\) ，就是這個現象。

t 統計量與 p 值#

\(t = \frac{b_1 - 0}{\mathrm{SE}(b_1)}\)

在「誤差為常態、變異數固定」的假設下，這個 t 統計量服從自由度 \(n-2\) 的 t 分配（見 5.4 節）。所以 p 值就是

\(p = P\left(|T_{n-2}| \ge |t|\right)\)

而 95% 信賴區間（概念，完整公式在 9.5）是

\(b_1 \pm t_{0.975, n-2} \cdot \mathrm{SE}(b_1)\)

三件事是同一件事的三種說法：檢定（係數是否為 0）、信賴區間（係數可能落在哪裡）、p 值（在「係數為 0」的假設下，會看到這麼極端的資料的機率）。

動手：讀懂 summary(lm)#

資料用 R 內建的 women（15 位女性的體重與身高），小到能手算。

import rstatistics as r

HEIGHT = [58.0, 59.0, 60.0, 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 67.0,
          68.0, 69.0, 70.0, 71.0, 72.0]
WEIGHT = [115.0, 117.0, 120.0, 123.0, 126.0, 129.0, 132.0, 135.0, 139.0, 142.0,
          146.0, 150.0, 154.0, 159.0, 164.0]

m = r.lm(WEIGHT, [HEIGHT], names=["height"])
sm = m.summary()
print(sm)

# ---------- 手算：殘差變異數、標準誤、t 值、p 值 ----------
n = len(HEIGHT)
res = [yi - fi for yi, fi in zip(WEIGHT, m.fitted_values)]
rss = sum(e * e for e in res)
s2 = rss / (n - 2)
sd_res = s2 ** 0.5
xb = r.mean(HEIGHT)
yb = r.mean(WEIGHT)
sxx = sum((xi - xb) ** 2 for xi in HEIGHT)
se_b = sd_res / sxx ** 0.5
se_a = sd_res * (1.0 / n + xb * xb / sxx) ** 0.5
t_b = m.coefficients[1] / se_b
p_b = 2 * r.pt(abs(t_b), n - 2, lower_tail=False)

print("手算對照：")
print("  n = %d, 自由度 = %d" % (n, n - 2))
print("  RSS = %.4f" % rss)
print("  s² = RSS/(n-2) = %.4f, s = %.4f   (summary: %.4f)" % (s2, sd_res, sm.sigma))
print("  SE(b) = s/sqrt(Sxx) = %.5f   (summary: %.5f)" % (se_b, sm.std_errors[1]))
print("  SE(a) = s*sqrt(1/n+xbar^2/Sxx) = %.5f   (summary: %.5f)" % (se_a, sm.std_errors[0]))
print("  t(b) = %.4f   (summary: %.4f)" % (t_b, sm.t_values[1]))
print("  p(b) = %.4g   (summary: %.4g)" % (p_b, sm.p_values[1]))

# ---------- R² 與 cor 的關係 ----------
rho = r.cor(HEIGHT, WEIGHT)
print()
print("R² = %.6f" % sm.r_squared)
print("cor(x, y) = %.6f, cor² = %.6f" % (rho, rho * rho))
print("R² == cor² ? %s" % (abs(sm.r_squared - rho * rho) < 1e-12))

# ---------- 係數的信賴區間（概念，完整版在 9.5）----------
tcrit = r.qt(0.975, n - 2)
b = m.coefficients[1]
a = m.coefficients[0]
print()
print("95%% CI for 斜率 b: [%.4f, %.4f]" % (b - tcrit * se_b, b + tcrit * se_b))
print("95%% CI for 截距 a: [%.4f, %.4f]" % (a - tcrit * se_a, a + tcrit * se_a))

輸出：

Call:
lm(formula = y ~ height)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
height        3.45000    0.09114   37.85 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared:  0.991, Adjusted R-squared:  0.9903 
F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

手算對照：
  n = 15, 自由度 = 13
  RSS = 30.2333
  s² = RSS/(n-2) = 2.3256, s = 1.5250   (summary: 1.5250)
  SE(b) = s/sqrt(Sxx) = 0.09114   (summary: 0.09114)
  SE(a) = s*sqrt(1/n+xbar^2/Sxx) = 5.93694   (summary: 5.93694)
  t(b) = 37.8553   (summary: 37.8553)
  p(b) = 1.091e-14   (summary: 1.091e-14)

R² = 0.991010
cor(x, y) = 0.995495, cor² = 0.991010
R² == cor² ? True

95% CI for 斜率 b: [3.2531, 3.6469]
95% CI for 截距 a: [-100.3427, -74.6907]

逐行解讀輸出的每一個欄位#

第一段 Residuals: — 殘差的五數摘要。Min -1.7333、Max 3.1167 是最小的與最大的殘差（本例資料幾乎是一條直線，所以殘差很小）；Median -0.3833 表示殘差的中位數稍微偏負，也就是資料略微呈現曲率（下一節的診斷主題）。

Coefficients: 表格 — 整張表的重點：

欄位 意義 本例（height 那一列）
Estimate 係數的點估計 3.45：身高每多 1 公分，體重平均多 3.45 公斤
Std. Error 係數的標準誤 0.09114：斜率的估計誤差約 ±0.09
t value 檢定統計量 \(t = \hat\beta / \mathrm{SE}\) 37.85
Pr(>|t|) 雙尾 p 值 1.09e-14：極不可能是巧合
*** 顯著性標記 p < 0.001

Residual standard error: 1.525 on 13 degrees of freedom — 這就是 \(s = \sqrt{\mathrm{RSS}/(n-2)}\) 。它的意思是「扣掉模型能解釋的部分之後，資料還剩下約 1.5 公斤的典型誤差」。

Multiple R-squared: 0.991, Adjusted R-squared: 0.9903 — 模型解釋了 99.1% 的變異。調整 \(R^2\) ** 則把「用掉一個參數」的成本扣掉：在多元迴歸中多加一個變數幾乎一定會讓 \(R^2\) 變大（哪怕它一點用處都沒有），所以比較不同變數個數**的模型時，要看調整過的版本。

F-statistic: 1433 on 1 and 13 DF — 這是「整個模型是否比『只有截距』的模型好」的檢定。在只有一個預測變數時， \(F = t^2\) （本例 \(37.85^2 \approx 1433\) ），所以 F 檢定和 t 檢定是同一件事的兩種寫法；多元迴歸時 F 檢定才有獨立意義。

\(R^2 = r^2\) 這個好用的關係#

在簡單（單預測變數）迴歸中，有一個可以直接背下來的等式：

\(R^2 = r^2\)

其中 \(r = \mathrm{cor}(x, y)\) 是皮耳森相關係數。程式也驗證了這件事：

R² = 0.991010
cor(x, y) = 0.995495, cor² = 0.991010
R² == cor² ? True

為什麼會這樣？ 因為在單元迴歸中， \(R^2 = 1 - \frac{\mathrm{RSS}}{\mathrm{TSS}}\) ，而 \(\mathrm{RSS} = \mathrm{TSS}(1-r^2)\) ，所以 \(R^2 = r^2\) 。這也解釋了為什麼 \(R^2\) 有個上限 1：相關係數的絕對值不超過 1。

這條等式有個實際用途：當你只想要「線性關係的強度」而不關心截距時，相關係數就是 \(R^2\) 的平方根，而且相關係數可以在任何一對變數之間計算，不一定要做迴歸。

信賴區間：估計值 ± 一個Margin#

看最後兩行：

95% CI for 斜率 b: [3.2531, 3.6469]
95% CI for 截距 a: [-100.3427, -74.6907]

斜率的區間不包含 0，所以「身高對體重有影響」的檢定會被拒絕；截距的區間則離 0 非常遠（因為身高不可能是 0），這個檢定在實務上毫無意義——這正是為什麼截距的 p 值通常沒有人看。

與 R 對照#

# 係數的估計與標準誤：R 對照（資料：women，身高體重）
data(women)

fit <- lm(weight ~ height, data = women)
sm <- summary(fit)
cat("=== summary(lm) ===\n")
print(sm)

cat("\n=== 手算：殘差變異數與標準誤 ===\n")
rss <- sum(residuals(fit)^2)
s2 <- rss / (nrow(women) - 2)
cat(sprintf("RSS = %.4f, s2 = %.4f, s = %.4f\n", rss, s2, sqrt(s2)))
xb <- mean(women$height); sxx <- sum((women$height - xb)^2)
cat(sprintf("SE(b) = %.5f, SE(a) = %.5f\n",
            sqrt(s2) / sqrt(sxx), sqrt(s2) * sqrt(1 / nrow(women) + xb^2 / sxx)))

cat("\n=== 係數的信賴區間 confint() ===\n")
print(confint(fit))

cat("\n=== R squared 與 cor ===\n")
cat(sprintf("R-squared = %.6f, cor = %.6f, cor^2 = %.6f\n",
            sm$r.squared, cor(women$height, women$weight),
            cor(women$height, women$weight)^2))

cat("\n=== 變異數共變異矩陣 vcov() ===\n")
print(vcov(fit))

輸出：

=== summary(lm) ===

Call:
lm(formula = weight ~ height, data = women)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
height        3.45000    0.09114   37.85 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared:  0.991, Adjusted R-squared:  0.9903 
F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14


=== 手算：殘差變異數與標準誤 ===
RSS = 30.2333, s2 = 2.3256, s = 1.5250
SE(b) = 0.09114, SE(a) = 5.93694

=== 係數的信賴區間 confint() ===
                  2.5 %     97.5 %
(Intercept) -100.342655 -74.690679
height         3.253112   3.646888


=== R squared 與 cor ===
R-squared = 0.991010, cor = 0.995495, cor^2 = 0.991010

=== 變異數共變異矩陣 vcov() ===
            (Intercept)       height
(Intercept)   35.247305 -0.539880952
height        -0.539881  0.008305861

vcov() 印出的是係數的變異數共變異矩陣。它的對角線就是標準誤的平方（0.008305861 = 0.09114²），非對角線元素則是斜率與截距的共變異：這個例子裡是負的（-0.54），意思是「 \(x\) 的抽樣誤差讓斜率偏大時，截距會偏小」。多元迴歸中，這個矩陣的結構會讓係數的解讀變複雜（9.5 會再遇到）。

小結#

  • 標準誤是估計值的精度：沒有標準誤的係數是無法評估的。
  • \(s^2 = \mathrm{RSS}/(n-2)\) 、 \(s = \sqrt{s^2}\) 是殘差變異數與殘差標準差；除以 \(n-2\) 是因為估了兩個參數。
  • \(\mathrm{SE}(b_1) = s/\sqrt{S_{xx}}\) ； \(x\) 越分散、樣本越大，斜率估得越準。
  • summary(lm) 的每一欄都有意義：Estimate（點估計）、Std. Error（精度）、t value（係數是幾個標準誤 away from 0）、Pr(>|t|)（巧合的機率）。
  • 單元迴歸中 \(R^2 = r^2\) ；報告模型時要一併給出 \(R^2\) 與殘差標準差 \(s\) 。
  • 信賴區間、檢定、p 值是同一件事的三種說法。

練習#

  1. 用 r.lm 擬合 weight ~ height 但設 intercept=False，比較係數與 summary 的輸出（ \(t value\) 那一列會變成什麼？為什麼？）。提示：intercept=False 時正規方程少了一條， \(t\) 檢定的自由度也跟著變。
  2. 手算 \(\mathrm{SE}(b_1)\) 時，把 \(s\) 換成 \(\sqrt{\mathrm{RSS}/n}\) （少除一個自由度），看看 \(t\) 值變大多少。這解釋了為什麼自由度要用 \(n-2\) 嗎？
  3. 把 women 的資料點減半（例如只留身高 <= 65 的），重做擬合。 \(R^2\) 變大還是變小？ \(\mathrm{SE}(b_1)\) 變大還是變小？兩者為什麼可能往不同方向走？
  4. 用 r.cor 算出 women 的相關係數，再用 \(R^2 = r^2\) 反推 \(R^2\) ，和 summary 印出的 Multiple R-squared 比較（應該差在浮點誤差內）。
  5. 想一想：為什麼截距的 p 值在實務上通常沒有人關心？在什麼情況下截距反而是有意義的（例如沒有原點的物理量、或做「校正」用的模型）？