線性迴歸的統計性質與假設檢定

tags: Stats

回想到古典統計假設二：對於

$i=1,2,3,\cdots ,n$

$$y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_K{x}_{iK}+{\epsilon }_i$$
其中

$\{x_{i1},x_{i2},\dots ,x_{iK}\}$ 是非隨機的，

$\mathbb{E}({\epsilon }_i)=0$、

$\mathrm{var}({\epsilon }_i)={\sigma }^2$，且

$\mathrm{cov}({\epsilon }_i,{\epsilon }_j)=0$。而根據統計古典假設一：

$$\mathbb{E}(\hat{\beta })=\beta ,\mathrm{cov}(\hat{\beta })={\sigma }^2{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}$$
且

$\mathbb{E}(s^2)={\sigma }^2$，其中

$$s^2=\frac{1}{n-K-1}\sum _{i=1}^n{\widehat{{\epsilon }_i}}^2=\frac{1}{n-K-1}{\hat{\epsilon }}^{\mathrm{\top }}\hat{\epsilon }$$

而分母為

$n-K-1$ 的統計意涵解釋為：「有

$n$ 個變數，其中

$K$ 為解釋變數個數，作為懲罰項；

$1$ 則為截距項個數」。不過，我們關心的是

$\hat{\beta }$ 的分配是什麼？我們可以做一個最基本的假設，就是令其服從常態分配，即

$${\epsilon }_i\stackrel{i.i.d.}{\sim }\mathcal{N}(0,{\sigma }^2).$$
也就是說，

$$\mathbf{y}\sim \mathcal{N}(\mathbf{X}\beta ,{\sigma }^2\mathbf{I})$$
且

$$\hat{\beta }={\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}\sim \mathcal{N}(\beta ,{\sigma }^2{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1})$$

簡單線性迴歸的假設檢定

對於經濟學家來說，當我們想要了解一些經濟現象時，便會建構一些模型，試圖去驗證、解釋它。假設一個 Cobb-Douglas 函數

$$Q=A{K}^{\alpha }{L}^{\beta }$$
對其取對數並加上誤差項，便可以得到實證模型：

$$q_t=\mu +\alpha k_t+\beta {\ell }_t+{ϵ}_t,t=1,2,\dots ,T$$
其中

$q_t=\ln Q_t$、

$k_t=\ln K_t$，且

${\ell }_t=\ln L_t$。我們可以檢測其是否為固定規模報酬(constant return to scale, CRTS)，即檢定：

$$\begin{array}{r}{H}_0:\alpha +\beta =1\\ H_1:\alpha +\beta <1\end{array}$$
另一個有趣經濟學（或財務金融上）實證問題：效率市場假說。

$$r_t=\mu +l_{t-1}^{\mathrm{\top }}\gamma +{ϵ}_t$$
其中

$r_t$ 代表在

$t-1$ 至

$t$ 期間資產所得之報酬，而

$l_{t-1}$ 則是在時點

$t-1$ 時的公開資訊。如果我們要驗證它，就必須檢測

$\gamma$ 是否為

$0$。

Pizza 的例子

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

經理關心的事情是：外送員是否會「準時」的送達客戶家中，且沒有送出免費的酒。我們可以做一個時間(

$\mathtt{\text{time}}$)與外送員(

$\mathtt{\text{operator}}$)

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] 所以我們要盡可能避免它。回到上面的例子，

$\mathtt{\text{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{array}{r}{H}_0:{\beta }_{om}=0\\ H_1:{\beta }_{om}>0\end{array}$$
我們已經知道

$\hat{\beta }$ 矩陣服從常態分配：

$$\hat{\beta }\sim \mathcal{N}(\beta ,{\sigma }^2{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1})$$

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]

solve(t(pizza.lm$x) %*% pizza.lm$x)
##                  (Intercept)   operatorMelissa
## (Intercept)      0.001567398      -0.001567398
## operatorMelissa -0.001567398       0.003159755

在統計古典假設二之下，

$${\hat{\beta }}_{om}\sim \mathcal{N}(\beta ,{\sigma }^2\times 0.003159755)$$
而根據我們上面設定的虛無假設：

$H_0:{\beta }_{om}=0$，

$${\hat{\beta }}_{\mathrm{om}}\sim \mathcal{N}(0,{\sigma }^2\times 0.003159755)$$

事實是，我們不知道

${\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}}^{\mathrm{\top }}\beta =\gamma$，而因為

$$\hat{\beta }\sim \mathcal{N}(\beta ,{\sigma }^2{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1})$$

所以我們得到

$${\mathbf{c}}^{\mathrm{\top }}\hat{\beta }\sim \mathcal{N}({\mathbf{c}}^{\mathrm{\top }}\beta ,{\sigma }^2{\mathbf{c}}^{\mathrm{\top }}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}\mathbf{c})$$

舉例來說，如果我們令

${\mathbf{c}}^{\mathrm{\top }}=[10,\cdots ,0]$，以及

$\gamma =0$，那麼虛無假設則為

$H_0:{\beta }_0=0$。而在虛無假設下

$${\mathbf{c}}^{\mathrm{\top }}\hat{\beta }\sim \mathcal{N}(\gamma ,{\sigma }^2{\mathbf{c}}^{\mathrm{\top }}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}\mathbf{c})$$

$Z$ 值等於

$$Z=\frac{{\mathbf{c}}^{\mathrm{\top }}\hat{\beta }-\gamma }{\sigma \sqrt{{\mathbf{c}}^{\mathrm{\top }}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}\mathbf{c}}}\sim \mathcal{N}(0,1)$$

但實際情況

${\sigma }^2$ 通常是未知的，因此我們必須透過樣本變異數

$s^2$進行估計。故得到

$T$ 值為

$$T=\frac{{\mathbf{c}}^{\mathrm{\top }}\hat{\beta }-\gamma }{s\sqrt{{\mathbf{c}}^{\mathrm{\top }}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}\mathbf{c}}}\sim t(n-K-1)$$

單個樣本的

$T$ 檢定

我們可以用 Monte Carlo 模擬實際的情況。假設一個簡單迴歸式為

$\mathbf{y}=\mathbf{X}\beta +\epsilon$，其中

$$\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{ 且 }\epsilon \sim \mathcal{N}(0,0.25\times \mathbf{I})$$

# 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

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

我們算出跟上面一樣的結果：

$\mathtt{\text{opeartorMelissa}}$ 的係數為

$0.2551$，標準差為

$0.3663$，

$T$ 值為

$0.702$，而在給定顯著水準為

$\alpha =0.05$ 之下，

$0.702<1.64606$，代表我們沒有足夠證據拒絕虛無假設。那麼

$p$ 值呢？

$p$ 值代表

$\mathbb{P}(T>t)$，

1 - pt(0.702, 1264)

## [1] 0.2414042

因為

$0.2414042>0.05$，也告訴我們沒有足夠證據拒絕虛無假設。對於對立假設，我們可以計算其拒絕域

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|)$，

pt(-0.702, 1264) + (1 - pt(0.702, 1264))

## [1] 0.4828084

因為

$0.4828084>0.05$，我們不拒絕虛無假設。

兩個樣本的

$T$ 檢定

給定

$\{X_i{\}}_{i=1}^{n_X}\stackrel{i.i.d}{\sim }\mathcal{N}({\mu }_X,{\sigma }^2)$ 以及

$\{Y_i{\}}_{i=1}^{n_Y}\stackrel{i.i.d}{\sim }\mathcal{N}({\mu }_Y,{\sigma }^2)$。我們設定虛無假設為：

$H_0:{\mu }_X-{\mu }_Y=0$，母體變異數可被

$s^2$估計為

$$\begin{array}{rl}{s}^2& =\frac{(n_X-1)s_X^2+(n_Y-1)s_Y^2}{n_X+n_Y-2}\\ & =\frac{\sum _{i=1}^{n_X}{(X_i-\overline{X})}^2+\sum _{i=1}^{n_Y}{(Y_i-\overline{Y})}^2}{n_X+n_Y-2}\end{array}$$
因此在虛無假設下：

$$T=\frac{\overline{X}-\overline{Y}}{\sqrt{\frac{s^2}{n_X}+\frac{s^2}{n_Y}}}\sim t(n_X+n_Y-2)$$

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$ 下，我們可以拒絕虛無假設。注意到

$\text{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{array}{r}{H}_{0,1}:{\beta }_{\text{be}}=0\\ H_{0,2}:{\beta }_{\text{bw}}=0\end{array}$$

而當我們拒絕其中一個虛無假設

$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}(H_{0,2}\text{ is rejected}\mid H_{0,1}\text{ is not rejected})>0.05.$$

除非

$\mathbb{P}(H_{0,2}\text{ is rejected}\mid H_{0,1}\text{ is not rejected})=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 }_{be}\\ {\beta }_{bw}\\ {\beta }_{bill}\\ {\beta }_{om}\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}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}\sim \mathcal{N}(\beta ,{\sigma }^2{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1})$$

因此

$$\mathbf{R}\hat{\beta }-\mathbf{r}\sim \mathcal{N}(\mathbf{R}\beta -\mathbf{r},{\sigma }^2\mathbf{R}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{R}}^{\mathrm{\top }})$$

根據我們的虛無假設，

$$\mathbf{R}\hat{\beta }-\mathbf{r}\sim \mathcal{N}(\mathbf{0},{\sigma }^2\mathbf{R}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{R}}^{\mathrm{\top }})$$

那麼

$$\frac{1}{\sigma }{\left(\mathbf{R}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{R}}^{\mathrm{\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,\dots ,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}{)}^{\mathrm{\top }}{\left(\mathbf{R}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{R}}^{\mathrm{\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}{)}^{\mathrm{\top }}{\left(\mathbf{R}{\left({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\right)}^{-1}{\mathbf{R}}^{\mathrm{\top }}\right)}^{-1}(\mathbf{R}\hat{\beta }-\mathbf{r})/q}{{\hat{\epsilon }}^{\mathrm{\top }}\hat{\epsilon }/(n-K-1)}$$

其中，

${\hat{\epsilon }}^{\mathrm{\top }}\hat{\epsilon }/(n-K-1)$ 就是

$s^2$。則，

$$F\sim F(q,n-K-1)$$

考慮一個 Monte Carlo 模擬實驗。令

$\mathbf{y}=\mathbf{X}\beta +\epsilon$，其中

$$\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],\epsilon \sim \mathcal{N}(\mathbf{0},0.25\times \mathbf{I})$$

# 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

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}$，則

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 }_{be}\\ {\beta }_{bw}\\ {\beta }_{bill}\\ {\beta }_{om}\end{array}\right],\mathbf{r}=\left[\begin{array}{l}0\\ 0\\ 0\\ 0\end{array}\right]$$

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

代表我們可以拒絕虛無假設，其實際值小於

$\text{verb}|R|$ 的寬容值

$2.2\times {10}^{-16}$。一個簡單的計算方式是：

$$F=\frac{\sum _{i=1}^n{({\hat{y}}_i-\overline{y})}^2/K}{\sum _{i=1}^n{(y_i-{\hat{y}}_i)}^2/(n-K-1)}=\frac{SSE}{SSR/(n-K-1)}$$

可以看出

$F$ 檢定統計量與

$R^2$ 有一定程度的（非線性）關係：

$R^2$ 越大，

$F$ 檢定統計量也越大。而我們可以建構一個多元迴歸模型的 ANOVA（變異數分析）表進行計算分類後的資料是否能夠良好的詮釋模型。

Source	SS	df	MS	$\varphi$
Factor	SSR	$k$	MSR = $\frac{SSR}{k}$	${\varphi }_F=\frac{MSR}{MSE}$
Error	SSE	$n-k-1$	MSE = $\frac{SSE}{n-k-1}$
Total	SST	$n-1$

其中檢定統計量

${\varphi }_F=\frac{MSR}{MSE}\sim \mathrm{F}(1,n-2)$。而在

$\text{verb}|R|$ 中，我們可以用 anova() 進行變異數分析。

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$ 的不偏估計式，而只要殘差與誤差項夠接近，就可以印證上述的假設為真。首先我們先將誤差項進行標準化，得到

$${\epsilon }_i\stackrel{i.i.d}{\sim }\mathcal{N}(0,{\sigma }^2)$$

並由小排到大。如果我們有

$n$ 個服從常態分配的殘差與

$n$ 個經過標準化後服從常態分配的誤差項，分別將其畫在

$x$ 軸與

$y$ 軸，此兩者應該要分布在

${45}^{\circ }$ 線上，形成 QQPlot。^[3]舉例來說，

X <- runif(1000)
Y <- 1 + X + 0.2 * rnorm(1000)
lm <- lm(Y~X)
plot(pizza.lm, which = 2)

如果誤差項不是服從常態分配，而是服從均勻分配的話，QQPlot 則會長得像下面那樣，並非在

${45}^{\circ }$ 線上。

X <- runif(1000)
Y <- 1 + X + 0.2 * runif(1000) - 0.5
lm <- lm(Y~X)
plot(pizza.lm, which = 2)

最後，我們可以畫出 pizza 中誤差項的 QQPlot。

pizza.lm <- lm(time~branch + bill + operator, data = pizza)
plot(pizza.lm, which = 2)

但是，如果誤差項真的不是服從常態分配，我們還能夠繼續進行上述的分析嗎？答案是可以的！事實上，利用漸進理論──只要樣本數夠大，我們還是能夠維持原本的假設。

---

1. 當2個（或以上）的自變數互不獨立(即彼此相關)，就是具有「線性重合」，將會使迴歸模型中存在著重複的自變數，提高某一自變數的解釋力與預測力，使得理論的建構不正確。 ↩︎

2. 如果我們要解出

   $\mathbf{Ax}\mathbf{=}\mathbf{b}$ 的話，我們可以輸入 solve(A, b)，便會算出
   $\mathbf{x}\mathbf{=}{\mathbf{A}}^{\mathbf{-}\mathbf{1}}\mathbf{b}$；如果我們僅輸入 solve(A)，
   $\mathtt{\text{R}}$便會自動認定我們要計算
   $\mathbf{Ax}\mathbf{=}\mathbf{I}$，得到
   $\mathbf{x}\mathbf{=}{\mathbf{A}}^{\mathbf{-}\mathbf{1}}\mathbf{I}\mathbf{=}{\mathbf{A}}^{\mathbf{-}\mathbf{1}}$。當然，看到上面的計算過程，對於線性代數理解夠的人就可以知道， solve() 可以拿來取反矩陣(inverse matrix)。 ↩︎

3. plot() 中的指令 which = 2 代表要畫出 QQPlot。 ↩︎