# 線性迴歸的統計性質與假設檢定
{%hackmd @themes/orangeheart %}
###### tags: `Stats`
回想到古典統計假設二：對於 $i = 1,2,3, \cdots, n$
$$
y_{i}=\beta_{0}+\beta_{1} x_{i 1}+\beta_{2} x_{i 2}+\cdots+\beta_{K} x_{i K}+\varepsilon_{i}
$$
其中 $\left\{x_{i 1}, x_{i 2}, \ldots, x_{i K}\right\}$ 是非隨機的，$\mathbb{E}(\varepsilon_i)= 0$、$\operatorname{var}(\varepsilon_i) = \sigma^2$，且 $\operatorname{cov}(\varepsilon_i, \varepsilon_j)=0$。而根據統計古典假設一：
$$
\mathbb{E}(\hat{\beta})=\beta, \operatorname{cov}(\hat{\beta})=\sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}
$$
且 $\mathbb{E}(s^2) = \sigma^2$，其中
$$
s^2 = \frac{1}{n-K-1}\sum^n_{i=1}\hat{\varepsilon_i}^2= \frac{1}{n-K-1}\hat{\varepsilon}^{\top}\hat{\varepsilon}
$$
而分母為 $n - K -1$ 的統計意涵解釋為：「有 $n$ 個變數，其中 $K$ 為解釋變數個數，作為懲罰項；$1$ 則為截距項個數」。不過，我們關心的是 $\hat{\beta}$ 的分配是什麼？我們可以做一個最基本的假設，就是令其服從常態分配，即
$$
\varepsilon_{i} \stackrel{i . i . d .}{\sim} \mathcal{N}\left(0, \sigma^{2}\right) .
$$
也就是說，
$$
\mathbf{y} \sim \mathcal{N}\left(\mathbf{X} \beta, \sigma^{2} \mathbf{I}\right)
$$
且
$$
\hat{\beta}=\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{X}^{\top} \mathbf{y} \sim \mathcal{N}\left(\beta, \sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right)
$$
## 簡單線性迴歸的假設檢定
對於經濟學家來說，當我們想要了解一些經濟現象時，便會建構一些模型，試圖去驗證、解釋它。假設一個 Cobb-Douglas 函數
$$
Q = AK^{\alpha}L^{\beta}
$$
對其取對數並加上誤差項，便可以得到實證模型：
$$
q_{t}=\mu+\alpha k_{t}+\beta \ell_{t}+\epsilon_{t}, t=1,2, \ldots, T
$$
其中 $q_t = \ln Q_t$、$k_t = \ln K_t$，且 $\ell_t = \ln L_t$。我們可以檢測其是否為[固定規模報酬(constant return to scale, CRTS)](https://www.economics.utoronto.ca/osborne/2x3/tutorial/RTSFRM.HTM)，即檢定：
$$
\begin{aligned}
H_0: \alpha + \beta = 1\\
H_1: \alpha + \beta < 1
\end{aligned}
$$
另一個有趣經濟學（或財務金融上）實證問題：[效率市場假說](https://www.stockfeel.com.tw/%E9%97%9C%E6%96%BC%E6%95%88%E7%8E%87%E5%B8%82%E5%A0%B4%E5%81%87%E8%AA%AA/)。
$$
r_{t}=\mu+l_{t-1}^{\top} \gamma+\epsilon_{t}
$$
其中 $r_t$ 代表在 $t-1$ 至 $t$ 期間資產所得之報酬，而 $l_{t-1}$ 則是在時點 $t-1$ 時的公開資訊。如果我們要驗證它，就必須檢測 $\gamma$ 是否為$0$。
### Pizza 的例子
```r
options(width = 70)
Pizza <- read.csv("pizza_delivery.csv")
head(Pizza)
## day date time operator branch driver temperature
## 1 Thursday 01-May-14 35.12837 Laura East Bruno 68.28772
## 2 Thursday 01-May-14 25.20307 Melissa East Salvatore 70.99779
## 3 Thursday 01-May-14 45.64340 Melissa West Salvatore 53.39415
## 4 Thursday 01-May-14 29.37430 Melissa East Salvatore 70.30660
## 5 Thursday 01-May-14 29.99461 Melissa West Salvatore 71.50169
## 6 Thursday 01-May-14 40.25432 Melissa Centre Bruno 60.75950
## bill pizzas free_wine got_wine discount_customer
## 1 58.4 4 0 0 1
## 2 26.4 2 0 0 0
## 3 58.1 3 1 0 0
## 4 35.2 3 0 0 0
## 5 38.4 2 0 0 0
## 6 61.8 4 1 1 0
```
經理關心的事情是：外送員是否會「準時」的送達客戶家中，且沒有送出免費的酒。我們可以做一個時間($\texttt{time}$)與外送員($\texttt{operator}$)
```r
lm(time~operator, data = Pizza)
##
## Call:
## lm(formula = time ~ operator, data = Pizza)
##
## Coefficients:
## (Intercept) operatorMelissa
## 34.1030 0.2551
```
我們可以令一個**二元變數(dummy variable)**，$0$ 代表外送員 Laura，$1$ 代表 Melissa。那為什麼我們不單純就建立兩個分別給這兩位外送員的解釋變數呢？當然，很明顯這會造成**線性重合**的問題，[^1] 所以我們要盡可能避免它。回到上面的例子，$\texttt{operatorMelissa}$ 係數檢定出來的結果是 $0.2551$。
箱型圖可以告訴我們一件事： Melissa可能只是運氣不好！多花了大約15秒。
> Laura: "I am a better operator than Melissa, since the time is shorter when the order is operated by me."
> Melissa: "The time difference is only 0.2551 minutes (15.31 seconds)! I don’t think 0.2551 is very different from zero $\cdots$""
那我們要怎麼用假設檢定去驗證 Laura 與 Melissa 他們各自的論述呢？沒錯！我們可以進行假設檢定！令 $\beta_{om}$ 代表 Melissa，且
$$
\begin{aligned}
H_0: \beta_{om} = 0\\
H_1: \beta_{om} > 0
\end{aligned}
$$
我們已經知道 $\hat{\beta}$ 矩陣服從常態分配：
$$
\hat{\beta} \sim \mathcal{N}\left(\beta, \sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right)
$$
```r
Pizza <- read.csv("pizza_delivery.csv")
pizza.lm <- lm(time~operator,data=Pizza,x=TRUE)
sum(pizza.lm$residualsˆ2)/pizza.lm$df.residual
## [1] 41.76744
```
接著我們可以利用 `solve()` 這個指令解線性關係。[^2]
```r
solve(t(pizza.lm$x) %*% pizza.lm$x)
## (Intercept) operatorMelissa
## (Intercept) 0.001567398 -0.001567398
## operatorMelissa -0.001567398 0.003159755
```
在統計古典假設二之下，
$$
\hat{\beta}_{o m} \sim \mathcal{N}\left(\beta, \sigma^{2} \times 0.003159755\right)
$$
而根據我們上面設定的虛無假設：$H_0:\beta_{om} = 0$，
$$
\hat{\beta}_{\mathrm{om}} \sim \mathcal{N}\left(0, \sigma^{2} \times 0.003159755\right)
$$
事實是，我們不知道 $\sigma^2$，但我們可以知道其不偏估計式 $s^2$，
$$
s^2 = 41.76744
$$
則 $41.76744 \times 0.003159755 = 0.1319749$，這就代表我們可以說 $\beta_{om} \sim \mathcal{N}(0,0.1319749)$ 了嗎？考慮一個簡單線性迴歸的虛無假設：$H_0 = \mathbf{c}^{\top}\beta = \gamma$，而因為
$$
\hat{\beta} \sim \mathcal{N}\left(\beta, \sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right)
$$
所以我們得到
$$
\mathbf{c}^{\top}\hat{\beta} \sim \mathcal{N}\left(\mathbf{c}^{\top}\beta, \sigma^2 \mathbf{c}^{\top}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\mathbf{c}\right)
$$
舉例來說，如果我們令 $\mathbf{c}^{\top} = [10, \cdots, 0]$，以及 $\gamma = 0$，那麼虛無假設則為$H_0:\beta_0 = 0$。而在虛無假設下
$$
\mathbf{c}^{\top}\hat{\beta} \sim \mathcal{N}\left(\gamma, \sigma^2 \mathbf{c}^{\top}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\mathbf{c}\right)
$$
$Z$ 值等於
$$
Z = \frac{\mathbf{c}^{\top}\hat{\beta} - \gamma }{\sigma\sqrt{ \mathbf{c}^{\top}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\mathbf{c}}} \sim \mathcal{N}(0, 1)
$$
但實際情況 $\sigma^2$ 通常是未知的，因此我們必須透過樣本變異數$s^2$進行估計。故得到 $T$ 值為
$$
T = \frac{\mathbf{c}^{\top}\hat{\beta} - \gamma }{s\sqrt{ \mathbf{c}^{\top}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\mathbf{c}}} \sim t(n-K-1)
$$
### 單個樣本的 $T$ 檢定
我們可以用 Monte Carlo 模擬實際的情況。假設一個簡單迴歸式為 $\mathbf{y} = \mathbf{X}\beta + \varepsilon$，其中
$$
\mathbf{X}=\left[\begin{array}{ll}
1 & 1 \\
1 & 2 \\
1 & 3
\end{array}\right], \beta=\left[\begin{array}{l}
2 \\
1
\end{array}\right], \text { 且 } \varepsilon \sim \mathcal{N}(0, 0.25 \times \mathbf{I})
$$
```r
# Generate 5000 zeros
Z = T = numeric(5000)
# Create a matrix of X
X <- cbind(rep(1,3), (1:3))
# Solve (XᵀX)⁻¹
Xinv <- solve(t(X) %*% X)
# Create a matrix of Y
for (r in 1:5000) {
Y <- X %*% c(2, 1) + 0.5 * rnorm(3)
beta.h <- Xinv %*% t(X) %*% Y
s <- sqrt(sum((Y - X %*% beta.h)^2))
Z[r] <- (beta.h[2] - 1) / (0.5 * sqrt(Xinv[2,2]))
T[r] <- (beta.h[2] - 1) / (s * sqrt(Xinv[2,2]))
}
# Print out the results
c(mean(Z<qnorm(0.025)),mean(Z<qnorm(0.05)),mean(Z<qnorm(0.1)))
## [1] 0.0268 0.0532 0.1018
c(mean(T<qnorm(0.025)),mean(T<qnorm(0.05)),mean(T<qnorm(0.1)))
## [1] 0.1478 0.1746 0.2124
c(mean(T<qt(0.025,1)),mean(T<qt(0.05,1)),mean(T<qt(0.1,1)))
## [1] 0.0254 0.0498 0.1000
```
```r
summary(lm(time~operator,data=Pizza))
##
## Call:
## lm(formula = time ~ operator, data = Pizza)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.0921 -4.1931 0.2245 4.2840 18.9933
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.1030 0.2559 133.286 <2e-16 ***
## operatorMelissa 0.2551 0.3633 0.702 0.483
## ---
## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
##
## Residual standard error: 6.463 on 1264 degrees of freedom
## Multiple R-squared: 0.00039, Adjusted R-squared: -0.0004009
## F-statistic: 0.4931 on 1 and 1264 DF, p-value: 0.4827
```
我們算出跟上面一樣的結果： $\texttt{opeartorMelissa}$ 的係數為$0.2551$，標準差為 $0.3663$，$T$ 值為 $0.702$，而在給定顯著水準為 $\alpha = 0.05$ 之下，$0.702< 1.64606$，代表我們沒有足夠證據拒絕虛無假設。那麼 $p$ 值呢？$p$ 值代表 $\mathbb{P}(T>t)$，
```r
1 - pt(0.702, 1264)
## [1] 0.2414042
```
因為 $0.2414042 > 0.05$，也告訴我們沒有足夠證據拒絕虛無假設。對於對立假設，我們可以計算其拒絕域
```r
c(qt(0.025, 1264), qt(0.975, 1264))
## [1] -1.961843 1.961843
```
因為 $0.702 \nless −1.961843$ 且 $0.702 \ngtr 1.961843$：不拒絕虛無假設；最後 $p$ 值則是 $\mathbb{P}(T > |t|)$，
```r
pt(-0.702, 1264) + (1 - pt(0.702, 1264))
## [1] 0.4828084
```
因為 $0.4828084 > 0.05$，我們不拒絕虛無假設。
### 兩個樣本的 $T$ 檢定
給定 $\{X_i\}^{n_X}_{i=1} \overset{i.i.d}{\sim} \mathcal{N}(\mu_X, \sigma^2)$ 以及 $\{Y_i\}^{n_Y}_{i=1} \overset{i.i.d}{\sim} \mathcal{N}(\mu_Y, \sigma^2)$。我們設定虛無假設為：$H_0:\mu_X - \mu_Y = 0$，母體變異數可被$s^2$估計為
$$
\begin{aligned}
s^{2} &=\frac{\left(n_{X}-1\right) s_{X}^{2}+\left(n_{Y}-1\right) s_{Y}^{2}}{n_{X}+n_{Y}-2} \\
&=\frac{\sum_{i=1}^{n_{X}}\left(X_{i}-\bar{X}\right)^{2}+\sum_{i=1}^{n_{Y}}\left(Y_{i}-\bar{Y}\right)^{2}}{n_{X}+n_{Y}-2}
\end{aligned}
$$
因此在虛無假設下：
$$
T=\frac{\bar{X}-\bar{Y}}{\sqrt{\frac{s^{2}}{n_{X}}+\frac{s^{2}}{n_{Y}}}} \sim t\left(n_{X}+n_{Y}-2\right)
$$
```r
Pizza <- read.csv("pizza_delivery.csv")
summary(lm(time~branch+bill+operator,data=Pizza))
##
## Call:
## lm(formula = time ~ branch + bill + operator, data = Pizza)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.8794 -3.9774 -0.3761 3.7712 17.6560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.19138 0.78752 33.258 < 2e-16 ***
## branchEast -3.03606 0.42330 -7.172 1.25e-12 ***
## branchWest -0.50339 0.38907 -1.294 0.196
## bill 0.21319 0.01535 13.885 < 2e-16 ***
## operatorMelissa 0.15973 0.31784 0.503 0.615
## ---
## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
##
## Residual standard error: 5.653 on 1261 degrees of freedom
## Multiple R-squared: 0.2369, Adjusted R-squared: 0.2345
## F-statistic: 97.87 on 4 and 1261 DF, p-value: < 2.2e-16
```
令 $\beta_\text{bill}$ 代表帳單的係數，我們欲檢測虛無假設 $H_0 : \beta_\text{bill} = 0$ 是否成立。根據上面的結果，在顯著水準 $\alpha = 0.05$ 下，我們可以拒絕虛無假設。注意到 $\verb|R|$ 在進行迴歸分析的假設檢定時，會透過星號(asterisk)告訴我們顯著水準：
- `.`：代表顯著水準為 $0.1$
- `*`：代表顯著水準為 $0.05$
- `**`：代表顯著水準為 $0.01$
- `***`：代表顯著水準為 $0.001$
## 多元線性迴歸的假設檢定
有時候我們想要同時檢定多個係數。令 $\beta_\text{be}$ 為東區分店的係數，$\beta_\text{bw}$ 為西區分店的係數，我們想要檢定這兩個係數是否為 $0$，即 $H_0 : \beta_\text{be} = \beta_\text{bw} = 0$，也就是檢定由哪家分店處理訂單，與披薩從接線員到送達客戶手中那一刻前所花的時間，是否存在一定的關係。最簡單的想法是：對於這兩個係數，我們個別寫下其 $t$ 值，然後按次序地(sequentially)進行 $T$ 檢定。也就是：
$$
\begin{aligned}
H_{0,1} : \beta_\text{be} = 0\\
H_{0,2} : \beta_\text{bw} = 0
\end{aligned}
$$
而當我們拒絕其中一個虛無假設 $H_{0,1}$ 或 $H_{0,2}$ 後，就說我們可以拒絕 $H_0$。但這個想法跟作法其實不太好，因為這牽涉到假設檢定的基本邏輯。
### 假設檢定的本質：避免扭曲 size (Size Distortion)
我們是先選定一個夠小的機率，當然目前都是承襲 Ronald Fisher 於 100年前的 $\alpha$ 值，即 $0.05$。在虛無假設是正確的前提下，我們寫下一個拒絕域，如果虛無假設是正確的，我們有 $0.05$ 的機率拒絕虛無假設，也就是犯下所謂的**型一錯誤(Type I error)**。但如果我們將虛無假設拆分成兩個（也就是上面的想法）並依序檢定，檢定第一個虛無假設而犯下型一錯誤的機率為 $0.05$；檢定第二個虛無假設而犯下型一錯誤的機率亦為 $0.05$。但是！**連續兩次不犯下型一錯誤的機率**絕對不會是 $0.05$，也就是
$$
\mathbb{P}\left(H_{0,1} \text{ is rejected}\right) + \mathbb{P}\left(H_{0,2} \text{ is rejected} \mid H_{0,1} \text{ is not rejected} \right) > 0.05.
$$
除非 $\mathbb{P}\left(H_{0,2} \text{ is rejected} \mid H_{0,1} \text{ is not rejected} \right) = 0$。一切的麻煩都源自於**我們寫下了兩個統計量**！
### $F$ 檢定
我們可以把虛無假設寫成 $H_0 : \mathbf{R}\beta = \mathbf{r}$，其中 $\mathbf{R}$ 代表一個 $q \times (K+1)$ 的矩陣，$K+1$ 表示我們有幾個解釋變數（包含截距項），$q$ 則是有幾個虛無假設要檢定；而 $\mathbf{r}$ 則是一個 $q \times 1$的向量。回到上面的例子，令
$$
\mathbf{R}=\left[\begin{array}{lllll}
0 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0
\end{array}\right], \;
\beta=\left[\begin{array}{c}
\beta_{0} \\
\beta_{b e} \\
\beta_{b w} \\
\beta_{b i l l} \\
\beta_{o m}
\end{array}\right], \;
\mathbf{r}=\left[\begin{array}{c}
0 \\
0
\end{array}\right]
$$
因此 虛無假設 $H_0 : \beta_\text{be} = \beta_\text{bw} = 0$ 就可以寫成 $H_0 : \mathbf{R}\beta = \mathbf{r}$。在統計古典的假設二之下，
$$
\hat{\beta} = \left(\mathbf{X}^{\top}\mathbf{X}\right)^{-1}\mathbf{X}^{\top}\mathbf{y} \sim \mathcal{N}\left(\beta, \sigma^{2}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1}\right)
$$
因此
$$
\mathbf{R} \hat{\beta}-\mathbf{r} \sim \mathcal{N}\left(\mathbf{R} \beta-\mathbf{r}, \sigma^{2} \mathbf{R}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{R}^{\top}\right)
$$
根據我們的虛無假設，
$$
\mathbf{R} \hat{\beta}-\mathbf{r} \sim \mathcal{N}\left(\mathbf{0}, \sigma^{2} \mathbf{R}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{R}^{\top}\right)
$$
那麼
$$
\frac{1}{\sigma}\left(\mathbf{R}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{R}^{\top}\right)^{-1 / 2}(\mathbf{R} \hat{\beta}-\mathbf{r}) \sim \mathcal{N}(\mathbf{0}, \mathbf{I})
$$
但一樣的問題 $\sigma$ 仍是未知的，而且這是一個 $q \times 1$ 的向量，而不是一個統計量。不過有個方法可以解決上述的問題：因為上面的形式有根號，我們何不把它透過平方消除然後相加總呢？這個想法就是我們之前學過 $\chi^2$ 分配的由來：當 $X_{1}, X_{2}, \ldots, X_{k} \stackrel{\text { i.i.d. }}{\sim} \mathcal{N}(0,1)$，則 $Y=X_{1}^{2}+X_{2}^{2}+\cdots+X_{k}^{2} \sim \chi^{2}(k)$。令
$$
X^{2}=\frac{1}{\sigma^{2}}(\mathbf{R} \hat{\beta}-\mathbf{r})^{\top}\left(\mathbf{R}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{R}^{\top}\right)^{-1}(\mathbf{R} \hat{\beta}-\mathbf{r})
$$
則，
$$
X^{2} \sim \chi^{2}(q)
$$
到了這一步看起來算是可以了，但 $\sigma^2$ 還是未知的，不過很簡單就可以解決：利用樣本變異數的不偏性取代母體變異數，
$$
F=\frac{(\mathbf{R} \hat{\beta}-\mathbf{r})^{\top}\left(\mathbf{R}\left(\mathbf{X}^{\top} \mathbf{X}\right)^{-1} \mathbf{R}^{\top}\right)^{-1}(\mathbf{R} \hat{\beta}-\mathbf{r}) / q}{\hat{\varepsilon}^{\top} \hat{\varepsilon} /(n-K-1)}
$$
其中，$\hat{\varepsilon}^{\top} \hat{\varepsilon} /(n-K-1)$ 就是 $s^2$。則，
$$
F \sim F(q, n-K-1)
$$
考慮一個 Monte Carlo 模擬實驗。令 $\mathbf{y}=\mathbf{X} \beta+\varepsilon$，其中
$$
\mathbf{X}=\left[\begin{array}{ll}
1 & 1 \\
1 & 2 \\
1 & 3
\end{array}\right],\; \beta=\left[\begin{array}{l}
2 \\
1
\end{array}\right], \; \varepsilon \sim \mathcal{N}(\mathbf{0}, 0.25 \times \mathbf{I})
$$
```r
# Generate 5000 zeros
C = F = numeric(5000)
# Create a matrix of X
X <- cbind(rep(1,3), (1:3))
# Solve (XᵀX)⁻¹
Xinv <- solve(t(X) %*% X)
# Create a matrix of Y
for (r in 1:5000) {
Y <- X %*% c(2, 1) + 0.5 * rnorm(3)
beta.h <- Xinv %*% t(X) %*% Y
s2 <- sum((Y - X %*% beta.h)ˆ2)
C[r] <- t(beta.h-c(2,1)) %*% t(X) %*% X %*% (beta.h-c(2,1)) / 0.25
F[r] <- t(beta.h-c(2,1)) %*% t(X) %*% X %*% (beta.h-c(2,1)) / (2*s2)
}
# Print out the results
c(mean(C<qchisq(0.05,2)),mean(C<qchisq(0.1,2)),mean(C<qchisq(0.5,2)))
## [1] 0.0552 0.1006 0.4976
c(mean(F<qf(0.05,2,1)),mean(F<qf(0.1,2,1)),mean(F<qf(0.5,2,1)))
## [1] 0.0532 0.1022 0.5020
```
```r
Pizza <- read.csv("pizza_delivery.csv")
summary(lm(time~branch+bill+operator,data=Pizza))
##
## Call:
## lm(formula = time ~ branch + bill + operator, data = Pizza)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.8794 -3.9774 -0.3761 3.7712 17.6560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.19138 0.78752 33.258 < 2e-16 ***
## branchEast -3.03606 0.42330 -7.172 1.25e-12 ***
## branchWest -0.50339 0.38907 -1.294 0.196
## bill 0.21319 0.01535 13.885 < 2e-16 ***
## operatorMelissa 0.15973 0.31784 0.503 0.615
## ---
## Signif. codes:
## 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
##
## Residual standard error: 5.653 on 1261 degrees of freedom
## Multiple R-squared: 0.2369, Adjusted R-squared: 0.2345
## F-statistic: 97.87 on 4 and 1261 DF, p-value: < 2.2e-16
```
回到原本的例子，虛無假設為 $H_0 : \beta_\text{be} = \beta_\text{bw} = 0$，也就是 $H_0 : \mathbf{R} \beta = \mathbf{r}$，則
```r
Pizza <- read.csv("pizza_delivery.csv")
pizza.lm <- lm(time~branch+bill+operator,data=Pizza)
b <- pizza.lm$coefficients
V <- vcov(pizza.lm)
R <- matrix(c(0,0,1,0,0,1,0,0,0,0),nrow=2)
F <- t(R %*% b) %*% solve(R %*% V %*% t(R)) %*% R %*% b / 2
F
## [1,] 29.32826
1 - pf(F,2,1261)
## [1,] 3.550493e-13
```
如果說我們令虛無假設為 $H_{0}: \beta_\text{be}=\beta_\text{bw}=\beta_\text{bill}=\beta_\text{om}=0$
$H_0 : \mathbf{R} \beta = \mathbf{r}$，其中，
$$
\mathbf{R}=\left[\begin{array}{lllll}
0 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 1
\end{array}\right], \;
\beta=\left[\begin{array}{c}
\beta_0\\
\beta_{b e} \\
\beta_{b w} \\
\beta_{b i l l} \\
\beta_{o m}
\end{array}\right], \;
\mathbf{r}=\left[\begin{array}{l}
0 \\
0 \\
0 \\
0
\end{array}\right]
$$
```r
Pizza <- read.csv("pizza_delivery.csv")
pizza.lm <- lm(time~branch+bill+operator,data=Pizza)
b <- pizza.lm$coefficients
V <- vcov(pizza.lm)
F <- t(b[2:5]) %*% solve(V[2:5,2:5]) %*% b[2:5] / 4
F
## [1,] 97.87125
1 - pf(F,4,1261)
## [1,] 0
```
代表我們可以拒絕虛無假設，其實際值小於 $\verb|R|$ 的寬容值 $2.2 \times 10^{-16}$。一個簡單的計算方式是：
$$
F=\frac{\sum_{i=1}^{n}\left(\hat{y}_{i}-\bar{y}\right)^{2} / K}{\sum_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2} /(n-K-1)} = \frac{SSE}{SSR/(n-K-1)}
$$
可以看出 $F$ 檢定統計量與 $R^2$ 有一定程度的（非線性）關係：$R^2$ 越大，$F$ 檢定統計量也越大。而我們可以建構一個多元迴歸模型的 ANOVA（變異數分析）表進行計算分類後的資料是否能夠良好的詮釋模型。
| Source | SS | df | MS | $\phi$ |
| -------- | -------- | -------- | -------- | -------- |
| Factor | SSR | $k$ | MSR = $\frac{SSR}{k}$ | $\phi_F = \frac{MSR}{MSE}$|
| Error | SSE | $n-k-1$ | MSE = $\frac{SSE}{n-k-1}$|
| Total | SST | $n-1$ |
其中檢定統計量 $\phi_F = \frac{MSR}{MSE} \sim \operatorname{F}(1, n-2)$。而在 $\verb|R|$ 中，我們可以用 `anova()` 進行變異數分析。
```r
Pizza <- read.csv("pizza_delivery.csv")
anova(lm(time~branch,data=Pizza))
## Analysis of Variance Table
##
## Response: time
## Df Sum Sq Mean Sq F value Pr(>F)
## branch 2 6334 3166.8 86.05 < 2.2e-16 ***
## Residuals 1263 46481 36.8
## ---
## Signif. codes:
## 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
```
## 檢證誤差項服從常態分配的假設
我們前面都假設誤差項服從常態分配：如果不是服從常態分配的話，前面做的任何分析都是徒勞無功。而我們在看不到誤差項的情況之下，又要如何檢測其是否服從常態分配呢？由於我們已經知道 $\hat{\beta}$ 是 $\beta$ 的不偏估計式，而只要殘差與誤差項夠接近，就可以印證上述的假設為真。首先我們先將誤差項進行標準化，得到
$$
\varepsilon_i \overset{i.i.d}{\sim} \mathcal{N}\left(0, \sigma^2\right)
$$
並由小排到大。如果我們有 $n$ 個服從常態分配的殘差與 $n$ 個經過標準化後服從常態分配的誤差項，分別將其畫在 $x$ 軸與 $y$ 軸，此兩者應該要分布在 $45^{\circ}$ 線上，形成 QQPlot。[^3]舉例來說，
```r
X <- runif(1000)
Y <- 1 + X + 0.2 * rnorm(1000)
lm <- lm(Y~X)
plot(pizza.lm, which = 2)
```
如果誤差項不是服從常態分配，而是服從均勻分配的話，QQPlot 則會長得像下面那樣，並非在 $45^{\circ}$ 線上。
```r
X <- runif(1000)
Y <- 1 + X + 0.2 * runif(1000) - 0.5
lm <- lm(Y~X)
plot(pizza.lm, which = 2)
```
最後，我們可以畫出 `pizza` 中誤差項的 QQPlot。
```r
pizza.lm <- lm(time~branch + bill + operator, data = pizza)
plot(pizza.lm, which = 2)
```
但是，如果誤差項真的不是服從常態分配，我們還能夠繼續進行上述的分析嗎？答案是可以的！事實上，利用**漸進理論**──只要樣本數夠大，我們還是能夠維持原本的假設。
[^1]: 當2個（或以上）的自變數互不獨立(即彼此相關)，就是具有「線性重合」，將會使迴歸模型中存在著重複的自變數，提高某一自變數的解釋力與預測力，使得理論的建構不正確。
[^2]: 如果我們要解出 $\mathbf{Ax = b}$ 的話，我們可以輸入 `solve(A, b)`，便會算出 $\mathbf{x = A^{-1}b}$；如果我們僅輸入 `solve(A)`，$\texttt{R}$便會自動認定我們要計算 $\mathbf{Ax = I}$，得到 $\mathbf{x = A^{-1}I = A^{-1}}$。當然，看到上面的計算過程，對於線性代數理解夠的人就可以知道， `solve()` 可以拿來取反矩陣(inverse matrix)。
[^3]:`plot()` 中的指令 `which = 2` 代表要畫出 QQPlot。
Sign in with Wallet
Connect another wallet