---
title: 【R 語言與統計資料分析】Ch14：OLS 迴歸的統計診斷與多元線性迴歸模型
image: https://ppt.cc/fnT4xx@.png
---

# 【R 語言與統計資料分析】Ch14：OLS 迴歸的統計診斷與多元線性迴歸模型

在前一個章節，我們成功建構出一個簡單線性迴歸模型，並且也對其係數進行 t 檢定，但這個**檢定是準確的嗎**，會不會有犯下什麼錯誤呢?

另外，簡單線性迴歸只能描述兩個變數間的關係，但真實世界是更加複雜，我們也得增加線性模型的複雜度，也就是所謂**多元迴歸模型**。

# 14.1 迴歸分析的統計診斷

在進行假設檢定時，我們自然是需要了解統計量的統計分佈，而關於線性迴歸分析有一些經典的假設得要先去確認。

讓我們接續上一章節的模型，繼續作為本章的例子:

# 讀取網路上的 CSV 檔
engel <- read.csv("https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/quantreg/engel.csv")

# 移除自動產生的索引欄位 
engel <- engel[, -1]


# 建立 OLS 模型
# foodexp: 食物支出 (Y), income: 收入 (X)
ols_model <- lm(foodexp ~ income, data = engel)

# 查看初步結果
summary(ols_model)

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/Skq5CnudZx.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">簡單線性迴歸報表</span>
</div>

## 14.1.1 簡單線性迴歸的古典統計假設

我們假設模型為：$Y_i = \beta_0 + \beta_1 X_i + e_i, \quad i = 1, \dots, n$
1. **隨機抽樣假設** (Random Sampling)：
$\{(Y_i, X_i)\}_{i=1}^n$ 是來自母體的獨立同分佈（i.i.d.）隨機樣本。這保證了觀測值之間的獨立性。
3. **零條件均值** (Zero Conditional Mean)： $E[e_i | X_i] = 0$
這是迴歸分析最核心的假設。它意味著所有未觀察到的影響因素（誤差項 $e$）與我們觀察到的 $X$ 無關。若此假設成立，則 $\hat{\beta}_1$ 是 $X$ 對 $Y$ 迴歸的不偏估計量，也就表示其不存在內生性的問題。
4. **變異數齊一性** (Homoscedasticity)：$Var[e_i | X_i] = \sigma^2$
表示誤差項的變異數不隨 $X$ 的大小而改變。無論 $X$ 是大是小，$Y$ 繞著迴歸線波動的幅度是一樣的。
5. **無自相關** (No Autocorrelation)：$Cov[e_i, e_j] = 0, \forall i \neq j$
第 $i$ 筆資料的誤差跟第 $j$ 筆資料的誤差沒有關係。在橫斷面資料通常成立，但在時間序列資料中常被違反。

5. 誤差符合**常態分配**：
$$\begin{aligned}
e_i | X_i \sim N(0, \sigma^2) &\Rightarrow (Y_i - \beta_0 - \beta_1 X_i) | X_i \sim N(0, \sigma^2) \\
&\Rightarrow Y_i | X_i \sim N(\beta_0 + \beta_1 X_i, \sigma^2)
\end{aligned}$$

## 14.1.2 同質變異數與異質變異數

<table style="margin-left: auto; margin-right: auto; width: 90%; border-collapse: collapse; border: none; text-align: center;">
<tr style="border: none;">
<td style="border: none; width: 45%; padding: 10px;">
<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/b/bf/Homocedastico.JPG/250px-Homocedastico.JPG" style="width: 100%; max-width: 450px;">
<br>
<span style="color: gray; font-size: 0.9em;">同質變異數 (Homoscedasticity)</span>
</td>
<td style="border: none; width: 45%; padding: 10px;">
<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/6/68/Heterocedastico.JPG/250px-Heterocedastico.JPG" style="width: 100%; max-width: 450px;">
<br>
<span style="color: gray; font-size: 0.9em;">異質變異數 (Heteroscedasticity)</span>
</td>
</tr>
</table>

在線性迴歸的概念中，迴歸直線表示當 $X=x$ 時，$Y$ 的平均取值為何。所謂「**同質變異數**」的意思是：「當 $X$ 的數值在變化時，其預測誤差 (殘差) 服從之常態分配變異數是一致的。」反之，所謂「**異質變異數**」即為殘差分配的變異數會根據 $X$ 的變化而改變。

**異質變異數**所導致的問題就是會讓我們在進行係數估計時使用的變異數會有所偏誤，其中一個解決方式即為使用更為穩健 (robust) 的變異數估計。

## 14.1.3 檢測是否存在異質變異數問題

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/ryK-zpuObe.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">簡單線性迴歸與殘差</span>
</div>

在這張圖中，我們可以看到當收入越高時，其花在食物消費上的變異數越來越大，殘差的波動就有點像是一個喇叭狀的圖隨著收入而擴大。

在 `lm()` 函數中，有內建殘差分析的視覺化圖形，我們再一一與各位做說明。

#### 殘差分析圖形

# 繪製殘差分析圖

plot(ols_model, which = 1) 
# 取出第一張圖，如果沒有設置，會產生連續 4 張圖形

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/HyyCID3OWx.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">殘差分析圖，觀察殘差與配適值間的分佈</span>
</div>

這張圖形表示殘差的分布狀況，若殘差平均分佈在 0 的兩側，並且離散程度沒有隨著配適值 ($\hat{Y}$) 而變化，則可以認為此迴歸模型符合同質變異數的假設。然而，我們發現 `lm(foodexp ~ income)` 的殘差分析呈現漏斗狀，這就有可能是存在異質變異數的現象。

# Scale-Location 圖
plot(ols_model, which = 3) 
# 繪出第三張圖

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/B1v2hDn_bl.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">簡單線性迴歸與殘差</span>
</div>

此圖透過分析 $\sqrt{|\text{Standardized Residuals}|}$ 與配適值間的關係，若模型符合同質變異數，則此處的殘差應跟會呈現水平均勻分佈。

#### 檢定是否存在異質變異數

**Breusch–Pagan test** 是檢測異質變異數的一種方法，其虛無假設與對立假設如下：

- $H_0$：**Homoscedasticity**
- $H_1$：**Heteroscedasticity**

# 若尚未安裝，請先執行 install.packages("lmtest")
library(lmtest)

# 執行 Breusch-Pagan 檢定
bptest(ols_model)

<div style="text-align: left;">
<img src="https://hackmd.io/_uploads/SJzdDdndZg.png" width="250">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

此結果讓我們拒絕同質變異數的虛無假設，因此此模型存在**異質變異數**。

## 14.1.4 檢測是否符合常態假設

#### Q-Q plot

這張圖用來檢查「**殘差常態性假設**」，若模型符合常態假設，線性迴歸假設我們的預測誤差應該呈現鐘型曲線，大多數誤差很小，不會有極端誤差。

- X 軸 (Theoretical Quantiles / 理論分位數)：
電腦模擬出一個完美的「標準常態分佈（鐘型曲線）」，計算出對應的分位數值，把這些理論上的數值放在 X 軸。

- Y 軸 (Sample Quantiles / 樣本分位數 或 Standardized Residuals)：
把我們模型算出來的實際殘差，由小到大排序，以分位數的形式，放在 Y 軸。

若模型符合常態分配，其 X 軸 與 Y 軸分位數間的位置應該是一一對應，完美呈現 45 度。

# qqplot
plot(ols_model, which = 2) 
# 繪出第二張圖

<div style="text-align: left;">
<img src="https://hackmd.io/_uploads/SJBbEt3dbg.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

#### 檢定分配是否為常態

我們可以使用 **Shapiro–Wilk test** 來檢驗分配是否為常態分配。

- $H_0$：**Normality**
- $H_1$：**Not Normality**

# 假設 model 是你已經建立好的線性迴歸模型
# 提取模型的殘差
res <- residuals(ols_model)

# 執行 Shapiro-Wilk 常態性檢定
shapiro.test(res)

<div style="text-align: left;">
<img src="https://hackmd.io/_uploads/ryeOoi2Obe.png" width="250">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

我們會拒絕虛無假設，轉而接受此模型殘差非常態的對立假設。

> 我們仍要注意，若檢定無法拒絕虛無假設 (即 $p>0.05$)，並不代表我們肯定認為殘差就是服從常態分配。

# 14.2 內生性的問題

假設我們現在有一個迴歸模型：

$$Y_i = \beta_0 + \beta_1 X_i + e_i, \quad i = 1, \dots, n$$

若我們想分析 $X$ 與 $Y$ 之間的關係，尤其是「**因果**」關係時，我們就必須面對到所謂「**內生性**」的問題。

所謂內生性的問題即為模型中的 $e$ 與迴歸變數 $X$ 存在相關性，也就是：
$$Cov(X, e) \neq 0$$

內生性問題的來源基本上可以分為以下幾種：

1. **遺漏變數偏誤** (Omitted Variable Bias)
若存在一個遺漏變數 $Z_i$ ，跟原本的 $X_i$ 有所關係，並且會影響到 $Y_i$ ，那麼誤差項 $u_i$ 也自然會跟 $X_i$ 有關係。
- 例：天賦高（如：智商）的人通常容易讀到較高的學歷（影響 $X$），同時天賦高的人即使不讀書，工作表現通常也較好、薪水較高（進入 $\epsilon$ 影響 $Y$），此時 $Cov(X, \epsilon) > 0$。

2. **互為因果關係** (Reverse Causality)
$X$ 會影響到 $Y$ ，但 $Y$ 同時又會反過來影響 $X$，這時 $X$ 和 $Y$ 是在一個系統內互相決定的狀態。
- 例：增加警察能降低犯罪率（$X \to Y$）；但犯罪率越高的危險地區，政府通常會派駐更多警察（$Y \to X$）。如果沒有仔細思考，甚至可能得出「警察越多，犯罪率越高」的荒謬結論。

3. **測量誤差** (Measurement Error)
若我們對 $X_i$ 的資料搜集的過程有所偏誤，也就是搜集到的資料並不是真實的 $X_i$ 這會導致模型誤差項 $u_i$ 和錯誤的 $X_i$ 有相關性。

內生性問題除了影響因果關係之外，也會導致係數估計出現偏誤，這就會使得我們的研究過程出現錯誤。

# 14.3 多元迴歸分析 Multiple Regression

先前，我們曾經討論過簡單線性迴歸的架構，我們希望用**一個**解釋變數 $X$ 去解釋果變數 $Y$。模型設定如下：

- 簡單線性迴歸模型：
$Y_i=\underbrace{\hat{\beta_0}+\hat{\beta_1}}_{estimators}X_i + \underbrace{u_i}_{residual}\ ,i=1,2,\cdots,n$

而現在我們要擴張這個架構，將**多個解釋變數**納入模型，而這也很符合直覺，因為事物的成因往往都不是那麼單一，同時也避免了「**遺漏變數偏誤**」。
除此之外，多元迴歸的架構也容許我們在模型中放入「**控制變數**」，來確保模型正確的**因果關係**。一般而言，**多元迴歸模型**設定如下：

$$Y_i=\hat{\beta_0}+\hat{\beta_1}X_{1i} +\hat{\beta_2}X_{2i}+\cdots+\hat{\beta_k}X_{ki}+ u_i\ ,i=1,2,\cdots,n$$

> 我們可以將簡單線性迴歸當作是多元迴歸的 $k=1$ 時的一種特例。

若我們將多元線性迴歸模型寫成**矩陣**的形式：

多元線性迴歸模型(矩陣版)：
$\textbf{Y} = \textbf{X}\beta + \textbf{e}$ (母體)
$\textbf{Y} = \textbf{X}\hat{\beta} + \hat{e}$ (估計)

其估計係數之解為：$\hat{\beta}=(\textbf{X}^T\textbf{X})^{-1}\textbf{X}^T\textbf{Y}$

:::spoiler 多元線性迴歸模型(矩陣版)
多元線性迴歸模型(矩陣版)：
$\textbf{Y} = \textbf{X}\beta + \textbf{e}$ (母體)
$\textbf{Y} = \textbf{X}\hat{\beta} + \hat{e}$ (估計)
其中，

$\textbf{Y}_{n\times1} =\begin{bmatrix}Y_1 \\ Y_2 \\\vdots\\Y_n\\ \end{bmatrix}$

$\textbf{X}_{n\times (k+1)} =\begin{bmatrix}1 &X_{11}   &X_{21} & \cdots &X_{k1} \\ 1 & X_{12}  &X_{22}  &\cdots & X_{k2}\\\vdots& \vdots  & \vdots& \vdots & \vdots \\1&X_{1n}     &X_{2n} &  \cdots &X_{kn}\\ 
\end{bmatrix}$

$\bf{\beta}_{(k+1)\times1} =\begin{bmatrix}\beta_0 \\ \beta_1\\\beta_2 \\\vdots\\\beta_k\\ \end{bmatrix}$

$\textbf{e}_{n\times1} =\begin{bmatrix}e_1 \\ e_2 \\\vdots\\e_n\\ \end{bmatrix}$

我們所求出之 $\hat{\beta}$ 如下所示：

$\hat{\beta}=(\textbf{X}^T\textbf{X})^{-1}\textbf{X}^T\textbf{Y}$

若要求出這個答案，我們會要求 $(\textbf{X}^T\textbf{X})$ 的**反矩陣存在**，也就是 $(\textbf{X}^T\textbf{X})$ 是**滿秩 full rank**。

:::

## 14.3.1 多元迴歸模型的解釋

一般而言，我們可以設定這個模型如下：

$$Y =f(\textbf{X})= {\beta_0}+{\beta_1}X_{1} + {\beta_2}X_{2}+\cdots+ {\beta_k}X_{k}+ e$$

對特定變項的變動可以這樣子看：

$\frac{\partial f(\textbf{X})}{\partial X_k} = \beta_k$

也就是說，$\beta_k$ 是在「**其他條件不變之下」**$X_k$ 變動對於 $Y$ 的影響程度。

在二維的情況之下，我們希望用一條「**線**」去預測資料，而在三維的情況下，我們則是希望用一個平面去進行預測，而在更高的維度，我們已無法運用圖形去呈現，若想進行資料視覺化，就需要進行高維度往低維度的投影。

<div style="text-align: center;">
<img src="https://media.geeksforgeeks.org/wp-content/uploads/20240627122545/Hyperplane.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">迴歸模型用超平面來預測資料，在二維平面上就是直線，在三維空間中就是平面</span>
</div>

<div style="text-align: center;">
<img src="https://statsandr.com/blog/multiple-linear-regression-made-simple/images/multiple-linear-regression-plane.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">多元迴歸模型中實際資料與預測平面間的關係</span>
</div>

## 14.3.2 多元迴歸模型的統計假設

1. $\{(Y_i,X_i)\}_{i=1}^n$ 是來自 random sample
2. $\textbf{X}_{n \times k}$ 是 full column rank：為了能做最小平方法，也就是**沒有完全多重共線性**的問題
3. $\textbf{E}[\textbf{e}|\textbf{X}]=0$：沒有內生性問題
4. **變異數齊一性** (Homoscedasticity)：
- $Var[\textbf{e}|\textbf{X}]=\sigma^2I \Leftrightarrow Var[e_i|X_{1i},X_{2i},\cdots,X_{ki}]=\sigma^2$
- $Cov[e_i,e_j]=0, \forall i \neq j$

5. 誤差 error 服從**常態分配**：$\textbf{e}|\textbf{X} \sim N(0,\sigma^2I)\Leftrightarrow e_i|\textbf{X} \sim N(0,\sigma^2)$
透過這個假設，舉例來說，我們可以將 $\textbf{e}$ 改寫如下：

$$\begin{aligned}e_i|\textbf{X} \sim N(0,\sigma^2)&\Rightarrow (Y_i-\beta_0-\beta_1 X_1-\beta_2 X_2)|\textbf{X}  \sim N(0,\sigma^2)\\&\Rightarrow Y_i|\textbf{X} \sim N(\beta_0+\beta_1 X_1+\beta_2 X_2,\sigma^2)
\end{aligned}$$

我們就可以看到 $Y$ 服從一個常態分配，而他的中心點 (mean) 就是 $X$ 的線性組合。

## 14.3.3 完全多重共線性 Perfect Multicollinearity

我們一樣看到多元回歸模型：$Y =f(\textbf{X})= {\beta_0}+{\beta_1}X_{1} + {\beta_2}X_{2}+\cdots+ {\beta_k}X_{k}+ e$

在進行多元迴歸分析時，我們最希望的情況是 $X_i$ 與 $X_j$ 是沒有關係的，代表我們將形成 $Y$ 的所有變因都被分門別類的衡量出來。

在數學上，就是前述假設中的 $\textbf{X}_{n \times k}$ 是 **full column rank**，這是因為在估計 $\beta$ 時，最小平方法的公式為：$\hat{\beta} = (X^T X)^{-1} X^T Y$，若 $\textbf{X}_{n \times k}$ 不為 full column rank，$(X^T X)$ 就不會是 **full rank**，會導致**反矩陣無法計算，線性迴歸模型會直接無法使用**，而這正是所謂**完全多重共線性**。

:::success
💡 **資料的滿秩 Full Rank**
假設我們把 5 個自變數（$X_1$ 到 $X_5$）丟進模型，而這 5 個變數各自都帶有別人無法取代的獨立資訊，那麼這個矩陣的 Rank 就是 5。當「有效變數數量」等於「我們放入的變數總數」時，這就叫做 Full Rank。
:::

而在實務上，$X_i$ 與 $X_j$ 之間的相關係數幾乎不可能剛好等於 0，因此，我們是可以試圖控制 $X_i$ 與 $X_j$ 「**不要過度相關**」即可。

倘若 $X_i$ 與 $X_j$ 的相關性過高，雖然不會導致反矩陣無法運算，但是會導致估計的標準誤不準確，會影響模型的顯著性，此即「**不完全多重共線性**」(imperfect multicollinearity)。

:::success
💡 **標準差 Standard Deviation vs. 標準誤 Standard Error**
- **標準差** (Standard Deviation, SD)： 衡量的是「**實際資料**」之間的變異程度。
- **標準誤** (Standard Error, SE)： 衡量的是「**統計估計值**」的精準度/誤差範圍，比如說對**平均數**、$\beta_1$ 的估計。
:::

# 14.4 多元迴歸模型實戰

我們現在對於多元迴歸模型已有基本認識，讓我們透過實際資料來進行演練，並且隨著資料一起展開說明其中的細節。

我們在此使用 R 內建的資料集 `mtcars`，此為 1974 年美國的 Motor Trend 雜誌對 32 種車款在 11 項設計與效能表現所做的資料整理。

| 欄位名稱 (`Column`) | 資料型態 | 中文解釋義 | 備註 |
| :--- | :--- | :--- | :--- |
| **`mpg`** | 數值 (num) | **油耗 (Miles/(US) gallon)** | 每加侖能跑幾英里。**數值越大越省油**。這是我們最常拿來當作預測目標 ($Y$) 的變數。 |
| **`cyl`** | 數值 (num) | **汽缸數量 (Cylinders)** | 包含 4、6、8 缸。通常汽缸越多，動力越強但也越耗油。常與排氣量產生共線性。 |
| **`disp`** | 數值 (num) | **排氣量 (Displacement)** | 單位為立方英吋 (cu.in.)。引擎大小的指標，數值越大通常越耗油。 |
| **`hp`** | 數值 (num) | **馬力 (Gross horsepower)** | 引擎輸出的功率。馬力越大，車子跑得越快，但也越耗油。 |
| **`drat`** | 數值 (num) | **後軸齒輪比 (Rear axle ratio)** | 數值較高通常代表起步加速較快，但極速較低且較耗油。 |
| **`wt`** | 數值 (num) | **車重 (Weight)** | 單位為「千磅 (1000 lbs)」。例如 3.325 代表 3325 磅。**影響油耗最致命的因素之一**。 |
| **`qsec`** | 數值 (num) | **零四加速 (1/4 mile time)** | 靜止直線加速跑完 1/4 英里所需的「秒數」。**數值越小，代表加速越快、性能越好**（注意方向性）。 |
| **`vs`** | 數值 (num) | **引擎排列 (V/S)** | **0 = V型引擎**，**1 = 直列引擎 (Straight)**。實務上跑模型前，建議轉換為 Factor。 |
| **`am`** | 數值 (num) | **變速箱 (Transmission)** | **0 = 自排 (Automatic)**，**1 = 手排 (Manual)**。 |
| **`gear`** | 數值 (num) | **前進檔位數 (Forward gears)** | 包含 3、4、5 檔。我們可以把它當作類別變數看待 |
| **`carb`** | 數值 (num) | **化油器數量 (Carburetors)** | 早期控制引擎進氣與燃油混合的零件，包含 1 到 8 個。 |

## 14.4.1 資料整理

在進行任何分析前，我們都必須要對資料有清楚的認識。

# 載入資料
data(mtcars)
head(mtcars)

Output:
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

#### 敘述統計

# 敘述統計
library(psych)
describe(mtcars) 

     vars  n   mean     sd median trimmed    mad   min    max  range  skew kurtosis    se
mpg     1 32  20.09   6.03  19.20   19.70   5.41 10.40  33.90  23.50  0.61    -0.37  1.07
cyl     2 32   6.19   1.79   6.00    6.23   2.97  4.00   8.00   4.00 -0.17    -1.76  0.32
disp    3 32 230.72 123.94 196.30  222.52 140.48 71.10 472.00 400.90  0.38    -1.21 21.91
hp      4 32 146.69  68.56 123.00  141.19  77.10 52.00 335.00 283.00  0.73    -0.14 12.12
drat    5 32   3.60   0.53   3.70    3.58   0.70  2.76   4.93   2.17  0.27    -0.71  0.09
wt      6 32   3.22   0.98   3.33    3.15   0.77  1.51   5.42   3.91  0.42    -0.02  0.17
qsec    7 32  17.85   1.79  17.71   17.83   1.42 14.50  22.90   8.40  0.37     0.34  0.32
vs      8 32   0.44   0.50   0.00    0.42   0.00  0.00   1.00   1.00  0.24    -2.00  0.09
am      9 32   0.41   0.50   0.00    0.38   0.00  0.00   1.00   1.00  0.36    -1.92  0.09
gear   10 32   3.69   0.74   4.00    3.62   1.48  3.00   5.00   2.00  0.53    -1.07  0.13
carb   11 32   2.81   1.62   2.00    2.65   1.48  1.00   8.00   7.00  1.05     1.26  0.29

接著，我們可以進行簡單的探索性資料分析。

# 散佈圖與 kde
library(GGally) # 記得執行
library(ggplot2) # 記得執行
ggpairs(mtcars, 
        aes(alpha = 0.5)
        ) +
  theme_bw() 


#### 相關係數熱力圖
我們先前已討論到共線性問題，那當然要來看一下變數之間的相關性啦！此時，相關係數熱力圖就是我們的好夥伴

# 1. 確保已安裝並載入套件
# if (!require("ggcorrplot")) install.packages("ggcorrplot")
library(ggplot2)
library(ggcorrplot)

# 2. 準備資料：篩選您指定的 mtcars 欄位
# mpg, cyl, disp, hp, drat, wt, qsec, carb
selected_mtcars <- mtcars[, c("mpg", "cyl", "disp", "hp", "drat", "wt", "qsec", "carb")]

# 3. 計算相關係數矩陣
corr_matrix <- cor(selected_mtcars)

# 4. 繪製熱力圖 (套用您的設定)
ggcorrplot(corr_matrix, 
           hc.order = TRUE,           # 使用分群演算法排序
           type = "lower",            # 只顯示下三角
           show.diag = TRUE,          # 顯示對角線
           lab = TRUE,                # 顯示相關係數數值
           lab_size = 4,              # 數值文字大小
           colors = c("#6D9EC1", "white", "#E46726"), # 低、中、高顏色
           title = "mtcars 汽車性能參數相關性熱力圖",
           ggtheme = theme_minimal()) +
           theme(text = element_text(family = "PingFang TC"), # 確保 Mac 中文字體正常
           plot.title = element_text(hjust = 0.5, face = "bold", size = 16)) # 標題置中加粗

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/H1btbjAdZx.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">資料間有些處於高度正負相關的類別</span>
</div>

比如我們可以看到汽缸數量 `cyl` 就跟馬力 `hp` 就高度相關，或是車重 `wt` 與排氣量 `disp` 也高度相關。

## 14.4.2 在 R 中執行多元迴歸模型

如同簡單線性迴歸模型一樣，在 R 中一樣是使用 `lm()` 函數來執行多元迴歸模型，我們**以預測油耗** `mpg` **為目標**建立迴歸模型。

讓我們先做一些基本的資料處理：

library(dplyr)
# 先把 gear 轉為正確的類別變數
mtcars_new <- mtcars %>%
                mutate(gear = as.factor(mtcars$gear))

glimpse(mtcars_new)

Output:
Rows: 32
Columns: 11
$ mpg  <dbl> 21.0, 21.0, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2, 17.8, 16.4,…
$ cyl  <dbl> 6, 6, 4, 6, 8, 6, 8, 4, 4, 6, 6, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 8, 8, 8,…
$ disp <dbl> 160.0, 160.0, 108.0, 258.0, 360.0, 225.0, 360.0, 146.7, 140.8, 167.6, 1…
$ hp   <dbl> 110, 110, 93, 110, 175, 105, 245, 62, 95, 123, 123, 180, 180, 180, 205,…
$ drat <dbl> 3.90, 3.90, 3.85, 3.08, 3.15, 2.76, 3.21, 3.69, 3.92, 3.92, 3.92, 3.07,…
$ wt   <dbl> 2.620, 2.875, 2.320, 3.215, 3.440, 3.460, 3.570, 3.190, 3.150, 3.440, 3…
$ qsec <dbl> 16.46, 17.02, 18.61, 19.44, 17.02, 20.22, 15.84, 20.00, 22.90, 18.30, 1…
$ vs   <dbl> 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0,…
$ am   <dbl> 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0,…
$ gear <fct> 4, 4, 4, 3, 3, 3, 3, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 4, 4, 4, 3, 3, 3, 3,…
$ carb <dbl> 4, 4, 1, 1, 2, 1, 4, 2, 2, 4, 4, 3, 3, 3, 4, 4, 4, 1, 2, 1, 1, 2, 2, 4,…

# 觀察資料
head(mtcars_new)

Output:
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

接下來，就讓我們正式進入配適迴歸模型的環節：

# 將 dataframe 內容全部匯入
fit_all <- lm(mpg ~ .,data = mtcars_new)

# mpg ~ .
# 這邊放一個點，是指將全部資料匯入當作Ｘ

# 閱讀報表
summary(fit_all)

<div style="text-align: left;">
<img src="https://hackmd.io/_uploads/HkzYd7JFZx.png" width="380">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

我們可以看到除了 `mpg` 之外的變數都被放入自變數 $\mathbf{X}$ 的位置了。
另外，我們可以看到 `gear` 的欄位因為為類別變數，因此新增為**虛擬變數**，而虛擬變數的處理，我們待未來章節做詳細說明。

若我們要指定特定欄位作為模型的自變數，比如說不要放入 `gear`，可以用以下寫法：

fit_xgear <- lm(mpg ~ . -gear,data = mtcars_new) # 不放 gear

summary(fit_xgear)

<div style="text-align: left;">
<img src="https://hackmd.io/_uploads/r1_j5XkF-g.png" width="380">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

我們除了用 `.` 去操控自變數之外，也可以透過**指定**特定欄位作為自變數：

# 指定欄位，不放 gear
fit <- lm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + carb, 
                data = mtcars)

# 查看詳細統計報表
summary(fit)

<div style="text-align: left;">
<img src="https://hackmd.io/_uploads/rykaNV1F-l.png" width="380">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

我們可以看到目前這個模型的模型的統計效果不太好，係數都不太顯著，那怎麼辦呢？

## 14.4.3 多元迴歸模型的檢定

#### t test
同樣地，我們現在來看到這個模型：

$\begin{aligned}Y_i &= {\beta_0}+{\beta_1}X_{1i} + {\beta_2}X_{2i}+\cdots+ {\beta_k}X_{ki}+ e_i\ ,i=1,2,\cdots+n\end{aligned}$

而我們可以對他做 OLS估計：

$\begin{aligned}Y_i &= \hat{\beta_0}+\hat{\beta_1}X_{1i} + \hat{\beta_2}X_{2i}+\cdots+ \hat{\beta_k}X_{ki}+ \hat{e_i}\ ,i=1,2,\cdots+n\end{aligned}$

而 $\beta_k$ 指的是 $X_k$對 $Y$ 的效果，若無效果的話，我們可以猜想在母體模型之中 $\beta_k=0$。

因此，我們會希望**檢驗 $\beta_k$ 是否顯著異於 0**。

我們可以建構 t-統計量來**檢定回歸係數是否與某數有顯著的差異**。

- **虛無假設與對立假設**
$H_0: {\beta_j} = \beta_{j,0}$
$H_1: {\beta_j} \neq \beta_{j,0}$
- **建構 t-statistic**
$\begin{aligned}t = \frac{\hat{\beta_j}-\beta_{j,0}}{SE(\hat{\beta_j})} \sim t(n-k-1)\end{aligned}$

通常，我們會設定 $\beta_{j,0}=0$，也就是說：

$\begin{aligned}t = \frac{\hat{\beta_j}-0}{SE(\hat{\beta_j})} \sim t(n-k-1)\end{aligned}$

#### F 檢定

若我們想要一次檢定很多個變數是否異於零，則必須使用 **F-統計量**(想法與 ANOVA 相同)。

- **虛無假設與對立假設**
$H_0: {\beta_1} = \beta_2=\beta_3=\cdots=\beta_k=0$
$H_1: 至少其中一個 \beta_j 不為0$
- **建構 F-statistic
$\begin{aligned}F=\frac{ESS/k}{SSR/(n-k-1)}\end{aligned}$**

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/SyTbR4yt-e.png" width="400">
<br>
<span style="color: gray; font-size: 0.8em;">在 summary() 報表的最下方可以看到 F 檢定</span>
</div>

# 14.5 模型的評估與選擇

我們在實務上，常常會遇到變數欄位有多個可能性的情況，也因此可以建立起多種不同的迴歸模型。那自然，如何選擇所謂「**最佳**」模型又是一門學問了。

## 14.5.1 判定係數 $R^2$

透過建立一個迴歸模型，只要給我 $X$，我們就可以對未知的 $Y$ 做估計，如下：

$\hat{Y} = \hat{\beta_0}+\hat{\beta_1}X_1+\hat{\beta_2}X_2+ \cdots +\hat{\beta_k} X_k$

另外，我們再定義一些名詞

- **TSS** (Total Sum of Squares)： $\sum^n_{i=1}( {Y_i}-\bar{Y})^2$
- **ESS** (Explained Sum of Squares)：$\sum^n_{i=1}( \hat{Y_i}-\bar{Y})^2$
- **SSR** (Sum of Squared Residual)：$\sum^n_{i=1}({Y_i}-\hat{Y_i})^2 = \sum^n_{i=1} u_i^2$

根據[**證明**](https://en.wikipedia.org/wiki/Explained_sum_of_squares)，我們可以得知 $TSS =ESS + SSR$。
也就是說，**「資料的總變異」等於「模型能解釋的變異」加上「模型沒解釋到的變異」**。

因此，我們定義判定係數 $R^2$ 來衡量模型的解釋力：

$$\begin{aligned}
R^2 = \frac{ESS}{TSS}=1-\frac{SSR}{TSS}
\end{aligned}$$

直觀來說，**當 $R^2$ 之值越大，表示模型解釋 $Y$ 的能力越好**。

然而，$R^2$ 指標存在一個很顯然的問題，讓我們考慮以下情境：

① $SSR=\sum_{i=1}^{n}(Y_i-\hat{\beta_0}-\hat{\beta_1}X_{1i})^2$

② $SSR'=\sum_{i=1}^{n}(Y_i-\tilde{\beta_0}-\tilde{\beta_1}X_{1i}-\tilde{\beta_2}X_{2i})^2$

由於可選擇的變數變多，$SSR' \leq SSR$ 必定成立，又在兩個情境中的 $TSS$ 相同，故可知：**「當變數增加時，$R^2$ 只會上升不會下降」**。

**但**只看這個**指標真的有意義嗎？**

#### 修正後的判定係數 $\text{Adjusted-}R^2$

當增加自變數時，$R^2$ 會隨之上升，故統計學家們透過加入**懲罰項**來修正 $R^2$，此即 $\text{Adjusted-}R^2$ 或稱 $\bar{R}^2$。

$\begin{aligned}\bar{R^2} = 1-\frac{\frac{SSR}{n-k-1}}{\frac{TSS}{n-1}}=1-\frac{n-1}{n-k-1}\times\frac{SSR}{TSS}\end{aligned}$ 

當變數數量 $k$ **增加**時： $SSR\downarrow \quad \frac{n-1}{n-k-1}\uparrow$。

兩者有**抵銷**的效果，而 $\bar{R^2}$ 也可能因此變成負數。$\bar{R^2}$ 會相較 ${R^2}$ 更能告訴我們，新加入的變數是否真的提高了模型的解釋力。

## 14.5.2 變異數膨脹因子 Variance Inflation Factor

我們知道當模型存在共線性問題時，會導致模型的變異增大，也會導致係數在檢定上不顯著的狀況，我們在此介紹**變異數膨脹因子 (VIF)** 作為衡量共線性問題的一個方法。

**VIF** 是針對每一個自變數分別計算的。對於自變數 $X_j$，其 **VIF** 計算方式如下：
$$VIF_j = \frac{1}{1 - R_j^2}$$

$R_j^2$ 並不是總模型的 $R^2$！而是將 $X_j$ 作為應變數，對其他所有自變數進行迴歸所得到的判定係數。
我們可以這樣想，如果其他變數可以把 $X_j$ 解釋得很完美（$R_j^2$ 很高），分母就會變得很小，**VIF** 就會爆炸。

[一般來說](https://en.wikipedia.org/wiki/Variance_inflation_factor)，我們會認為 **VIF 大於 10** (也有書本是設在 5)，即代表此變數存在共線性問題，而我們可以優先針對這些變數**進行剔除**。

> $VIF_j > 10 \Leftrightarrow R_j^2>0.9$

此外，對於這些要剔除的變數，我們可以往下去探討彼此之間的關聯性，比如說：汽缸數 `cyl` 與車重 `wt` 的 VIF 值都很高是因為汽缸數若增加確實會增加車重，但車重的原因並不盡然是來自於汽缸數，故我們也可以再研究這兩個變數油耗 `mpg` 之間的關聯性。

讓我們以剛剛的模型 `fit` 來作為範例：

# 指定欄位，不放 gear
fit <- lm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb, 
                data = mtcars)

library(car)  # 若沒下載過，請執行 install.packages("car")
vif_values <- sort(vif(fit), decreasing = TRUE) # VIF 由大排到小
print(vif_values)

# 視覺化 VIF 
barplot(vif_values, main = "各變數的 VIF 值", 
        col = "steelblue", family = 'PingFang TC')
abline(h = 10, col = "red", lty = 2) # 加入 10 的標準線

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/r1zgtL1Y-g.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">優先觀察 VIF > 10 的變數</span>
</div>

我們發現**排氣量** `disp` 的 VIF 值是最高的，讓我們優先將他進行移除，此時重新跑一個迴歸模型：

fit1 <- lm(mpg ~ cyl  + hp + drat + wt + qsec + vs + am + gear + carb, 
           data = mtcars)

vif_values <- sort(vif(fit1), decreasing = TRUE) # VIF 由大排到小
print(vif_values)

# 視覺化 VIF 
barplot(vif_values, main = "各變數的 VIF 值", 
        col = "steelblue", family = 'PingFang TC')
abline(h = 10, col = "red", lty = 2) # 加入 10 的標準線

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/Bkgar5UJY-e.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">當 VIF 高的變數配移除後，會連帶影響其他變數的 VIF</span>
</div>

我們發現，原本 `wt` 的 VIF 值也非常的高，但當我們移除 `disp` 後，其 VIF 就獲得大幅的下降。那我們可以繼續移除汽缸數 `cyl`，這代表我們選擇使用車重 `wt` 作為主要的代表變數。

# 再進行一次 VIF
fit2 <- lm(mpg ~   hp + drat + wt + qsec + vs + am + gear + carb, 
           data = mtcars)

vif_values <- sort(vif(fit2), decreasing = TRUE) # VIF 由大排到小
print(vif_values)

# 視覺化 VIF 
barplot(vif_values, main = "各變數的 VIF 值", 
        col = "steelblue", family = 'PingFang TC', ylim = c(0,11))
abline(h = 10, col = "red", lty = 2) # 加入 10 的標準線

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/BklRj8JFWg.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">當 VIF 高的變數配移除後，會連帶影響其他變數的 VIF</span>
</div>

# 匯出模型報表
summary(fit2)

<div style="text-align: left;">
<img src="https://hackmd.io/_uploads/S1i83LJtZl.png" width="380">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

## 14.5.3 向前向後與逐步迴歸

在此小節，我們以 p 值作為衡量變數是否適合被放入模型中的依據。

#### 向後選擇法 Backward selection method

**向後選擇法**是一種「由繁入簡」的變數選擇策略，其概念十分簡單，我們先把所有可能的變數都丟進模型，然後根據設定的指標 (如: p 值) 一個一個進行淘汰，直到萃取出最終模型。

讓我們用 R 語言撰寫一個向後選擇法的小迴圈，在此程式中，我們先把所有可用的變數放入模型，接著將模型內有最高 p 值的變數移除，以此類推，直到模型中 p 值都小於 0.05 才停止:

# 向後選擇法 Backward selection method
# alpha = 0.05

data(mtcars)
model <- lm(mpg ~ ., data = mtcars) #全部變數放入

while(TRUE) {
  # 取得當前模型的係數表
  stats <- summary(model)$coefficients
  
  # 排除截距項 (Intercept)，找出 p-value 最大的變數
  p_values <- stats[-1, 4] 
  max_p <- max(p_values)
  max_var <- names(which.max(p_values))
  
  # 判斷是否大於 0.05
  if(max_p > 0.05) {
    cat("移除變數:", max_var, "(p-value =", round(max_p, 4), ")\n")
    # 更新模型：移除該變數
    model <- update(model, as.formula(paste(". ~ . -", max_var)))
  } else {
    break
  }
}

# 顯示最終模型
summary(model)

#### 向前選擇法 Forward selection method

**向前選擇法**與我們剛剛提到的「向後」演算法不同的點是在於其為「由簡入繁」，我們先從一個**最顯著**的變數開始，再往下選擇可以「提升」此模型的第二個變數，以此類推，直到新增變數無法滿足我設定的優化條件時，我就停止。

以下程式為向前選擇法的實踐，

target <- "mpg"
all_vars <- setdiff(names(mtcars), target)
alpha <- 0.05

# 從空模型開始
model <- lm(mpg ~ 1, data = mtcars)

while(TRUE) {
  # 找出目前不在模型中的變數
  current_vars <- all.names(formula(model))[-1]
  remaining_vars <- setdiff(all_vars, current_vars)
   
  if (length(remaining_vars) == 0) break   # 如果所有變數都丟進去了，就停止
  
  # 暫存區：測試每個新變數加入後的狀況
  candidate_p <- c()
  
  for (v in remaining_vars) {
    test_formula <- update(formula(model), as.formula(paste(". ~ . +", v)))  # 加入新變數UPDATE
    test_model <- lm(test_formula, data = mtcars)
    test_summary <- summary(test_model)$coefficients
    
    # 條件 A：新變數自己必須顯著
    new_var_p <- test_summary[v, 4]
    
    # 條件 B：原本已經在裡面的變數 (除了截距) 也要維持顯著
    # 我們檢查模型中所有變數的 p-value 是否都 < alpha
    all_p_values <- test_summary[-1, 4] 
    is_everyone_sig <- all(all_p_values < alpha)
    
    if (is_everyone_sig) {
      candidate_p[v] <- new_var_p
    }
  }
  
  # 檢查是否有任何候選變數符合「兩項條件」
  if (length(candidate_p) > 0) {
    best_var <- names(which.min(candidate_p))
    min_p <- min(candidate_p)
    
    cat("加入變數:", best_var, "(p-value =", round(min_p, 4), ")\n")
    model <- update(model, as.formula(paste(". ~ . +", best_var)))
  } else {
    cat("停止迭代：剩餘變數若加入會導致既有變數不顯著，或新變數本身不顯著。\n")
    break
  }
}

# 最終結果
summary(model)

#### 逐步選擇法 Stepwise selection method

**逐步選擇法**（Stepwise Selection） 結合了向前選取與向後剔除的優點，克服了單向搜尋「不可逆」的缺陷。它在每加入一個新變數後，都會重新檢視模型內所有既有變數的顯著性。若新變數的加入導致舊變數的貢獻降至門檻以下(如: p 值 < 0.05)，則會果斷將其剔除，直到新加入的變數全都不顯著，無法改變模型就停止。

data(mtcars)
target <- "mpg"
all_vars <- setdiff(names(mtcars), target)
alpha_in <- 0.05   # 進入門檻：低於此值才准進來
alpha_out <- 0.05  # 剔除門檻：高於此值就得離開
max_iter <- 1000   # 設定總步數上限，防止無限迴圈
# 1. 從空模型開始 (只有截距項)
model <- lm(as.formula(paste(target, "~ 1")), data = mtcars)

cat("--- 啟動雙向逐步選擇法 (Stepwise Selection) ---\n")
iter = 1
while(iter < max_iter) {
  iter <- iter + 1
  changed <- FALSE
  
  # --- 【第一階段：向前選取】 ---
  current_vars <- all.names(formula(model))[-1]
  remaining_vars <- setdiff(all_vars, current_vars)
  
  if(length(remaining_vars) > 0) {
    candidate_p <- c()
    for (v in remaining_vars) {
      test_mod <- update(model, as.formula(paste(". ~ . +", v)))
      # 取得新加入變數的 p-value
      p_val <- summary(test_mod)$coefficients[v, 4]
      if (p_val < alpha_in) candidate_p[v] <- p_val
    }
    
    if (length(candidate_p) > 0) {
      best_var <- names(which.min(candidate_p))
      model <- update(model, as.formula(paste(". ~ . +", best_var)))
      cat("[+] 加入新變數:", best_var, "(p =", round(min(candidate_p), 4), ")\n")
      changed <- TRUE
      
      # --- 【第二階段：向後剔除 (核心步驟)】 ---
      # 每次加入新成員後，立即檢查所有舊成員是否依然顯著
      while(TRUE) {
        stats <- summary(model)$coefficients
        if (nrow(stats) <= 1) break # 只有截距項就不檢查
        
        # 排除截距項，找出目前模型中 p-value 最大的變數
        p_values_in_model <- stats[-1, 4, drop = FALSE]
        max_p <- max(p_values_in_model)
        max_var <- rownames(p_values_in_model)[which.max(p_values_in_model)]
        
        if (max_p > alpha_out) {
          cat("[-] 剔除變數:", max_var, "(p =", round(max_p, 4), "，受新變數影響變得不顯著)\n")
          model <- update(model, as.formula(paste(". ~ . -", max_var)))
          changed <- TRUE
        } else {
          break # 所有變數都顯著，停止剔除
        }
      }
    }
  }
  
  # 如果這輪沒有任何「加入」或「剔除」動作，代表模型達到最優狀態
  if (!changed) break
}

cat("--- 逐步選擇結束 ---\n")
summary(model)

## 14.5.4 使用資訊量準則評估模型

我們經過先前的學習已知道，若我們一直將入新的變數，就能提升模型的配適能力 (如: $R^2$)，但反而會加大模型的變異，此即**偏差與方差之權衡 (Bias-Variance Tradeoff)**。

**那麼，如何找到一個平衡的模型大小呢?**
其中一個方法即為使用訊息量準則 (Information criterion)，訊息量準則基本上是在關注模型適配度的同時，對於模型複雜度進行懲罰，我們希望可以在兩者中找到平衡，也就是找到使得訊息量準則之值最小的模型，在統計學領域中最常見的訊息量準則有 **AIC (赤池訊息量準則)** 和 **BIC (貝氏訊息量準則)** 兩者。

#### Akaike Information Criterion 赤池訊息量準則 (AIC)

$$\text{AIC}(k) = \ln(\frac{SSR(k)}{n})+ k \cdot \frac{2}{n}$$

#### Bayesian Information Criterion 貝氏訊息量準則 (BIC)

$$\text{BIC}(k) = \ln(\frac{SSR(k)}{n})+ k \cdot \frac{\ln n}{n}$$

**BIC** 與 **AIC** 不同的地方在於懲罰項的部分與**樣本數**有關，若 $n>8$，則 $\ln n>2$ ，因此 **BIC 對於模型複雜度的懲罰效果常會比 AIC 來得大**，因而會傾向挑選到較為精簡的模型。

在 R 中，我們可以使用 `step()` 函數針對 AIC、BIC 做設定後進行向前、向後與逐步選擇的模型設定。

# 假設全模型與空模型已建立
full_model <- lm(mpg ~ ., data = mtcars)
base_model <- lm(mpg ~ 1, data = mtcars)

# 使用 AIC 進行逐步選擇 (預設 k = 2)
model_aic <- step(base_model, 
                  scope = list(lower = base_model, upper = full_model), 
                  direction = "both",   # c("both", "backward", "forward")
                  trace = 1) # trace=0 關閉中間落落長的運算過程

# 使用 BIC 進行逐步選擇 (將 k 設為 log(n))
n <- nrow(mtcars)
model_bic <- step(base_model, 
                  scope = list(lower = base_model, upper = full_model), 
                  direction = "both", 
                  k = log(n), 
                  trace = 1)

# 比較兩者的變數數量
cat("AIC 選出的變數:", names(coef(model_aic))[-1], "\n")
cat("BIC 選出的變數:", names(coef(model_bic))[-1], "\n")

# 查看報表
summary(model_aic)
summary(model_bic)

此外，亦有其他訊息量準則或是類似的方法，如 HQC (Hannan-Quinn criterion) 或是 $\text{Mallow's C}_p$，可見 [Guerrier et al. (2019)](https://smac-group.github.io/ts/the-family-of-autoregressive-moving-average-models.html#model-selection)。

# 14.6 變數的處理與極端值刪除

## 14.6.1 Box-Cox 轉換

在許多時候，資料並不符合常態假設，但是適度地對資料動刀，就有機會將資料轉為常態分配，如此就可以符合迴歸分析的古典假設。(其中一個最常見的案例就是對數常態分配)
我們在此介紹一個相當常用的資料轉換方法，即為 **Box-Cox transformation**，其公式如下:

$$y(\lambda) = \begin{cases} \frac{y^\lambda - 1}{\lambda} & \text{if } \lambda \neq 0 \\ \ln(y) & \text{if } \lambda = 0 \end{cases}$$

以下是常見 $\lambda$ 值對應的傳統轉換意義:
- $\lambda = 0$ 時，$y(\lambda)$ 為對數變換
- $\lambda = 0.5$ 時，$y(\lambda)$ 為根號變換
- $\lambda = -1$ 時，$y(\lambda)$ 為倒數變換

而接下來的問題就是，我們該如何選擇 $\lambda$

讓我們在 R 中可以很簡單的透過套件計算 `log-likelihood` 值，進而找出最適合的 $\lambda$。

library(MASS)

# 1. 呼叫 Box-Cox 尋找最佳 Lambda，並繪出 log-likelihood 圖
bc_engel <- boxcox(ols_model, lambda = seq(-1, 1, by = 0.05), 
                   plotit = TRUE)
title(main = "Box-Cox for Engel")

# 2. 抓出最佳的 Lambda 值
best_lambda <- bc_engel$x[which.max(bc_engel$y)]
cat("建議的 Lambda 值約為:", round(best_lambda, 3), "\n")

建議的 Lambda 值約為: 0.636 

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/HJwQf6gKbl.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">找出使 log-likelihood 最大的 λ</span>
</div>

我們接下來可以進行資料轉換，由於 0.636 接近 0.5，我們可以嘗試將變數開根號，以獲得更好的解釋性。

# 使用最佳 lambda
engel$foodexp_bc <- (engel$foodexp^best_lambda - 1) / best_lambda
bc_model <- lm(foodexp_bc ~ income, data = engel)

# 開根號
engel$foodexp_sqrt <- sqrt(engel$foodexp)
sqrt_model <- lm(foodexp_sqrt ~ income, data = engel)

# 查看報表
summary(bc_model)
summary(sqrt_model)


# 畫出 qqplot
par(mfrow = c(1, 3))

# 原模型
plot(ols_model, which = 2,
     main = "Original dataset")

# 使用 Box-Cox 轉換最佳 lambda    
plot(bc_model, which = 2,
     main = "Box-Cox transformation, λ =0.636")

# 開根號
plot(sqrt_model, which = 2,
     main = "Squared root transformation")


我們可以看到資料在進行轉換後越來越貼近常態分配，但還是有遇到一些離群值的問題，我們可以後續再進行處理。

#### 將資料取自然對數

在財務領域或是總體經濟領域，我們常常會要去計算成長率、報酬率等指標，而這些指標也往往會是我們分析的重要變數。在迴歸分析的架構之中當我們想要計算變數「變化率」對 $Y$ 的影響時，我們會對其取自然對數 $ln$。舉例來說，若股價是我們的X，對股價資料取 $ln$，迴歸係數 $\beta$ 即表示股價變動 1% 所帶來的平均效果。我們將這類模型分成三類(都以 population model 呈現)：

**自然對數的近似**：

$$\ln(1+x) \approx x, \ as\ x\rightarrow 0$$

1. **Linear-log** model
$\begin{aligned}Y = \beta_0 + \beta_1 \ln (X)+e\end{aligned}$

> $\beta_1$ **如何解釋？**
>

考慮 $X$ 微小變動對 $Y$ 的影響。

$\begin{aligned}Y { + \Delta y}  = \beta_0 + \beta_1 \ln (X{+ \Delta x})+e\end{aligned}$

上面二式相減：

$\begin{aligned} { \Delta y} &= \beta_1 (\ln (X{+ \Delta x})-\ln X)\\\Delta y &= \beta_1 \ln ({ \frac{X + \Delta x}{X}})=\beta_1 \ln(1+\frac{\Delta x}{X})\\\Delta y &\approx \beta_1 \frac{\Delta x}{X} \\
{ \beta_1} &{ \approx }{ \frac{\Delta y}{\frac{\Delta x}{X}}}\end{aligned}$

> **$\beta_1$：當 $X$ 變動 1% 時，$Y$ 會隨之變動 $0.01\beta_1$ 單位。**

2. **Log-linear** model
$\begin{aligned}\ln Y = \beta_0 + \beta_1 X+e \end{aligned}$

> $\beta_1$ 如何解釋？
>

考慮 $X$ 微小變動對 $Y$ 的影響。

$\begin{aligned}\begin{aligned}\ln (Y { + \Delta y})  = \beta_0 + \beta_1 (X{+ \Delta x})+e\end{aligned}\end{aligned}$

上面二式相減：

$$\begin{aligned}
\ln (Y { + \Delta y}) - \ln Y &= \beta_1 \Delta X \\
\ln(1+\frac{\Delta y}{Y}) &= \beta_1 \Delta X \\
\frac{\Delta y}{Y} &\approx \beta_1 \Delta X \\
{ \beta_1} &{ \approx} { \frac{\frac{\Delta y}{Y}}{\Delta X}}
\end{aligned}$$

> **$\beta_1$：當 $X$ 變動 1單位時，$Y$ 會隨之變動 $100\beta_1 \%$ 單位。**

</aside>

1. **Log-log** model
$\begin{aligned}\ln Y = \beta_0 + \beta_1 \ln X+e \end{aligned}$

> $\beta_1$ 如何解釋？
>

考慮 $X$ 微小變動對 $Y$ 的影響。

$\begin{aligned}\begin{aligned}\ln (Y { + \Delta y}) = \beta_0 + \beta_1 \ln (X{+ \Delta x})+e\end{aligned}\end{aligned}$

上面二式相減：

$\begin{aligned}  \ln(Y{ + \Delta y})-\ln Y &= \beta_1[\ln(X{ +\Delta x})-\ln X] \\\ln(\frac{Y+ \Delta y}{Y}) &= \beta_1 \ln(\frac{X+ \Delta x}{X})\\\ln(1+\frac{ \Delta y}{Y}) &= \beta_1 \ln(1+ \frac{ \Delta x}{X}) \\\frac{ \Delta y}{Y} & \approx \beta_1 \frac{ \Delta x}{X}\\{ \beta_1} & { \approx} { \frac{\frac{ \Delta y}{Y}}{\frac{ \Delta x}{X}}}
\end{aligned}$

> **$\beta_1$：當 $X$ 變動 1% 時，$Y$ 會隨之變動 $\beta_1 \%$ 。這也就是經濟學中所說的「彈性」。**

## 14.6.2 極端值處理

我們知道，線性迴歸分析是希望找出一條能夠代表資料「**平均**」的線，那既然是跟平均的概念有關，那就會與極端值處理有所關係。

#### 庫克距離 Cook's Distance

**庫克距離** (Cook's Distance) 是由 R. Dennis Cook 於 1977 年提出的一種基於「刪除法 (Deletion Diagnostics)」的全局影響力測度 (Global Influence Measure)。其核心概念為: 「衡量當我們從樣本集中移除第 $i$ 個觀測值時，模型所有適配值 (Fitted Values) 所發生的綜合性變動幅度。」
我們定義，第 $i$ 個觀測值的庫克距離記為 $D_i$。

我們以 `engel` 的原模型為例，畫出變數的庫克距離:

plot(ols_model, which = 4)

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/S1cuWRgtZl.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;"></span>
</div>

而根據資料分析實務的[經驗法則](https://en.wikipedia.org/wiki/Cook%27s_distance#cite_note-8)，我們可以優先刪除**庫克距離**$D_i > 1$ 的資料點。若我們採取[嚴格標準](https://taylorandfrancis.com/knowledge/Engineering_and_technology/Engineering_support_and_special_topics/Cook%27s_distance/)，則可以採用$D_i > 4/n$ ( $n$ 為樣本數)。

我們以原資料以及先前進行過轉換的資料為例，來觀察移除極端值對資料的影響:

# 宣告空清單與門檻
model_list <- list()                # 先準備一個空的清單來裝模型
threshold <- 4/nrow(engel)          # 設定 Cook's Distance 門檻

for (y in c('foodexp', 'foodexp_sqrt', 'foodexp_bc')){
  
  cat("\n========== 現在處理變數:", y, "==========\n")
  
  # 利用動態公式建立模型
  formula_str <- paste0(y, " ~ income")
  model <- lm(as.formula(formula_str), data = engel)
  
  # 算出 Cook's Distance 並抓取極端值
  cooks_d <- cooks.distance(model)
  bad_apples <- which(cooks_d > threshold)
  
  cat("抓到", length(bad_apples), "筆極端值！\n")
  
  # 防呆設置
  if (length(bad_apples) > 0) {
    print(bad_apples)
    engel_clean <- engel[-bad_apples, ] # 有bad_apples，執行剔除
  } else {
    engel_clean <- engel                # 沒bad_apples，沿用原資料
  }
  
  cat("清理前樣本數:", nrow(engel), "清理後樣本數:", nrow(engel_clean), "\n")
  
# 將模型存入樣本
  model_list[[y]] <- lm(as.formula(formula_str), data = engel_clean)
}

觀察移除極端值後的模型，是否有更好的表現

# 畫出 qqplot
par(mfrow = c(1, 3))

# 原模型
plot(model_list[['foodexp']], which = 2,
     main = "Original dataset")

# 使用 Box-Cox 轉換最佳 lambda    
plot(model_list[['foodexp_bc']], which = 2,
     main = "Box-Cox transformation, λ =0.636")

# 開根號
plot(model_list[['foodexp_sqrt']], which = 2,
     main = "Squared root transformation")


# 畫出 Residuals vs Fitted
par(mfrow = c(1, 3))

# 原模型
plot(model_list[['foodexp']], which = 1,
     main = "Original dataset")

# 使用 Box-Cox 轉換最佳 lambda    
plot(model_list[['foodexp_bc']], which = 1,
     main = "Box-Cox transformation, λ =0.636")

# 開根號
plot(model_list[['foodexp_sqrt']], which = 1,
     main = "Squared root transformation")


# 14.7 解決異質變異數問題

我們透過前面的學習，已知道在進行迴歸分析時，常常會遇到異質變異數的問題，其中一個作法是對資料進行轉換，但這可能會使資料失去經濟意涵，因此我們為何不去直接接受異質變異數呢?

**Robust 穩健變異數矩陣公式 (HC0)：**$$\text{Var}(\hat{\beta}) = (X^T X)^{-1} (X^T \Omega X) (X^T X)^{-1}$$

# 載入需要的套件，請事先下載
library(lmtest)   # 為了使用 coeftest()
library(sandwich) # 為了使用 vcovHC() 也就是 HC0/White 檢定

ols_model <- lm(foodexp ~ income, data = engel)

# 此為原本的 OLS 報表，方便等一下對比
cat("--- 原始 OLS 報表 (未校正) ---\n")
printCoefmat(summary(ols_model)$coefficients)

cat("\n--- 穩健標準誤報表 (HC0 校正後) ---\n")
# 2. 執行 Robust
robust_result <- coeftest(ols_model, 
                          df = Inf, # 樣本數夠大，直接用常態分配 (z分配) 逼近
                          vcov = vcovHC(ols_model, type = "HC0"))
print(robust_result)

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/HJyN6pxtZx.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">使用原版的假設在異質變異數下的檢定的顯著性值得存疑</span>
</div>

<div style="text-align: center;">
<img src="https://hackmd.io/_uploads/BJek9aTgYbe.png" width="450">
<br>
<span style="color: gray; font-size: 0.8em;">使用穩健的變異數，可以讓我們在進行係數檢定時更有信心</span>
</div>

# 參考資料
1. 陳旭昇(2024)，資料分析的統計學基礎:使用R語言，東華書局
2. 林建甫 Jeff Lin(2020)，[R 資料科學與統計](https://bookdown.org/jefflinmd38/r4biost)
3. Gareth James, Daniela Witten, Trevor Hastie, Robert Tibshirani. (2017). An Introduction to Statistical Learning: With Applications in R. New York: Springer.
4. 陳基國(2024). 基礎統計與R語言. 台北：五南圖書出版股份有限公司
5. [geeksforgeeks(2025) Hyperplane, Subspace and Halfspace](https://www.geeksforgeeks.org/engineering-mathematics/hyperplane-subspace-and-halfspace/)