# 多元線性迴歸
{%hackmd @themes/orangeheart %}
###### tags: `Stats`
回顧一下簡單線性迴歸，我們觀測到有一組 $y$ 與 $x$：
$$
y_{i}=\beta_{0}+\beta_{1} x_{i}+\epsilon_{i}, i=1,2, \ldots, n
$$
若要估計 $\beta_0$ 與 $\beta_1$，我們會利用**最小平方法(ordinary least squares, OLS)：**
$$
\min _{\beta_{0}, \beta_{1}} \sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\beta_{1} x_{i}\right)^{2}
$$
根據一階條件，我們可以得到：
$$
\hat{\beta}_{1}=\frac{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}}\; \text {, and } \hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1} \bar{x}
$$
#### 🖋️ 飢寒起盜心？2001年至2018年的失業率與竊盜案的關係
:open_file_folder: [theft.csv](https://cool.ntu.edu.tw/courses/14061/files/1843652/download?download_frd=1)
```julia
Dat<-read.csv("theft.csv")
lm(theft~umemploy_ppl, data = Dat)
##
## Call:
## lm(formula = theft ~ umemploy_ppl, data = Dat)
##
## Coefficients:
## (Intercept) umemploy_ppl
## 19149.284 -7.783
```
在上面的例子中，我們得到了一個關於竊盜案與失業率的線性關係：
$$
\texttt{theft}_t = \beta_0 + \beta_1 \texttt{unemploy_ppl}_t + \varepsilon_t,\; t= 1,2,\ldots, 220
$$
但這個線性關係似乎無法良好解釋這個現象。
## 多元線性迴歸的基本觀念
針對上述的例子，另一個想法是我們可以增加其他的解釋變數：
$$
y_{i}=\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{K} x_{i K}+\epsilon_{i},\; i=1,2, \ldots, n
$$
同樣利用最小平方法，我們也可以估計$\beta_{0}, \beta_{1}, \beta_{2}, \ldots, \beta_K$，若要解釋 $\beta_1$，代表以平均來說，假設其他解釋變數不變的情況之下，$X_1$每增加一單位，$Y$ 會增加 $\beta_1$ 單位。考慮以下條件，
$$
\min _{\beta_{0}, \beta_{1}, \beta_{2}, \ldots, \beta_K} \sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-\beta_{2} x_{i 2}-\cdots-\beta_{K} x_{i K}\right)^{2}
$$
而 $K+1$ 的一階條件如下：
$$
\begin{aligned}&-2 \sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-\beta_{2} x_{i 2}-\cdots-\beta_{K} x_{i K}\right)=0 \\&-2 \sum_{i=1}^{n} x_{i 1}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-\beta_{2} x_{i 2}-\cdots-\beta_{K} x_{i K}\right)=0 \\&-2 \sum_{i=1}^{n} x_{i 2}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-\beta_{2} x_{i 2}-\cdots-\beta_{K} x_{i K}\right)=0 \\&\qquad\;\vdots \\&-2 \sum_{i=1}^{n} x_{i K}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-\beta_{2} x_{i 2}-\cdots-\beta_{K} x_{i K}\right)=0\end{aligned}
$$
假設 $K =2$，我們對於 $\beta_0$、$\beta_1$ 與 $\beta_2$ 進行估計，根據一階條件，可以得到：
$$
\begin{gathered}
{\left[\sum_{i}\left(x_{i 2}-\bar{x}_{2}\right)^{2}\right]\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)\left(y_{i}-\bar{y}\right)\right]-} \\
\hat{\beta}_{1}=\frac{\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)\left(x_{i 2}-\bar{x}_{2}\right)\right]\left[\sum_{i}\left(x_{i 2}-\bar{x}_{2}\right)\left(y_{i}-\bar{y}\right)\right]}{\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)^{2}\right]\left[\sum_{i}\left(x_{i 2}-\bar{x}_{2}\right)^{2}\right]} \\
-\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)\left(x_{i 2}-\bar{x}_{2}\right)\right]^{2} \\
\hat{\beta}_{2}=\frac{\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)\left(x_{i 2}-\bar{x}_{2}\right)\right]\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)\left(y_{i}-\bar{y}\right)\right]}{\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)^{2}\right]\left[\sum_{i}\left(x_{i 2}-\bar{x}_{2}\right)^{2}\right]} \\
-\left[\sum_{i}\left(x_{i 1}-\bar{x}_{1}\right)\left(x_{i 2}-\bar{x}_{2}\right)\right]^{2} \\
\hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1} \bar{x}_{1}-\hat{\beta}_{2} \bar{x}_{2} .
\end{gathered}
$$
但會發現到，這樣的結果似乎有點複雜，甚至難以理解。此時，我們可以透過線性代數幫我們解決這個~~世紀難題~~複雜的公式。令 $\mathbf{y}$、$\mathbf{X}$ 與 $\varepsilon$如下：
$$
\mathbf{y}=\left[\begin{array}{c}y_{1} \\y_{2} \\\vdots \\y_{n}\end{array}\right], \mathbf{X}=\left[\begin{array}{ccccc}1 & x_{11} & x_{12} & \cdots & x_{1 K} \\1 & x_{21} & x_{22} & \cdots & x_{2 K} \\\vdots & \vdots & \vdots & & \vdots \\1 & x_{n 1} & x_{n 2} & \cdots & x_{K}\end{array}\right], \epsilon=\left[\begin{array}{c}\epsilon_{1} \\\epsilon_{2} \\\vdots \\\epsilon_{n}\end{array}\right]
$$
而斜率與截距可以寫為
$$
\beta=\left[\begin{array}{c}\beta_{0} \\\beta_{1} \\\vdots \\\beta_{K}\end{array}\right]
$$
因此，我們可以寫出由 $\mathbf{y}$、$\mathbf{X}$ 與 $\varepsilon$ 組成的線性方程：
$$
\mathbf{y} = \mathbf{X} \beta + \varepsilon.
$$
我們可以把 $\mathbf{X} \beta$ 移到左邊，得到
$$
\begin{gathered}\mathbf{y}-\mathbf{X} \beta=\left[\begin{array}{c}y_{1} \\y_{2} \\\vdots \\y_{n}\end{array}\right]-\left[\begin{array}{ccccc}1 & x_{11} & x_{12} & \cdots & x_{1 K} \\1 & x_{21} & x_{22} & \cdots & x_{2 K} \\\vdots & \vdots & \vdots & & \vdots \\1 & x_{n 1} & x_{n 2} & \cdots & x_{K}\end{array}\right]\left[\begin{array}{c}\beta_{0} \\\beta_{1} \\\vdots \\\beta_{K}\end{array}\right] \\=\left[\begin{array}{c}y_{1}-\beta_{0}-\beta_{1} x_{11}-\beta_{2} x_{12}-\cdots-\beta_{K} x_{1 K} \\y_{2}-\beta_{0}-\beta_{1} x_{21}-\beta_{2} x_{22}-\cdots-\beta_{K} x_{2 K} \\\vdots \\y_{n}-\beta_{0}-\beta_{1} x_{n 1}-\beta_{2} x_{n 2}-\cdots-\beta_{K} x_{n K}\end{array}\right]\end{gathered}
$$
而為了做最小平方法，我們可以對於 $\mathbf{y} - \mathbf{X}\beta$ 這個矩陣進行轉置並相乘，就可以得到平方的形式：
$$
(\mathbf{y}-\mathbf{X} \beta)^{\top}(\mathbf{y}-\mathbf{X} \beta)=\sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-\beta_{2} x_{i 2}-\cdots-\beta_{K} x_{i K}\right)^{2}\\
$$
經過計算，可以得到
$$
\mathbf{y}^{\top} \mathbf{y}-\beta^{\top} \mathbf{X}^{\top} \mathbf{y}-\mathbf{y}^{\top} \mathbf{X} \beta+\beta^{\top} \mathbf{X}^{\top} \mathbf{X} \beta
$$
若要估計 $\beta$ ，我們則須對以下項目取一階條件：
$$
\nabla_{\beta} \beta^{\top} \mathbf{X}^{\top} \mathbf{y}=\nabla_{\beta} \mathbf{y}^{\top} \mathbf{X} \beta=\mathbf{X}^{\top} \mathbf{y}, \text { and } \nabla_{\beta} \beta^{\top} \mathbf{X}^{\top} \mathbf{X} \beta=2 \mathbf{X}^{\top} \mathbf{X} \beta
$$
得到
$$
-2 \mathbf{X}^{\top} \mathbf{y}+2 \mathbf{X}^{\top} \mathbf{X} \hat{\beta}=\mathbf{0}
$$
若 $\mathbf{X}^{\top} \mathbf{X}$ 是可反轉的，我們可以得到估計後的 $\hat{\beta}$：
$$
\hat{\beta}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y}
$$
當 $K=1$ 時，
$$
\mathbf{X}^{\top} \mathbf{X}=\left[\begin{array}{cc}n & \sum_{i} x_{i} \\\sum_{i} x_{i} & \sum_{i} x_{i}^{2}\end{array}\right], \text { and } \mathbf{X}^{\top} \mathbf{y}=\left[\begin{array}{c}\sum_{i} y_{i} \\\sum_{i} x_{i} y_{i}\end{array}\right]
$$
因為
$$
\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}=\frac{1}{n \sum_{i} x_{i}^{2}-\left(\sum_{i} x_{i}\right)^{2}}\left[\begin{array}{cc}\sum_{i} x_{i}^{2} & -\sum_{i} x_{i} \\-\sum_{i} x_{i} & n\end{array}\right]
$$
我們可以得到
$$
\hat{\beta}=\left[\hat{\beta}_{0} \;\hat{\beta}_{1}\right]^{\top}=\frac{1}{n \sum_{i} x_{i}^{2}-\left(\sum_{i} x_{i}\right)^{2}}\left[\begin{array}{c}\left(\sum_{i} x_{i}^{2}\right)\left(\sum_{i} y_{i}\right)-\left(\sum_{i} x_{i}\right)\left(\sum_{i} x_{i} y_{i}\right) \\n\left(\sum_{i} x_{i} y_{i}\right)-\left(\sum_{i} x_{i}\right)\left(\sum_{i} y_{i}\right)\end{array}\right]
$$
因此
$$
\hat{\beta}_{1}=\frac{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}}, \text { and } \hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1} \bar{x}
$$
注意到，在 $K=1$ 時，我們僅要求 $x_i$ 的每個值都不相同，或者說 $\sum_{i}\left(x_{i}-\bar{x}\right)^{2}=0$；總的來說，若 $K= k$，我們必須確定 $\hat{\beta} = \mathbf{X}^{\top}\mathbf{X}$ 是可反轉的。我們會說若$(\mathbf{X}^{\top}\mathbf{X})^{-1}$存在的條件是其必須為滿秩(full rank)的，也就是 $\mathbf{X}$ 的每個欄位**不存在線性關係**，每個解釋變數都有各自的解釋力，彼此互不重合。$\mathbf{X}$ 不是一個滿秩的矩陣，我們稱其為一個共線性(multicollinearity)或**線性重合**的矩陣，此時 $\mathbf{X}^{\top}\mathbf{X}$ 不可進行反轉，最小平方法也就失去效力了。若當 $K=1$ 時，
$$
\hat{y} = \hat{\beta_0}+\hat{\beta_1}x
$$
我們稱其為**迴歸線(regression line)**。若 $K=2$ 時，
$$
\hat{y}=\hat{\beta}_{0}+\hat{\beta}_{1} x_{1}+\hat{\beta}_{2} x_{2}
$$
其稱為**迴歸平面(regression plane)**。若 $K$ 為任何數，則
$$
\hat{y}=\hat{\beta}_{0}+\hat{\beta}_{1} x_{1}+\hat{\beta}_{2} x_{2}+\cdots+\hat{\beta}_{K} x_{K}
$$
其稱為**超迴歸平面(regression hyperplane)**，無法在目前所知的三維空間繪圖。
## 代數上的性質
在 $K +1$ 個一階條件的情況下，
$$
\mathbf{X}^{\top}(\mathbf{y}-\mathbf{X} \hat{\beta})=\mathbf{X}^{\top} \hat{\epsilon}=\left[\begin{array}{c}\sum_{i=1}^{n} \hat{\epsilon}_{i} \\\sum_{i=1}^{n} x_{i 1} \hat{\epsilon}_{i} \\\sum_{i=1}^{n} x_{i 2} \hat{\epsilon}_{i} \\\vdots \\\sum_{i=1}^{n} x_{i K} \hat{\epsilon}_{i}\end{array}\right]=\mathbf{0}
$$
也就是說，$\sum_{i=1}^{n} \hat{\epsilon}_{i}=\sum_{i=1}^{n} x_{i 1} \hat{\epsilon}_{i}=\cdots=\sum_{i=1}^{n} x_{i K} \hat{\epsilon}_{i}=0$。又因為 $\hat{\mathbf{y}}=\mathbf{X} \hat{\beta}$，因此
$$
\sum_{i=1}^{n} \hat{y}_{i} \hat{\epsilon}_{i}=\hat{\mathbf{y}}^{\top} \hat{\epsilon}=\hat{\beta}^{\top} \mathbf{X}^{\top} \hat{\epsilon}=0.
$$
## 配適值
與簡單線性迴歸一樣，我們也可以利用 $S S T=S S R+S S E$ 計算配適值：
$$
\sum_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}=\sum_{i=1}^{n}\left(\hat{y}_{i}-\bar{y}+\hat{\epsilon}_{i}\right)^{2}=\sum_{i=1}^{n}\left(\hat{y}_{i}-\bar{y}\right)^{2}+\sum_{i=1}^{n} \hat{\epsilon}_{i}^{2}
$$
- **SST: $\sum_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}$**
- **SSR: $\sum_{i=1}^{n} \hat{\varepsilon}_{i}^{2}$**
- **SSE: $\sum_{i=1}^{n}\left(\hat{y}_{i}-\bar{y}\right)^{2}$**
$$
0\leq R^{2}=\frac{S S E}{S S T}=1-\frac{S S R}{S S T} \leq 1
$$
當 $SSE$ 佔 $SST$ 的比例越高，代表其配適度越好。一個值得注意的性質是，當我們放越多解釋變數，$R^2$ 會越高，但實際上是這樣嗎？考慮以下兩個模型：
$$
\begin{aligned}&y_{i}=\beta_{0}+\beta_{1} x_{i 1}+\epsilon_{i} \\&y_{i}=\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\epsilon_{i}\end{aligned}
$$
則
$$
\begin{aligned}&SSR_{1}=\min _{\beta_{0}, \beta_{1}} \sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-0 \cdot x_{i 2}\right)^{2} \\&\geq \min _{\beta_{0}, \beta_{1}, \beta_{2}} \sum_{i=1}^{n}\left(y_{i}-\beta_{0}-\beta_{1} x_{i 1}-\beta_{2} x_{i 2}\right)^{2}=SSR_{2} .\end{aligned}
$$
在第一個模型中，我們特別規定 $\beta_2 = 0$，第二個模型則無，這個限制將會讓第一個模型的 $SSR$ 大於第二個模型的 $SSR$（可以想像成是讓它們打架，但第一個模型被綁手綁腳，僅能在 $0$ 或 $1$ 之間進行選擇）。但這就變成研究者必須進行取捨：要放很少的解釋變數讓模型變簡單，抑或是放很多解釋變數讓 $R^2$ 變大，但模型變複雜。有個解決上述問題的方法，利用調整後的 $R^2$(adjusted $R^2$)：
$$
\bar{R}^{2}=1-\frac{S S R /(n-K-1)}{S S T /(n-1)}
$$
或寫作
$$
\bar{R}^{2}=R^{2}-\frac{K}{n-K-1}\left(1-R^{2}\right)
$$
其中 $\frac{K}{n-K-1}$ 稱為懲罰項， $K$ 則代模型的複雜度。這個方法某種程度上就可以對於模型的配適度與複雜度進行平衡。白話來說，$\bar{R^2}$ 為控制解數變數個數與樣本數下，迴歸模型中可解釋變異占總變異的比 例，即考慮了模型中放太多**垃圾變數（沒有解釋力）**的情形，在計算時給予懲罰項，其下界可能小於0，但上界還是1。迴歸分析的主要目的不是極大化 $R^2$、$\bar{R^2}$，而是探討 $X$ 和 $Y$ 之間是否具有特定關係，而使用多元迴歸的模型，是為了避免遺漏變數（違反外生性）等相關問題而引發許多統計推論上的不便。
## 統計上的性質：古典條件二
被解釋變數 $\{y_i\}_{i=1,2,\ldots,n}$ 使得
- $\mathbb{E}\left(y_{i}\right)=\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{K} x_{i K}$，其中 $\left\{x_{i 1}, x_{i 2}, \ldots, x_{i K}\right\}$為非隨機的
- $\operatorname{var}(y_i) = \sigma^2_0$，且 $\operatorname{cov}(y_i,y_j)=0,\; \forall i \neq j$。
上述假設等同於
$$
y_{i}=\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{k} x_{i K}+\epsilon_{i}
$$
其中 $\left\{x_{i 1}, x_{i 2}, \ldots, x_{i K}\right\}$ 為非隨機的、 $\mathbb{E}(\varepsilon_i) = 0$ 、$\operatorname{var}\left(\varepsilon_{i}\right)=\sigma_{0}^{2}$ 且 $\operatorname{cov}\left(\epsilon_{i}, \epsilon_{j}\right)=0,\; \forall i \neq j$。在上述的假設條件下，$\hat{\beta}$ 是不偏的，也就是 $\mathbb{E}(\hat{\beta}) = \beta$。因為$\mathbb{E}\left(y_{i}\right)=\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{K} x_{i K}$，得到 $\mathbb{E}(\mathbf{y})=\mathbf{X} \beta$，
$$
\begin{aligned}\mathbb{E}(\hat{\beta}) &=\mathbb{E}\left[\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y}\right]=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbb{E}(\mathbf{y}) \\&=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{X} \beta=\beta\end{aligned}
$$
又
$$
\begin{aligned}\hat{\beta}-\beta &=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y}-\beta=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top}(\mathbf{X} \beta+\varepsilon)-\beta \\&=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{X} \beta+\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon-\beta=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon\end{aligned}
$$
在條件二之下，
$$
\mathbb{E}(\varepsilon\varepsilon^{\top}) = \mathbb{E}\left\{\left[\begin{array}{cccc}
\varepsilon_{1}^{2} & \varepsilon_{1} \varepsilon_{2} & \cdots & \varepsilon_{1} \varepsilon_{n} \\
\varepsilon_{1} \varepsilon_{2} & \varepsilon_{2}^{2} & \ddots & \vdots \\
\vdots & \ddots & \ddots & \varepsilon_{n-1} \varepsilon_{n} \\
\varepsilon_{1} \varepsilon_{n} & \cdots & \varepsilon_{n-1} \varepsilon_{n} & \varepsilon_{n}^{2}
\end{array}\right]\right\}=\left[\begin{array}{cccc}
\sigma^{2} & 0 & \cdots & 0 \\
0 & \sigma^{2} & \ddots & \vdots \\
\vdots & \ddots & \ddots & 0 \\
0 & \cdots & 0 & \sigma^{2}
\end{array}\right]
$$
而 $\hat{\beta}$ 的變異數矩陣可以寫作：
$$
\begin{aligned}\operatorname{cov}(\hat{\beta}) &=\mathbb{E}\left[(\hat{\beta}-\beta)(\hat{\beta}-\beta)^{\top}\right] \\&=\mathbb{E}\left[\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon \varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right] \\&=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbb{E}\left(\varepsilon \epsilon^{\top}\right) \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \\&=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \sigma^{2} \mathbf{I}_{n} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \\&=\sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \\&=\sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\end{aligned}
$$
當 $K=1$ 時，
$$
\operatorname{cov}(\hat{\beta})=\sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}=\frac{\sigma^{2}}{n \sum_{i} x_{i}^{2}-\left(\sum_{i} x_{i}\right)^{2}}\left[\begin{array}{cc}\sum_{i} x_{i}^{2} & -\sum_{i} x_{i} \\-\sum_{i} x_{i} & n\end{array}\right]
$$
也就是
$$
\operatorname{var}\left(\hat{\beta}_{1}\right)=\frac{n \sigma^{2}}{n \sum_{i} x_{i}^{2}-\left(\sum_{i} x_{i}\right)^{2}}=\frac{\sigma^{2}}{\sum_{i}\left(x_{i}-\bar{x}\right)^{2}}
$$
$\sigma^2$ 的 OLS 代表殘差平方的平均：
$$
s^{2}=\frac{1}{n-K-1} \sum_{i=1}^{n} \hat{\varepsilon}_{i}^{2}=\frac{1}{n-K-1} \hat{\varepsilon}^{\top} \hat{\varepsilon}
$$
但為什麼是 $n-K-1$ 呢？因為
$$
\hat{\varepsilon}=\mathbf{y}-\mathbf{X} \hat{\beta}=\mathbf{y}-\mathbf{X} \beta-\mathbf{X}(\hat{\beta}-\beta)=\varepsilon-\mathbf{X}(\hat{\beta}-\beta)
$$
因此，
$$
\begin{aligned}\mathbb{E}\left(\hat{\varepsilon}^{\top} \hat{\varepsilon}\right)=& \mathbb{E}\left\{[\varepsilon-\mathbf{X}(\hat{\beta}-\beta)]^{\top}[\varepsilon-\mathbf{X}(\hat{\beta}-\beta)]\right\} \\=& \mathbb{E}\left(\varepsilon^{\top} \varepsilon\right)-2 \mathbb{E}\left[\varepsilon^{\top} \mathbf{X}(\hat{\beta}-\beta)\right] \\&+\mathbb{E}\left[(\hat{\beta}-\beta)^{\top} \mathbf{X}^{\top} \mathbf{X}(\hat{\beta}-\beta)\right] .\end{aligned}
$$
由於
$$
\varepsilon^{\top} \mathbf{X}(\hat{\beta}-\beta)=\varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon
$$
且
$$
\begin{aligned}&(\hat{\beta}-\beta)^{\top} \mathbf{X}^{\top} \mathbf{X}(\hat{\beta}-\beta) \\&=\varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon \\&=\varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon,\end{aligned}
$$
我們得到
$$
\mathbb{E}\left(\hat{\varepsilon}^{\top} \hat{\varepsilon}\right)=\mathbb{E}\left(\varepsilon^{\top} \varepsilon\right)-\mathbb{E}\left[\varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon\right]
$$
而 $\mathbb{E}\left(\epsilon^{\top} \epsilon\right)=n \sigma^{2}$，
$$
\begin{aligned}& \mathbb{E}\left[\varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon\right]=\mathbb{E}\left\{\text { trace }\left[\varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \varepsilon\right]\right\} \\=& \mathbb{E}\left\{\operatorname{trace}\left[\mathbf{X}^{\top} \varepsilon \varepsilon^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right]\right\}=\operatorname{trace}\left[\mathbf{X}^{\top} \mathbb{E}\left(\varepsilon \varepsilon^{\top}\right) \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right.\\=& \operatorname{trace}\left[\mathbf{X}^{\top} \sigma^{2} \mathbf{I} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right]=\sigma^{2} \operatorname{trace}\left[\mathbf{X}^{\top} \mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right] \\=& \sigma^{2} \text { trace }\left(\mathbf{I}_{K+1}\right)=(K+1) \sigma^{2} .\end{aligned}
$$
綜上，$\mathbb{E}\left(\hat{\varepsilon}^{\top} \hat{\varepsilon}\right)=(n-K-1) \sigma^{2}$，
$$
\mathbb{E}\left(\frac{1}{n-K-1} \hat{\varepsilon}^{\top} \hat{\varepsilon}\right)=\sigma^{2}
$$
```r
Dat <- read.csv("theft.csv")
summary(lm(theft~umemploy_ppl+clearance, data = Dat))
##
## Call:
## lm(formula = theft ~ umemploy_ppl + clearance, data = Dat)
##
## Residuals:
## Min 1Q Median 3Q Max
##-12613.4 -2285.3 -246.2 2013.6 18106.6
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65790.145 2813.975 23.380 <2e-16 ***
## umemploy_ppl -7.067 4.833 -1.462 0.145
## clearance -613.196 20.959 -29.256 <2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 4184 on 206 degrees of freedom
## Multiple R-squared: 0.8065, Adjusted R-squared: 0.8046
## F-statistic: 429.3 on 2 and 206 DF, p-value: < 2.2e-16
```
上述的所有性質可以合而為一成為一個定理，稱作Gauss-Markov定理：$\hat{\beta}$ 是**最好的線性不偏估計式(best linear unbiased estimator, BLUE)**。$\hat{\beta}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y}$，假設其他的線性不偏估計式 $\mathbf{Ay}$，即 $\mathbb{E}(\mathbf{Ay}) = \beta$。令 $\mathbf{A}-\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top}=\mathbf{C}$，也就是 $\mathbf{A}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top}+\mathbf{C}$，而 $\mathbf{A y}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y}+\mathbf{C y}=\hat{\beta}+\mathbf{C y}$。由於 $\mathbb{E}(\mathbf{Ay}) = \beta$，且 $\mathbf{A y}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y}+\mathbf{C y}=\hat{\beta}+\mathbf{C y}$，可以得到
$$
\mathbb{E}(\mathbf{C y})=\mathbb{E}(\mathbf{C X} \beta)+\mathbb{E}(\mathbf{C} \epsilon)=\mathbf{C X} \beta=\mathbf{0}
$$
隱含了 $C X=0$。$\mathbf{Ay}$ 的變異數矩陣為
$$
\operatorname{cov}(\hat{\beta}+\mathbf{C y})=\operatorname{cov}(\hat{\beta})+\operatorname{cov}(\mathbf{C y})+2 \mathbb{E}\left[(\hat{\beta}-\beta)(\mathbf{C y}-\beta)^{\top}\right]
$$
其中
$$
\begin{aligned}& \mathbb{E}\left[(\hat{\beta}-\beta)(\mathbf{C} \mathbf{y}-\beta)^{\top}\right]=\mathbb{E}\left[\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \epsilon \mathbf{y}^{\top} \mathbf{C}^{\top}\right] \\=&\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbb{E}\left(\epsilon \mathbf{y}^{\top}\right) \mathbf{C}^{\top}=\sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{C}^{\top}=\mathbf{0}\end{aligned}
$$
其證明了 $\operatorname{cov}(\mathbf{A y})-\operatorname{cov}(\hat{\beta})$ 為一個**半正定矩陣**，代表 $\hat{\beta}$ 之有效性較其他的不偏估計式還要來得高。令 $\omega = \omega_1$，其中，
$$
\omega_1 = \begin{bmatrix}0\\1\\0\\\vdots\\0\end{bmatrix}
$$
我們也可以用幾何學的角度來看這件事。對於兩個 $n \times 1$ 的向量 $\mathbf{X}$ 與 $\mathbf{z}$，其內積(inner product)為
$$
\mathbf{X}^{\top} \mathbf{z}=\sum_{i=1}^{n} x_{i} z_{i}
$$
另外，提供一個有用的工具：Euclidean norm，其能夠衡量向量的長度：
$$
\|\mathbf{X}\|=\sqrt{\mathbf{X}^{\top} \mathbf{X}}=\sqrt{\sum_{i=1}^{n} x_{i}^{2}}
$$
如果$\mathbf{x}^{\top}\mathbf{z}=0$，則我們稱這兩個矩陣相互垂直。此外，我們要利用投影矩陣(projection matrix)具有**冪等(idempotent)**的性質，也就是所謂的自乘不變，即
$$
\mathbf{A} \mathbf{A}=\mathbf{A}.
$$
假設 $\mathbf{A}$ 是一個可逆矩陣，即 $\mathbf{A}^{-1}$存在，則
$$
\begin{aligned}&\mathbf{A}\mathbf{A}= \mathbf{A}\\&\Rightarrow \mathbf{A}^{-1}\mathbf{A}\mathbf{A}=\mathbf{A}^{-1}\mathbf{A}\\&\Rightarrow \mathbf{A} = \mathbf{I}.\\\end{aligned}
$$
如果說一個矩陣的轉置矩陣與原本的矩陣相同，那我們便會說這個矩陣是**垂直的投影矩陣(orthogonal projection matrix)**，即
$$
\mathbf{A}= \mathbf{A}^{\top}
$$
而我們可利用這個性質，證明 $\mathbf{AX}$ 與 $\mathbf{(I-A)X}$ 是相互垂直的。對於兩個 $n \times n$ 的矩陣 $\mathbf{A}$ 與 $\mathbf{B}$ 而言，對於所有的 $\mathbf{X}$，若 $\mathbf{}^{\top} (\mathbf{A} - \mathbf{B})\mathbf{X} \geq 0$，則稱 $\mathbf{A} - \mathbf{B}$ 為一個半正定矩陣；而對於所有非零的 $\mathbf{X}$來說，若 $\mathbf{X}^{\top} (\mathbf{A} - \mathbf{B})\mathbf{X} > 0$，則稱 $\mathbf{A} - \mathbf{B}$ 為一個正定矩陣。令 $\mathbf{P} = \mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}$，我們可以將 $\mathbf{P}$ 視為是一個對稱(symmetric)、冪等(idemponent)，且為一個垂直投影的矩陣。因為 $\mathbf{\hat{y}=X}\hat{\beta}$，且$\hat{\beta}=\left(\mathbf{X}^{\top}\mathbf{X}\right)\mathbf{X}^{\top}\mathbf{y}$，那麼，配適值 $\mathbf{\hat{y}}$ 便可以寫成
$$
\hat{\mathbf{y}}=\mathbf{P y}=\mathbf{X}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y}
$$
其為$\mathbf{y}$ 的垂直投影向量，投影到由向量$\mathbf{X}$各欄位加總的集合，即 $\operatorname{span}(\mathbf{X})$。
而這個方法可以幫我們很好地去估計$\mathbf{y}$：考慮到殘差向量(residual vector) $\hat{\varepsilon}=(\mathbf{I}-\mathbf{P}) \mathrm{y}$ 的歐氏距離(Euclidean distance)比其他的 $(\mathbf{I}-\mathbf{A}) \mathbf{y}$ 都還要小，即
$$
\|\hat{\varepsilon}\|=\|(\mathbf{I}-\mathbf{P}) \mathbf{y}\| \leq\|(\mathbf{I}-\mathbf{A}) \mathbf{y}\|
$$
這正是OLS的本質阿！而殘差向量是 $\operatorname{span}(\mathbf{X})$ 的補集，記作 $\operatorname{span}(\mathbf{X})^{\perp}$，那麼便可得知其亦垂直於配適值向量：
$$
\hat{\mathbf{y}}^{\top} \hat{\epsilon}=\mathbf{y}^{\top} \mathbf{P}^{\top}(\mathbf{I}-\mathbf{P}) \mathbf{y}=0
$$
## ****Frisch-Waugh 定理****
在古老的統計學發展的時代，人們沒有電腦，也沒有計算機，要怎麼進行迴歸分析呢？有一個想法是把解釋變數與其係數分割成兩塊，將一塊進行迴歸後的結果拿去跑另一塊的迴歸。這聽起來有點複雜跟抽象，我們來看看他的運作。回憶到以下幾個性質：
$$
\begin{aligned}
\mathbf{y} =& \mathbf{X}\beta+\varepsilon\\
\hat{\beta} =& \left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{y}\\
\mathbf{P} =& \mathbf{X}\left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\\
\hat{\varepsilon} =& \mathbf{I}-\mathbf{P}=\mathbf{y}-\hat{\mathbf{y}}
\end{aligned}
$$
關於線性迴歸式，我們可以改寫成：
$$
\begin{bmatrix}y_1\\y_2\\y_3\\\vdots\\y_n\end{bmatrix} = \begin{bmatrix}x_{1,1} & x_{1,2} & \cdots & x_{1,k} & x_{1,k+1} & \cdots & x_{1,K}\\
x_{2,1} & x_{2,2} & \cdots & x_{2,k} & x_{2,k+1} & \cdots & x_{2,K}\\
\vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots
\\
x_{n,1} & x_{n,2} & \cdots & x_{n,k} & x_{n,k+1} & \cdots &x_{n,K}\end{bmatrix}
\begin{bmatrix}\beta_1\\\beta_2\\\vdots\\\beta_k\\\beta_{k+1}\\\vdots\\\beta_{K}\end{bmatrix}+\begin{bmatrix}\varepsilon_1\\\varepsilon_2\\\varepsilon_3\\\vdots\\\varepsilon_{K}\end{bmatrix}
$$
現在我們令 $\mathbf{X} = \left(\mathbf{x_1}, \mathbf{x_2}\right)$、$\mathbf{\beta} = \left(\mathbf{\beta_1}, \mathbf{\beta_2}\right)^{\top}$，而此時我們要重新對於迴歸式取OLS，即
$$
\min_{\beta}(\mathbf{y}-\mathbf{X} \beta)^{\top}(\mathbf{y}-\mathbf{X} \beta)
$$
根據一階條件，我們可以得到
$$
\mathbf{X}^{\top}\mathbf{X}\hat{\beta}=\mathbf{X}^{\top}\mathbf{y}
$$
代入我們上述的假設，可以將上式左邊寫作
$$
\left(\mathbf{x_1}, \mathbf{x_2}\right)^{\top}\left(\mathbf{x_1}, \mathbf{x_2}\right) = \begin{bmatrix}\mathbf{x_1}^{\top}\mathbf{x_1} & \mathbf{x_1}^{\top}\mathbf{x_2}\\\mathbf{x_2}^{\top}\mathbf{x_1} & \mathbf{x_2}^{\top}\mathbf{x_2}\end{bmatrix}
$$
右邊則可以寫作
$$
\mathbf{X}^{\top}\mathbf{y} = \begin{bmatrix}\mathbf{x_1}^{\top}\mathbf{y}\\\mathbf{x_2}^{\top}\mathbf{y}\end{bmatrix}
$$
故上式又可以寫成一個類似二元一次聯立方程式的形式：
$$
\mathbf{X}^{\top}\mathbf{X}\hat{\beta}=\mathbf{X}^{\top}\mathbf{y}\Rightarrow
\begin{cases}
\mathbf{x_1}^{\top}\mathbf{x_1}\hat{\mathbf{\beta_1}}+\mathbf{x_1}^{\top} \mathbf{x_2}\hat{\mathbf{\beta_2}}=\mathbf{x_1}^{\top}\mathbf{y}&\cdots\cdots(1)\\
\mathbf{x_2}^{\top}\mathbf{x_1}\hat{\mathbf{\beta_1}}+\mathbf{x_2}^{\top} \mathbf{x_2}\hat{\mathbf{\beta_2}}=\mathbf{x_2}^{\top}\mathbf{y}&\cdots\cdots(2)
\end{cases}
$$
根據 $(1)$ 我們可以得到 $\hat{\beta_1}$：
$$
\hat{\beta}_{1}=\left(\mathbf{x_1}^{\top} \mathbf{x_1}\right)^{-1}\left(\mathbf{x_1}^{\top} \mathbf{y}-\mathbf{x_1}^{\top} \mathbf{x_2} \hat{\beta}_{2}\right)=\left(\mathbf{x_1}^{\top} \mathbf{x_1}\right)^{-1} \mathbf{x_1}^{\top}\left(\mathbf{y}-\mathbf{x_2}\hat{\beta}_{2}\right) \;\;\; \cdots\cdots (3)
$$
我們可以將 $(3)$ 代入 $(2)$ 中，得到
$$
\left(\mathbf{x_2}^{\top} \mathbf{x_1}\right)\left(\mathbf{x_1}^{\top} \mathbf{x_1}\right)^{-1} \mathbf{x_1}^{\top}\left(\mathbf{y}-\mathbf{x_2} \hat{\beta}_{2}\right)+\mathbf{x_2}^{\top} \mathbf{x_2}\hat{\beta}_{2}=\mathbf{x_2}^{\top} \mathbf{y}
$$
接著可以發現到，我們可以令 $\mathbf{P} = \mathbf{x_1}\left(\mathbf{x_1}^{\top} \mathbf{x_1}\right)^{-1} \mathbf{x_1}^{\top}$，代表投影到 $\mathbf{x_1}$ 上。因此我們可以改寫上式：
$$
\begin{aligned}
&\mathbf{x_2}^{\top}\mathbf{P}\mathbf{y} - \mathbf{x_2}^{\top}\mathbf{P}\mathbf{x_2}\hat{\mathbf{\beta_2}} + \mathbf{x_2}^{\top}\mathbf{x_2}\hat{\mathbf{\beta_2}}=\mathbf{x_2}^{\top} \mathbf{y}\\
\Rightarrow \;&\left(\mathbf{x_2}^{\top}\mathbf{x_2}-\mathbf{x_2}^{\top}\mathbf{P}\mathbf{x_2}\right)\hat{\mathbf{\beta_2}}=\mathbf{x_2}^{\top} \mathbf{y}-\mathbf{x_2}^{\top}\mathbf{P} \mathbf{y}\\
\Rightarrow \;& \mathbf{x_2}^{\top}\left(\mathbf{x_2}-\mathbf{P}\mathbf{x_2}\right)\hat{\beta_2}=\mathbf{x_2}^{\top}\left(\mathbf{y}-\mathbf{P}\mathbf{y}\right)
\\
\Rightarrow \;& \mathbf{x_2}^{\top}\left(\mathbf{I}-\mathbf{P}\right)\mathbf{x_2}\hat{\mathbf{\beta_2}}=\mathbf{x_2}^{\top}\left(\mathbf{I}-\mathbf{P}\right)\mathbf{y}\\
\Rightarrow \;& \mathbf{x_2}^{\top}\varepsilon\mathbf{x_2}\hat{\mathbf{\beta_2}}=\mathbf{x_2}^{\top}\varepsilon\mathbf{y}
\end{aligned}
$$
根據殘差向量的垂直與冪等性質，即
$$
\begin{aligned}
\varepsilon=\varepsilon^{\top}\\
\varepsilon\varepsilon=\varepsilon
\end{aligned}
$$
可以得到
$$
\begin{aligned}
\hat{\beta_2} =& \left(\mathbf{x_2}^{\top}\varepsilon \mathbf{x_2}\right)^{-1}\mathbf{x_2}\varepsilon\mathbf{y}\\
=& \left(\mathbf{x_2}^{\top}\varepsilon\varepsilon\mathbf{x_2}\right)^{-1}\mathbf{x_2}\left(\varepsilon\varepsilon\right)\mathbf{y}\\
=& \big[\left(\varepsilon\mathbf{x_2}\right)^{\top}\varepsilon\mathbf{x_2}\big]\left(\varepsilon\mathbf{x_2}\right)^{\top}\varepsilon\mathbf{y}\\
=& \left(\tilde{X}^{\top} \; \tilde{X}\right)\tilde{X}\tilde{y}
\end{aligned}
$$
我們可以用下列這段原文做總結（不是我懶得翻譯，而是原文比較精確）。
> In the linear least squares regression of veotor $\mathbf{y}$ on two sets of var; ables, $\mathbf{X_{1}}$ and $\mathbf{X_{2}}$, the subvector $\hat{\beta}_{2}$ is the set of coefficients obtained when the residuals from a regression of $\mathbf{y}$ on $\mathbf{X_{1}}$ alone are regressed on the set of residuals obtained when each column of $\mathbf{X_{1}}$ is vegressed $\mathbf{X_{1}}$.
## 一些值得討論的迴歸分析方法
### 二元共變量(Binary Covariate)
假設我們的解釋變數不是 $0$就是 $1$，也就是僅有兩種類別呢？一個例子是，物理治療師(physiotherapist)對於他的 $12$ 名膝蓋受傷的病患進行追蹤，觀察他們術後是否能夠做超過 $40$ 分鐘的運動，若小於 $40$ 分鐘，則 $X = 0$；若大於 $40$ 分鐘則 $X = 1$，數據顯示如下：
```r
Y <- c(90, 65, 30, 60, 60, 80, 45, 45, 80, 35, 50, 45)
X <- c(0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0)
```
接著對其做迴歸，得到
```r
lm(Y ~ X)
##
## Call:
## lm(formula = Y ~ X)
##
## Coefficients:
## (Intercept) X
## 66.25 -27.50
```
我們可以寫出模型，
$$
Y = \beta_0 + \beta_1 X +\varepsilon
$$
其中 $X = 0$ 或 $1$。則
$$
\beta_{0}=\mathbb{E}(Y \mid X=0)=66.25, \; \beta_{1}=\mathbb{E}(Y \mid X=1)=38.75
$$
```r
mean(Y[X == 0])
## [1] 66.25
mean(Y[X == 1])
## [1] 38.75
```
接著繪圖可以發現：
```r
plot(X, Y, pch = 4, cex = 2, xlim = c(-1.0, 2.0))
abline(lm(Y~X), lwd = 3)
```
解釋變數都落在 $0$ 或 $1$ 上。如果我們對於竊盜率的那份資料也進行二元變數的分析呢？
```r
Dat<-read.csv("theft.csv")
summary(as.factor(Dat$recession))
## 0 1
## 163 46
```
得到了次數後，我們可以對其進行迴歸分析：
```r
lm(theft~recession, data = Dat)
##
## Call:
## lm(formula = theft ~ recession, data = Dat)
##
## Coefficients:
## (Intercept) recession
## 15463.41 -86.93
```
接著仿造上面的繪圖：
### 轉換共變量(Transformed Covariate)
不過，到現在我們可能還是沒搞清楚「線性」(linear)是什麼的線性。事實上，所謂線性，係指參數(parameters)對其（依變項）具有線性關係。例如
$$
Y =\beta_{0} + \beta_{1}X + \varepsilon
$$
就是線性模型。
$$
Y =\beta_{0} + \beta_{1}X^2 + \varepsilon
$$
也是線性模型。但是
$$
Y =\beta_{0} + \beta_{1}^2X + \varepsilon
$$
則是非線性模型(non-linear model)。同樣利用上面提到的物理治療師的例子，現在他改變追蹤方式，其觀察的解釋變數改為病人平均的運動時間（分鐘），被解釋變數改為觀察他們康復的時間。
```r
Y <- c(90, 65, 30, 60, 60, 80, 45, 45, 80, 35, 50, 45)
X <- c(24, 35, 64, 20, 33, 27, 42, 41, 22, 50, 36, 31)
```
考慮以下的模型：
$$
Y =\beta_{0} + \beta_{1}\sqrt{X} + \varepsilon
$$
而我們可以利用 `I()` 這個函數新建一個新變數：
```r
Y <- c(90, 65, 30, 60, 60, 80, 45, 45, 80, 35, 50, 45)
X <- c(24, 35, 64, 20, 33, 27, 42, 41, 22, 50, 36, 31)
lm(Y~I(sqrt(X)))
##
## Call:
## lm(formula = Y ~ I(sqrt(X)))
##
## Coefficients:
## (Intercept) I(sqrt(X))
## 145.77 -15.11
```
## 迴歸分析的相關議題
在美國，量測長度最常使用的單位是英寸、英尺（英制），而台灣則是以使用公分、公里（公制）居多。當我們想要對於 NBA 球員的身高、體重、練球年數$\cdots$ 與其在場上的表現進行迴歸分析時，會碰到一個狀況：我們雖說是蒐集完了資料，但由於這項報告是要給台灣的讀者以及 NBA 球迷看的，英制對於他們來說不是一個熟悉的單位制度，因此我們要對於量測單位進行轉換。給定迴歸模型 $Y = \beta_0 + \beta_1 X_i$，若我們對於原始資料進行線性轉換，我們就稱這樣的步驟稱為 **變數線性轉換**。
### 變數線性轉換基本概念
考慮一個簡單迴歸模型：
$$
Y = \beta_0 + \beta_1 X_i
$$
如果我們將原始資料 $X$ 與 $Y$ 進行線性轉換，那麼
$$
\begin{aligned}
X_i^* = a X_i + b\\
Y_i^* = cY_i + d
\end{aligned}
$$
則對於改寫過後的模型與原本的模型相比，其 OLS 估計式關係如下：
$$
\begin{aligned}
\hat{\beta_0} =&\; c\hat{\beta_0} + d - \frac{bc}{a}\hat{\beta_1} \\
\hat{\beta_1^*} =&\; \frac{c}{a}\hat{\beta}\\
\end{aligned}
$$
一般而言，線性轉換情況有(1)單位調整、(2)二元變數定義改變等等，例如：
- 去除平均(demean)：$(X_i^*, Y_i^*) = (X_i - \bar{X}, Y_i - \bar{Y} )$，則 $\hat{\beta_0} = 0$，而 $\hat{\beta_1^*} = \hat{\beta}$。
- 標準化(standardize)：$(X_i^*, Y_i^*) = (\frac{X_i - \bar{X}}{s_X}, \frac{Y_i - \bar{Y}}{s_Y} )$，則 $\hat{\beta_0} = 0$，而 $\hat{\beta_1^*} = \rho_{XY}$。
### 變數轉換技巧與步驟
以上述去除平均的例子說明。假設原本與新的模型之間的關係如下：
$$
\begin{aligned}
Y_i =& \; \beta_0 + \beta_1 X_i + \varepsilon_i\\
(Y_i - \bar{Y}) =& \; \beta_0 + \beta_1 (X_i - \bar{X})+ \varepsilon_i\\
\end{aligned}
$$
接下來，我們要進行兩件事：(1)**改寫模型**與(2)**比較兩模型 OLS 估計式**。首先，根據上面的描述，我們已經知道原本模型與新的模型變數之間的關係，所以我們先寫下兩者之間的轉換。
$$
\begin{aligned}
X_i^* =& X_i - \bar{X} \Rightarrow X_i = X_i^* + \bar{X}\\
Y_i^* =& Y_i - \bar{Y} \Rightarrow Y_i = Y_i^* + \bar{Y}
\end{aligned}
$$
然後進行第一步，注意到這邊我們要以原本模型的變數來表達，否則我們將無法進行計算；而我們的目標是求出
$$
Y = \square + \square X_i + \varepsilon_i
$$
的型態。
$$
(Y_i^* + \bar{Y}) = \beta_0 + \beta_1(X_i^* + \bar{X}) + \varepsilon_i
$$
由於 $\bar{X}$ 與 $\bar{Y}$ 為常數，因此我們可以將其移項，整理之後得到：
$$
Y_i^* = (\beta_0 + \beta_1 \bar{X} - \bar{Y}) + \beta_1 X_i^* + \varepsilon_i
$$
我們就得到了新的模型，且符合我們的要求：$Y = \square + \square X_i + \varepsilon_i$。求得新的模型之型態後，接著我們要分析其與原本模型的 OLS 估計式之間的關係。可以看到前面括號部分長得很像截距項，因此我們以 $\hat{\beta_0^*}$ 代表新模型的 OLS 估計式，其等於
$$
\hat{\beta_0^*} = \hat{\beta_0} + {\hat\beta_1} \bar{X} - \bar{Y}
$$
我們知道 $\hat{\beta_0} = \bar{Y} - \beta_1 \bar{X}$，因此上式可以改寫成
$$
\hat{\beta_0^*} = (\bar{Y} - \beta_1 \bar{X}) + {\hat\beta_1} \bar{X} - \bar{Y} = 0
$$
而斜率項則為 $\hat{\beta_1^*} = \hat{\beta_1} = \frac{\sum^n_{i=1}Y_i(X_i - \bar{X})}{\sum^n_{i=1}(X_i - \bar{X})^2}$。而這個方法也可以不僅可以用於簡單迴歸模型，多元迴歸模型亦可以在同樣的框架之下進行變數的線性轉換。
:::info
:pencil2: 練習
:::
考慮下列兩個迴歸模型：
$$
\begin{aligned}
\text{Model A: }& y_i = \beta_0 + \beta_1 x_i + \varepsilon_i,\\
\text{Model B: }& y_i = \gamma_0 + \gamma_1(x_i - \bar{x}) + \nu_i
\end{aligned}
$$
其中 $\bar{x} = n^{-1}\sum^n_{i=1}x_i$。
- 請求出 $\beta_0$ 與 $\gamma_0$ 的 OLS 估計式，並比較其兩者變異數之大小。
- 請求出 $\beta_1$ 與 $\gamma_1$ 的 OLS 估計式，並比較其兩者變異數之大小。
就 $\text{Model A}$ 而言，我們可以找出 $\varepsilon_i$ 之最小平方和。因此我們可以先定義：
$$
Q(\beta_0, \beta_1) = \sum^n_{i=1} \left(y_i - \beta_0 - \beta_1 x_i \right)^2
$$
接著我們要求其最小值，故在一階條件之下，
$$
\underset{\beta_0, \beta_1}{\arg \min }\; Q(\beta_0, \beta_1)
$$
接著我們可以對兩個變數分別取偏微分，得到
$$
\begin{equation}
\frac{\partial}{\partial\beta_0}Q(\beta_0, \beta_1) = -2 \sum^n_{i=1} \left(y_i - \beta_0 - \beta_1 x_i \right) = 0
\end{equation}
$$
$$
\begin{equation}
\frac{\partial}{\partial\beta_1}Q(\beta_0, \beta_1) = -2 x_i \sum^n_{i=1} \left(y_i - \beta_0 - \beta_1 x_i \right) = 0
\end{equation}
$$
根據上式$(1)$，我們得到
$$
\hat{\beta_0} = \frac{1}{n} \sum^n_{i=1}\left(y_i - \hat{\beta_1}x_i\right) = \bar{y} - \hat{\beta_1}\bar{x}.
$$
而利用 $(2) - \bar{x}(1)$，我們可以求出
$$
\sum^n_{i=1}\left(x_i - \bar{x}\right)\left(y_i - \hat{\beta_0}-\hat{\beta_1}x_i\right) = 0
$$
接著利用 $\hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{x}$ 的結果，
$$
\sum^n_{i=1}\left(x_i - \bar{x}\right)\left[\left(y_i - \bar{y}\right)-\hat{\beta_1}\left(x_i - \bar{x}\right)\right] = 0.
$$
移項之後可得到
$$
\begin{aligned}
&\hat{\beta_1}\sum^n_{i=1}\left(x_i - \bar{x}\right)^2 = \sum^n_{i=1}\left(x_i - \bar{x}\right)\left(y_i - \bar{y}\right)\\
&\Rightarrow \hat{\beta_1} = \frac{\sum^n_{i=1}\left(x_i - \bar{x}\right)\left(y_i - \bar{y}\right)}{\sum^n_{i=1}\left(x_i - \bar{x}\right)^2}
\end{aligned}
$$
而對於 $\text{Model B}$，我們可利用上述變數線性轉換的規則，我們發現新的解釋變數變成 $x_i^* = (x_i - \bar{x})$，$y_i^* = y_i$。接著改寫模型，
$$
\begin{aligned}
y_{i}^* =&\; \beta_0 + \beta_1\left(x_i^* + \bar{x}\right) + \varepsilon_i\\
=&\; \beta_0 + \beta_1x_i^* + \beta_1 \bar{x} + \varepsilon_i \\
=&\; \left(\beta_0+\beta_1 \bar{x}\right) + \beta_1 x_i^* + \varepsilon_i
\end{aligned}
$$
我們可以比較兩模型的參數：
$$
\hat{\gamma_0} = \left(\beta_0+\beta_1 \bar{x}\right) = \bar{y}, \; \hat{\gamma_1} = \beta_1 = \frac{\sum^n_{i=1}\left(x_i - \bar{x}\right)\left(y_i - \bar{y}\right)}{\sum^n_{i=1}\left(x_i - \bar{x}\right)^2}
$$
而我們可以計算 $\hat{\beta_1}$ 的變異數，
$$
\begin{aligned}
\operatorname{var}(\hat{\beta_1}) =&\; \operatorname{var}\left[\frac{\sum^n_{i=1}\left(x_i - \bar{x}\right)\left(y_i - \bar{y}\right)}{\sum^n_{i=1}\left(x_i - \bar{x}\right)^2}\right]\\
=&\; \operatorname{var}\left[\frac{\sum^n_{i=1}y_i\left(x_i - \bar{x}\right)}{\sum^n_{i=1}\left(x_i - \bar{x}\right)^2}\right]\\
=&\; \frac{1}{\left[\sum^n_{i=1}\left(x_i-\bar{x}\right)\right]^2}\operatorname{var}\left[\sum^n_{i=1}y_i\left(x_i - \bar{x}\right)\right]\\
=&\; \frac{1}{\left[\sum^n_{i=1}\left(x_i-\bar{x}\right)\right]^2}\operatorname{var}\left[\sum^n_{i=1}\left(\beta_0+\beta_1x_i+\varepsilon_i\right)\left(x_i-\bar{x}\right)\right]\\
=&\; \frac{1}{\left[\sum^n_{i=1}\left(x_i-\bar{x}\right)\right]^2}\left[\sum^n_{i=1}\left(x_i-\bar{x}\right)\right]^2\operatorname{var}(\varepsilon_i)\\
=&\; \frac{\sigma^2}{\sum^n_{i=1}\left(x_i-\bar{x}\right)^2}
\end{aligned}
$$
而 $\hat{\beta_0}$ 的變異數則為
$$
\begin{aligned}
\operatorname{var}(\hat{\beta_0}) =&\; \operatorname{var}\left(\bar{y}-\hat{\beta_0}\bar{x}\right)\\
=&\; \operatorname{var}\left(\beta_0 + \beta_1 \bar{x} + \bar{\varepsilon_i}-\hat{\beta_1}\bar{x}\right)\\
=&\; \bar{x}^2 \operatorname{var}(\hat{\beta_1}) + \operatorname{var}(\bar{\varepsilon_i})\\
=&\; \frac{\bar{x}^2\sigma^2}{\sum^n_{i=1}\left(x_i-\bar{x}\right)^2} + \operatorname{var}\left(\frac{1}{n}\sum^n_{i=1}\varepsilon_i\right)\\
=&\; \frac{\bar{x}^2\sigma^2}{\sum^n_{i=1}\left(x_i-\bar{x}\right)^2} + \frac{\sigma^2}{n}
\end{aligned}
$$
而就 $\text{Model B}$ 的參數而言，$\hat{\gamma_1} = \hat{\beta_1}$，而 $\hat{\gamma_0}$ 則是
$$
\operatorname{var}(\hat{\gamma_0}) = \operatorname{var}\left(\beta_0 + \hat{\beta_1}\bar{x} + \bar{\varepsilon_i}\right) = \frac{\sigma^2}{n}
$$
而因為 $\hat{\beta_0}$ 較 $\hat{\gamma_0}$ 多了 $\frac{\bar{x}^2\sigma^2}{\sum^n_{i=1}\left(x_i-\bar{x}\right)^2}$ 這一項，且其分子分母皆為平方，因此在假定樣本屬於實數( $x_i \in \mathbb{R}$ )之下，其必大於 $0$，因此 $\hat{\beta_0}$ 必定大於 $\hat{\gamma_0}$。
Sign in with Wallet
Connect another wallet