1 課程目標

本週上課將介紹迴歸的基本原理，首要目標是描述兩個變數之間的關係。例如馬拉松的參賽者完賽時間與年齡這兩個變數的散佈圖，加上迴歸線：

library(UsingR)
ggplot(nym.2002, aes(age, time)) +
       geom_point(col="#808000", size=2) +
       geom_smooth(method="lm")  +
       theme_bw() +
       theme(axis.title=element_text(size=12), 
             axis.text = element_text(size=12))

Figure 1.1: 迴歸線與散佈圖

從上面的迴歸線以及散佈圖1.1可看出，年齡越大，完成比賽所需的時間越長。
根據這樣的結果，我們可以預測某個年齡的平均完成時間，例如有一群45歲的人參賽，我們可以預測大概需要250分鐘跑完。

2 迴歸的意義

迴歸的目的是找到依變數或是反應變數(outcome)Y與一個或是一個以上的預測變數(predictor)X的函數關係。我們希望得到一個X與Y的關係，預測當某個事件發生或者某個數值出現時，另外一個變數發生的機率或者平均數。所以預測是迴歸模型的主要目的。

迴歸的基礎是條件機率，也就是迴歸線逼近在X成立的條件下，Y的平均值。所以模型的對象是\(E[Y|X]\)，也就是X變化時，Y的平均值變化大小。或者是說固定X情況下，Y的平均值大小。

假設我們全查某個母體的所有個案，例如所有OECD的會員國，母體迴歸函數(Population regression function)表示為：

\[\begin{align} Y=\beta_{0}+X\beta_{1}+\epsilon \tag{2.1} \end{align}\]

或者

\[ Y=f(X)+\epsilon \]

其中，\(Y\)是依變數，\(X=(X_{1},X_{2},\ldots,X_{p})\)為自變數。\(f\)代表固定但是未知的函數，也就是\(Y\)與\(X_{1},X_{2},\ldots,X_{p}\)之間的關係。\(\epsilon\)是隨機變數，稱為誤差項。可能包含我們沒有測量到的變數，或是無法測量的變數，並且獨立於\(X\)。例如，不同病情的病人對於某種藥劑有不同的反應，兩者之間的關係可以用\(Y=f(X)\)表示。但是這種藥劑製作過程中，可能受到某些因素的影響而有不同的藥效，也可能因為病人當天的生理狀況而有不同的反應，這些都可歸類為\(\epsilon\)。
當我們把\(Y\)改成\(E[Y|X]\)，就是改為計算\(Y\)的條件期望值，方程式 (2.1)可以改成條件期望值函數：

\[\begin{align} E[Y|X]=\beta_{0}+\beta_{1}X+E[\epsilon|X] \tag{2.2} \end{align}\]

我們假設\(E[\epsilon|X]=0\)，方程式(2.2)解釋成當X固定某一個值，Y的平均值是由\(\beta_{0}\)與\(\beta_{1}\)決定。而個別\(Y_{i}\)的預測值則寫成：

\[ \hat{Y_{i}}=\hat{\beta_{0}}+\hat{\beta_{1}}X_{i} \] - 如果\(E[\epsilon|X]\neq 0\)，表示誤差項和自變數之間有相關，這樣會造成\(\hat{\beta_{0}}\)與\(\hat{\beta_{0}}\)的估計有誤。

迴歸模型就是盡可能逼近X各類別下Y的平均數的函數。一旦找到這樣的迴歸函數，可以用來預測其他的觀察值。

近年來熱門的統計學習(statistical learning)包含許多估計\(f\)的方法。統計學習的目的是：

預測：假設X與Y的關係是\(Y=f(X)+\epsilon\)，並且假設誤差項 \(\epsilon\) 的平均值為0，並且定義\(\hat{Y}=\hat{f}(X)\)。

因為\(\hat{f}(x)\)只是\(f(x)\)的眾多預測之一，所以\(f(x)\)與\(\hat{f}(x)\)之間有誤差(bias)。
即使\(\hat{f}\)等於\(f\)，\(\hat{Y}=f(X)+\epsilon\)，因此讓\(E(\epsilon^2)=Var(\epsilon)\)越小，預測就越準確。不過，我們無法減少\(\epsilon\)，因此預測的正確性(accuracy)有上限，實際上無法得知上限有多少。

推論：我們關心的是\(X\)與\(Y\)之間的關係，\(Y\)如何隨著\(X\)而變動。例如，某家公司收集了許多消費者資料，他們想要知道那些背景的消費者對於他們的產品有興趣，也就是從許多變數中找到與\(Y\)相關的變數。此外，\(X\)與\(Y\)有可能是非線性的關係。線性模型的推論比較容易詮釋，但是不一定能預測所有觀察值。非線性模型的預測能力比較好，但是不容易詮釋。
預測與推論兩者並重：一方面估計哪些房屋的特徵會影響房價，一方面預測特定房屋的價格被低估或者高估。

迴歸模型常見的名詞

參數：\(E[Y|X]\)是觀察不到的母體參數，所以模型的對象\(E[Y|X]\)。因為我們只有一組樣本，所以得到的只是\(E[Y|X]\)的估計其中之一。

樣本統計表示為：\(\hat{E}[Y|X]\)，\(X\)成立的條件下，\(Y\)的期望值。

估計：某一個樣本所得到的\(\hat{E}(Y|X)\)。

統計上兩個重要概念

estimand：我們關心的統計對象，例如母體的平均數\(\mu=\mathbb{E}[X]\)，ATE=\(\mathbb{E}[Y(1)-Y(0)]\)。
estimator: 統計的方式，例如母體的平均數\(\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}X_{i}\)，ATE則是\(\hat{\tau}=\frac{1}{n_{1}}\sum_{i,D_{i}=1}Y_{i}-\frac{1}{n_{0}}\sum_{i,D_{i}=0}Y_{i}\)
而在迴歸模型，\(Y=\beta_{0}+\beta_{1}X+\epsilon\)，其中\(\beta\)是estimand，而estimator是\(\hat{\beta}=\frac{\sum (X_{i}-\bar{X})\sum (Y_{i}-\bar{Y})}{\sum (X_{i}-\bar{X})^2}\)。

3 有母數統計：線性迴歸

有母數分析方法或者檢定假設母群體為常態分布以進行檢定的方式，包括變異數分析、\(t\)檢定。無母數分析方法或者檢定並沒有假設，可用在未知的分佈，或是樣本數不夠大。無母數迴歸模型相對比較有彈性，但是不適用預測新的觀察值。

線性迴歸假設\(Y\)的條件平均值呈現線性關係，如同方程式(2.2)：

\[ \mathbb{E}[Y|X]=\beta_{0}+X\beta_{1} \]

上面的模型只有兩個係數。\(\beta_{0}\)稱為截距或常數，\(\beta_{1}\)稱為斜率係數。當X=0，\(E[Y|X=0]\)等於Y的期望值，也就是\(\beta_{0}\)。或者是當X減去自己的平均數，\(E(X)\)也會等於0，\(E[Y]=\beta_{0}\)。

因為平均數是一個常數，常數的期望值等於該常數，也就是：

\[ \mathbb{E}[\bar{X}]=\bar{X} \]

所以當\(X^{*}=X-\bar{X}\)，代入：

\[\begin{align*} E[Y|X^{*}] = \beta_{0}+\beta_{1}X^{*} \\ & = E[\beta_{0}]+E[\beta_{1}\times (X-\bar{X})] \\ & = E[\beta_{0}]+E[\beta_{1}\times (X-\bar{X})] \\ &= \beta_{0}+\beta_{1}\times E[(X-\bar{X})] \\ &= \beta_{0}+\beta_{1}\times (E[X]-E[\bar{X}]) \\ &= \beta_{0}+\beta_{1}\times (\bar{X}-\bar{X}) \\ & = \beta_{0} \tag{3.1} \end{align*}\]

方程式(3.1)顯示，當自變數減去自己平均值後，依變數的期望值等於截距\(\beta_{0}\)。這個結果顯示\(\beta_{1}=0\)，也就是依變數\(Y\)與\(Ｘ\)無關，這個時候\(Y\)的平均值等於\(\beta_0\)。圖??顯示\(X\)減去平均數之後，\(Y\)的期望值就是常數項。線性假設的其中之一是截距為固定。

使用觀察資料，或者訓練資料估計模型之後，可以得到模型的預測值，表示為：

\[\hat{y}=\hat{\beta}_{0}+\hat{\beta}_{1}x \]

\(y_{i}\)與\(\hat{y}_{i}\)之間的差異稱為殘差\(e_{i}\)，或者residual。

定義殘差的平方和為：

\[RSS=e_{1}^2+e_{2}^2+\ldots +e_{n}^2 \]

RSS越小，表示模型跟資料之間的差距越小，也就是模型符合資料。RSS除以觀察值\(n\)之後等於MSE。

3.1 什麼是E[Y|X]？

E[Y|X]表示固定X=x時，Y的期望值（平均值）。\(E[Y|X]\)稱為conditional expectation function(CEF)。例如，\(E(Y|X=1)=\beta_{0}+\beta_{1}X\)。

圖3.1表示CEF的Y的條件期望值落在迴歸線上（參考University of Memphis）：

Figure 3.1: Y的條件期望值

當X是隨機變數，E[Y|X=x]也是隨機變數，E[Y|X]隨著X變化而改變，可以寫成h(Y)。也就是知道X可以幫助我們預測Y。

因此，E(Y|X=x)就像Y一樣，有期望值也有變異數。當X是連續變數，期望值可用積分表示為：

\[E[X]=\int_{-\infty}^{\infty}xf(x)dx\]

當X是不連續變數，期望值可用總和表示為：

\[E[X]=\sum xf(x)=\sum x\cdot P(X=x)\]

變異數則是：

\[\rm{Var}(X)=E(X^2)-(E(X))^2\]

3.1.1 期望值

E[Y|X]的期望值為：

\[\boxed{E(Y|X=x)=\mu_{Y|X=x}=\sum yf(y|x)}\]

我們用以下的例子說明如何計算Y的條件期待值（見Rachel Fewster）：
假設X是一個隨機變數，有如下的機率分佈(3.2)：

\[\begin{equation} X = \begin{cases} 1 & \rm{with\hspace{.3em} probability}\quad 1/8\\ 2 & \rm{with\hspace{.3em}probability}\quad 7/8 \end{cases} \tag{3.2} \end{equation}\]

而且Y與X的關係如(3.3)：

\[\begin{equation} Y|X= \begin{cases} 2X & \rm{with\hspace{.3em}probability}\quad 3/4\\ 3X & \rm{with\hspace{.3em}probability}\quad 1/4 \end{cases} \tag{3.3} \end{equation}\]

如果X=1，

\[\begin{equation} Y|(X=1)= \begin{cases} 2 & \rm{with\hspace{.3em}probability}\quad 3/4\\ 3 & \rm{with\hspace{.3em}probability}\quad 1/4 \end{cases} \tag{3.4} \end{equation}\]

因此，E(Y|X=1)=\(2\times \frac{3}{4}+3\times\frac{1}{4}=\frac{9}{4}\)。\(E(Y|X=2)=\frac{18}{4}\)，依此類推。

3.1.2 變異數

conditional expectation function的變異數寫成：

\[\boxed{Var(Y|X=x)=E(Y^2|X)-(E(Y|X))^2}\]

根據(3.4)，計算\(\text{Var}(Y\mid X=1)\)需要求出\(E(Y^2|X)\)以及\((E(Y|X))^2\)，其中\(E(Y^2|X)\)等於：

\[\begin{equation} Y^2|(X=1)= \begin{cases} 4 & \rm{with\hspace{.3em}probability}\quad 3/4\\ 9 & \rm{with\hspace{.3em}probability}\quad 1/4 \end{cases} \tag{3.5} \end{equation}\]

\(E[Y^2\mid X=1]=4\times \frac{3}{4}+9\times \frac{1}{4}=\frac{21}{4}\)
\(E[Y\mid X=1]^2=\frac{9^2}{4^2}=\frac{81}{16}\)
因此，\(\text{Var}(Y\mid X=1)=\frac{21}{4}-\frac{81}{16}=\frac{3}{16}\)

3.1.3 Law of Iterated Expectations

條件期待函數有一個特性是條件期待值的期待值剛好等於該自變數的期待值。這是因為對於每一個自變數的值，都會有一個發生的機率。條件期望值乘上每一個X發生的機率加起來等於\(Y\)的期望值：

\[E[Y]=\sum E[Y|X=x]P(X=x)]\]

可以推導出：

\[E(Y)=E[E(Y|X)]\]

稱為Law of Iterated Expectations(LIE)。

3.1.4 範例

例如某個民調機構連續三天訪問民眾，詢問他們早餐喝什麼飲料，回答的分佈如表 3.1：

Table 3.1: Distribution of Drinks
	Tea	Coffee	Total
Monday	200	100	300
Tuesday	280	180	460
Wednesday	240	200	440

請問平均有多少人喝茶？喝茶的機率多少？
直覺來看，我們可以把喝茶的人加總起來，也就是表??的第一欄的總和700，然後除以1200，等於0.6。換句話說，有60%的人喝茶。
我們可以用邊際機率以及應用LIE計算平均喝茶的機率如方程式 (3.6)：

\[\begin{align} E[T] = E[E(Tea|Day)]\\ & = E(T|Day=\text{Monday})\cdot P(\text{Monday}) + E(T|Day=\text{Tuesday})\cdot P(\text{Tuesday}) + E(T|Day=\text{Wednesday})\cdot P(\text{Wednesday}) \\ & = \frac{200}{300}\cdot \frac{300}{1200}+ \frac{280}{460}\cdot \frac{460}{1200}+ \frac{240}{440}\cdot \frac{440}{1200}\\ & = 0.6 \tag{3.6} \end{align}\]

Table 3.2: Distribution of Drinks
	Tea	Coffee	Total
Monday	200	100	300
Tuesday	280	180	460
Wednesday	240	200	440
Total	720	480	1,200

LIE也可以想成是加權平均值的平均值，例如有兩班學生量身高，我們可以全部學生的身高總和除以兩個班學生人數得到全體學生的身高平均數，也可以先計算第一班的身高平均數乘以該班佔全部學生的比例，然後同樣方法計算第二個班，相加之後就是全部學生的平均身高。

\[ E[Y] = E[Y|X=1]\times P(X=1)+E[Y|X=2]\times P(X=2)=E[E[Y|X]] \]

如果我們只知道兩個班裡面男生與女生的平均身高，以及男生與女生在這兩個班的比例，我們可以根據LIE計算全部的平均身高。

我們模擬資料來說明LIE。首先創造兩個變數x, y。

set.seed(02138)
#X, Y
x<-rnorm(1000, 0, 0.1)
y<-rnorm(1000, 0, 1) + x

然後計算E[E(Y|X)]，方法是利用dplyr的\(\tt{group\_by}\)以及\(\tt{summarise}\)兩個功能。然後計算原始y的平均數，結果兩者一模一樣。

library(dplyr)
tmp <- data.table::data.table(X=x, Y=y) %>%
   group_by(X) %>%
   summarise(y.bar=mean(Y))
mean(tmp$y.bar)

[1] 0.04224

mean(y)

[1] 0.04224

4 最小平方法

OLS代表Ordinary Least Squares，也就是最小平方法。

假設\(\{\tilde{\beta}_{0},\tilde{\beta}_{1}\}\)是\(\{\beta_{0},\beta_{1}\}\)的可能值之一，OLS迴歸的目的是找到最小化以下模型的\(\{\hat{\beta}_{0},\hat{\beta}_{1}\}\)

：

\[\begin{align*} \{\tilde{\beta}_{0},\tilde{\beta}_{1}\}&=\arg\min\limits_{{\tilde{\beta}_{0},\tilde{\beta}_{1}}}\sum _{i=1}^n(y_{i}-\tilde{\beta}_{0}-x_{i}\tilde{\beta}_{1})^2\\ &= \arg\min\limits_{{\tilde{\beta}_{0},\tilde{\beta}_{1}}}\sum _{i=1}^n\hat{u_{i}}^2 \end{align*}\]

OLS迴歸的目標就是最小化殘差。因為殘差容易分析而且計算，並且在某些前提成立時可達到最佳化，因此我們分析最小化殘差。

Figure 4.1: 迴歸線圖

5 迴歸係數的意義

迴歸模型\(E[Y|X]=\hat{\beta}_{0}+\hat{\beta}_{1}X\) 中，\(\hat{\beta}_{1}\)的意義為何？

\(\beta_{1}>0\): \(E[Y|X]\) 隨著X增加
\(\beta_{1}<0\): \(E[Y|X]\) 隨著X減少
\(\beta_{1}＝0\): \(E[Y|X]\)與X無線性關係

如果用在預測：

估計樣本的\(\hat{\beta}_{0}\)以及\(\hat{\beta}_{1}\)。
根據\(\hat{E}[Y|X=x_{new}]=\hat{\beta}_{0}+\hat{\beta}_{1}x_{new}\)得到新的預測值\(\hat{E}[Y|X=x_{new}]\)
\(\hat{E}[Y|X=x_{new,1}]=\hat{\beta}_{0}+\hat{\beta}_{1}(x_{new,1})\)
\(\hat{E}[Y|X=x_{new,2}]=\hat{\beta}_{0}+\hat{\beta}_{1}(x_{new,2})\)
\(\hat{E}[Y|X=x_{new,1}]-\hat{E}[Y|X=x_{new,2}]=\hat{\beta}_{1}(x_{new,1}-x_{new,2})\)

5.1 迴歸模型：電視廣告與銷售量

電視廣告與銷售量之間的關係以迴歸表示為：

sales \(\approx {\beta}_{0}+{\beta}_{1}\cdot\) TV

用 \(\tt{lm()}\)估計迴歸方程式如5.1：

file<-here::here('data','advertising.csv')
Advertising<-read.csv(file, sep=',', header = TRUE)
OLS1<-lm(sales ~ TV, data=Advertising)

stargazer::stargazer(OLS1,  title='電視廣告與銷售量之間的關係',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 5.1: **電視廣告與銷售量之間的關係**

	Dependent variable:

	sales

TV	0.048^***
	(0.003)

Constant	7.033^***
	(0.458)


Observations	200
R²	0.612
Adjusted R²	0.610
Residual Std. Error	3.259 (df = 198)
F Statistic	312.100^*** (df = 1; 198)

Note:	p<0.1; p<0.05; p<0.01

資料來源：James et al., 2013
電視廣告與銷售量之間的關係可寫成sales = \(7.03+0.04\cdot\) TV
但是觀察到的sales與迴歸線預測的之間會有差異。
畫圖5.1如下：

library(ggplot2)
A1<-ggplot(Advertising, aes(x=TV, y=sales)) +
  labs(y = 'Sales', x = 'TV')  +
  geom_point(col="saddlebrown") +
  geom_smooth(method="lm", col='blue', se=F, size=1.5) ; A1

Figure 5.1: 電視廣告與銷售量之間的關係

加上預測值以及觀察值與預測值之間的殘差（參考BLOGR）：

m1 <- lm(sales ~ TV, Advertising)
Advertising$predicted <- predict(m1)   # Save the predicted values
Advertising$residuals <- residuals(m1) # Save the residual values

# Quick look at the actual, predicted, and residual values
#library(dplyr)
#Advertising %>% select(predicted, residuals) %>% head()
ggplot(Advertising, aes(x=TV, y=sales)) +
    geom_point(color="saddlebrown") +
    geom_point(aes(y = predicted), shape = 1, color="red") +
  geom_segment(aes(xend=TV, yend=predicted), color="lightgray") +
    theme_bw()

如果殘差的總和越小，預測值與觀察值就越接近，所得到的迴歸係數以及整個模型就越能預測未來的依變數。

6 最小平方法估計迴歸模型

假設有\(n\)個樣本，表示為\(y_{1},y_{2},\cdots,y_{n}\)。

有一個可以摘要說明\(y_{1},y_{2},\cdots,y_{n}\)的統計稱為\(\tilde{u}\)，用來估計\(\hat{u}\)。

\(\tilde{u}\)可以被定義成觀察值與預測值的差距：

\[\begin{align*} \tilde{u} & = y_{i}- \hat{y}_{i} \\ & = y_{i}- \tilde{\beta}_{0}-\tilde{\beta}_{1}x_{i} \end{align*}\]

OLS迴歸的目標是最小化殘差，殘差可以計算為\(y_{i}- \hat{y}_{i}\)，但是這會給所有的\({y}_{i}\)一樣的權重。改用殘差的平方和(sum of squared residuals, 才會給離迴歸線越遠的\({y}_{i}\)比較大的權重，這樣才能反映離差的大小。

SSR=\(\sum \tilde{u}^2\)，SSR越小越好。

命SSR為函數\(S(\tilde{u})\)。要找出最小化SSR的\(\tilde{u}\)

\[\begin{align} S(\tilde{u}) = {\sum\limits_{i=1}^n(y_{i}-\hat{y})^2} \\ L(y) = {\sum\limits_{i=1}^n(y_{i}-\hat{y})^2} \tag{6.1} \end{align}\]

我們稱\(L(y)\)為loss function，中文稱為損失函數。方程式@(eq:tildeu)又稱為Mean Squared Error (MSE)，如果開根號稱為Root Mean Squared Error (RMSE)，也稱為prediction loss。而在估計\(\beta\)的過程中，我們同樣地要盡量減少\(\hat{\beta}\)與\(\beta\)之間的差距，所以： \[ L(\hat{\beta}, \beta)=(\hat{\beta}-\beta)^2 \]

6.1 推導過程

對(6.1)進行微分：

\[\begin{align*} \frac{dL}{d\hat{y}}L(y) & = \frac{d}{d\hat{y}}{\sum\limits_{i=1}^n(y_{i}-\hat{y})^2} \\ & = \sum\limits_{i=1}^n(y_{i}^2-2y_{i}\hat{y}+\hat{y}^2) \\ & = \sum(-2y_{i}+2\hat{y}) \\ & = -2\sum(y_{i} - \hat{y}) \end{align*}\]

進行微分之後，因為要求出最小化的線性函數，所以設定

\[-2\sum(y_{i} - \hat{y})=0\]

\[\begin{align*} -2\sum(y_{i} - \hat{y}) & = -2(\sum y_{i})-n\hat{y}\\ n\hat{y} = \sum y_{i} \end{align*}\]

樣本平均數\(\frac{\sum y_{i}}{n}\)會讓殘差最小，也就是最小平方法的樣本統計(estimator)。換句話說，對特定的X，最小平方法找到的迴歸係數所產生的殘差最小，而該預測值剛好等於Y的平均值，參考Sophie Bennett。把這些平均值連起來就是迴歸線。

如何用模擬資料證明這個理論？如果我們隨機產生一筆資料以及一筆任意值，函數為兩者之間的平方差的和，如果是該筆資料的期望值最小化該函數，那麼我們就證明這個理論。

圖6.1顯示，當\(\tilde{u}=y-\hat{y}=\bar{y}\)=E(Y|X)時，y與其他觀察值之間的誤差平方和會最小。
我們稱這個損失函數為sum of squared of the errors (SSE)或者SSR。
我們可以應用SSR到模型適合度\(R^2\)，容易詮釋，所以廣被使用。進一步說明可見Dustin Stansbury的文章。

#create variable 1000 Y with mean=10
set.seed(666)
Y <- rnorm(1000, mean = 10, sd = 2)

#create 50 data points with mean = 10
obs <- rnorm(20, 10, 3)

# generate empty list
squared_residuals <- rep(NA, length(obs))
min_residual <- Inf
min_observation <- NULL

#calculate the sum of squared residuals between Y and obs
for (i in 1:length(obs)) {
  squared_residuals[i] = sum((Y - obs[i])^2)
  if (squared_residuals[i] < min_residual) {
         min_residual <- squared_residuals[i]
         min_observation <- obs[i]
}
}

cat('Minimization of Squared Residuals=', min_observation)
## Minimization of Squared Residuals= 9.893

#plot the sum of squared residuals along with the 50 observations to find the lowest sum of squared residuals
data.frame(obs_values = obs, squared_residuals = squared_residuals) %>% 
  ggplot(aes(obs_values, squared_residuals)) +
  stat_smooth(method="lm", 
              formula = y ~ poly(x, 2), 
              se = FALSE,
              colour = "#FCC3B0",
              linetype = "dashed") +
  geom_point(col = "#661E4F", size = 2.2, alpha = 0.7) +
  geom_vline(xintercept =  min_observation, linetype='dashed', colour='#EEAA11') +
  geom_hline(yintercept=min(squared_residuals), linetype='dashed', colour='#113BEA') +
  labs(title = "Sum of Squared Residuals (SSR) Loss Function",
     x = "Summary Value",
     y = "SSR") +
  theme_minimal() +
 # theme(text = element_text(family = ""),
 # plot.title = element_text(family = "", hjust = 0.5) +
  scale_y_continuous(labels = scales::comma)

Figure 6.1: 誤差平方和函數圖

6.1.1 最小平方法的係數

我們要找出可以最小化殘差的\(\{\tilde{\beta}_{0},\tilde{\beta}_{1}\}\)，做為\(\{\beta_{0},\beta_{1}\}\)的估計。也就是Linear Least Squares Estimator。方程式表示如下：

\[\begin{align*} S(\tilde{\beta}_{0},\tilde{\beta}_{1})= {\sum\limits_{i=1}^n(y_{i}-\tilde{\beta}_{0}-x_{i}\tilde{\beta}_{1})^2} \end{align*}\]

步驟：

1.首先對\(\{\tilde{\beta}_{0},\tilde{\beta}_{1}\}\)取微分

2.命兩者的微分等於0

3.求出\(\{\tilde{\beta}_{0},\tilde{\beta}_{1}\}\)的解。

因為要求出最小化殘差的線性函數，而殘差可轉換成：

\[y_{i}-\tilde{\beta}_{0}-x_{i}\tilde{\beta}_{1}\]

所以：

\[\begin{align*} S(\tilde{\beta}_{0},\tilde{\beta}_{1}) & = {\sum\limits_{i=1}^n(y_{i}-\tilde{\beta}_{0}-x_{i}\tilde{\beta}_{1})^2} \\ & = {\sum\limits_{i=1}^n(y_{i}^2-2y_{i}\tilde{\beta}_{0}-2y_{i}\tilde{\beta}_{1}x_{i}+ \tilde{\beta}_{0}^2+2\tilde{\beta}_{0}\tilde{\beta}_{1}x_{i}+ \tilde{\beta}_{1}^2x_{i}^2)} \end{align*}\]

\[\begin{align*} \frac{\partial S(\tilde{\beta}_{0},\tilde{\beta}_{1})}{\partial \tilde{\beta}_{0}}= {\sum\limits_{i=1}^n(-2y_{i}+2\tilde{\beta}_{0}+2\tilde{\beta}_{1}x_{i})} \end{align*}\]

\[\begin{align*} \frac{\partial S(\tilde{\beta}_{0},\tilde{\beta}_{1})}{\partial S(\tilde{\beta}_{1})}= {\sum\limits_{i=1}^n(-2y_{i}x_{i}+2\tilde{\beta}_{0}x_{i}+2\tilde{\beta}_{1}x_{i}^2)} \end{align*}\]

6.1.2 first order condition

進行微分之後，因為要求出最小化的線性函數，所以分別設定

\[\begin{align*} 0 & = {\sum\limits_{i=1}^n(-2y_{i}+2\tilde{\beta}_{0}+2\tilde{\beta}_{1}x_{i})} \tag{6.2} \end{align*}\]

\[\begin{align*} 0 & = {\sum\limits_{i=1}^n(-2y_{i}x_{i}+2\tilde{\beta}_{0}x_{i}+2\tilde{\beta}_{1}x_{i}^2)} \tag{6.3} \end{align*}\]

這兩個方程式被稱為 first order condition。

6.1.3 normal equations

根據(6.2)，

\[\begin{align*} \hat{\beta}_{0}n &= {(\sum\limits_{i=1}^ny_{i})- \hat{\beta}_{1}{(\sum\limits_{i=1}^nx_{i})}} \\ \tilde{\beta}_{0} &=\bar{Y}-\tilde{\beta}_{1} \bar{x} \tag{6.4} \end{align*}\]

根據(6.3)，

\[\begin{align*} \hat{\beta}_{1}\sum\limits_{i=1}^nx_{i}^2 &= (\sum\limits_{i=1}^nx_{i}y_{i})- \tilde{\beta}_{0}{(\sum\limits_{i=1}x_{i}}) \tag{6.5} \end{align*}\]

這兩個等式稱為normal equations。
從normal equations的第1式((6.4)可以得到：

\[\begin{align*} \hat{\beta}_{1} &=\frac{\sum y_{i}-\hat{\beta}_{0}n}{\sum x_{i}} \tag{6.6} \end{align*}\]

根據式(6.5)則可解\(\hat{\beta}_{0}\)為：

\[\begin{align*} \hat{\beta}_{0}=\frac{\sum x_{i}y_{i}-\hat{\beta}_{1}\sum x_{i}^2}{\sum x_{i}} \tag{6.7} \end{align*}\]

把方程式(6.7)代入(6.6):

\[\begin{align*} \hat{\beta}_{1} & =\frac{\sum y_{i}-n\frac{\sum x_{i}y_{i}-\hat{\beta}_{1}\sum x_{i}^2}{\sum x_{i}}}{\sum x_{i}}\\ & = \frac{n\sum x_{i}y_{i}-\sum x_{i}\sum y_{i}}{n\sum x_{i}^2-(\sum x_{i})^2} \tag{6.8} \end{align*}\]

在(6.8)的分子可以改寫為，

\[\begin{align*} {n\sum x_{i}y_{i}-\sum x_{i}\sum y_{i}} & = {n(\sum x_{i}y_{i}-\bar{y}\sum x_{i})}\\ & = {n(\sum x_{i}y_{i}-n\bar{x}\bar{y})} \\ & = {n(\sum x_{i}y_{i}-n\bar{x}\bar{y}-n\bar{x}\bar{y}+n\bar{x}\bar{y})} \\ & = {n(\sum x_{i}y_{i}-\bar{y}\sum x_{i}-\bar{x}\sum y_{i}+\sum \bar{x}\bar{y})} \\ & = {n\sum (x_{i}-\bar{x})(y_{i}-\bar{y})} \tag{6.9} \end{align*}\]

在(6.8)的分母可以改寫為，

\[\begin{align*} {n\sum x_{i}^2-(\sum x_{i})^2} & = {n\sum x_{i}^2-(\sum x_{i})^2+(\sum x_{i})^2-(\sum x_{i})^2}\\ & = {n\sum x_{i}^2-2n\bar{x}(\sum x_{i})+n^2(\bar{x})^2} \\ & = {n(\sum x_{i}^2-2\bar{x}\sum x_{i}+n(\bar{x})^2)}\\ & = {n(\sum x_{i}^2-2\bar{x}\sum x_{i}+\sum \bar{x}^2)}\\ & = {n\sum(x_{i}-\bar{x})^2} \tag{6.10} \end{align*}\] （因為\(\sum \bar{x}^2=\bar{x}^2+\bar{x}^2+\cdots +\bar{x}^2=n\bar{x}^2\)）

6.2 OLS與共變數

根據(6.9)以及(6.10)，式(6.8)可以寫成：

\[\begin{align*} {\hat{\beta}_{1}=\frac{\sum_{i=1}^n(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum_{i=1}^n(x_{i}-\bar{x})^2} }=\frac{\mathrm {Sample\hspace{.5em} covariance\hspace{.5em} of\hspace{.5em} X\hspace{.5em} and\hspace{.5em} Y}}{\mathrm {sample\hspace{.5em} variance\hspace{.5em} of\hspace{.5em} X}}\\ \end{align*}\]

X的變異數越大，\(\hat{\beta}_{1}\)越小
X, Y的共變量越大，\(\hat{\beta}_{1}\)越大。共變量可寫成\(\sum xy\)，其中 \(x=(X-\bar{X})\)，\(y=(Y-\bar{Y})\)。
另外，

\[ \hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1}\bar{x} \]

★如果使用圖書館人數為X，借出的書本（百本）為Y，而且知道過去25天的各項變數統計為\(\sum X=200\)、\(\sum Y=300\)、\(\sum X^2=16600\)、\(\sum Y^2=1696\)、\(\sum XY=2436\)，請求出\(\hat{Y}=\hat{\beta_{0}}+\hat{\beta_{1}X}\)。

\(\blacksquare\)把以上的數據代入(6.8)，分別求出\(\beta_{1}\)以及\(\beta_{0}\)：

n=25
#sum(X)
sum.x=200
#sum(Y)
sum.y=300
#sum(XY)
sum.xy <- 2436
#sum(x^2)
sumxsq <- 1660
#sum(y^2)
sumysq <- 1696
#(sumX)^2
sumx.sq <- sum.x^2
beta1 <- (n*sum.xy-sum.x*sum.y)/(n*sumxsq-sumx.sq)
beta1

[1] 0.6


beta0=sum.y/n-beta1*sum.x/n
beta0

[1] 7.2

6.3 從微分看迴歸係數的意義

假設\(y=f(x)\)，在\(\mathscr{Y}\)的每一個值，來自於\(\mathscr{X}\)的值。

理論上\(f\)函數可以是\(n\)度空間的實數對應\(\mathscr{X}\)的實數。但是我們限定在1度空間，也就是\(f:\Re\rightarrow \Re\)

圖6.2顯示，函數\(y=f(x)=x^2-x\)在\(x=(2,3)\)之間，\(y\)從2(\(2^2-2\))上升到6(\(3^2-3\))，那麼這個函數的上升程度可以用\(\frac{6-4}{3-2}=2\)來表示。圖6.2顯示一個一元二次函數，當x=2時，y=2，當x=3，y=6，以此類推。
我們可以從X=\(\xi\)延伸畫出許多條穿過\(f(x)\)的割線，只有一條線通過X=\(\xi\)但是沒有通過第二個點，稱為切線。換句話說，當許多函數上的點趨近於X=2時，割線也趨近於切線，所以切線的斜率發生在X=\(\xi\)，如圖6.3。
為什麼我們要求出切線的斜率？因為我們想知道X對應到這個曲線上的每一個點的Y變化速度(rate of change)。曲線上的頂點跟底點都代表沒有變化，以圖6.2而言，在到達頂點之前，X增加一個單位Y會減少一些單位，因為是曲線，所以每一個X增加的單位，Y減少的幅度會不同。同樣的，在到達底點之後，X增加一個單位Y會增加一些單位。切線與曲線交會的點，代表X增加或減少一個單位，Y同時有最大增加或減少的程度，我們如果知道切線的斜率，那麼就知道X=x時，對Y的最大影響。

## -x + x^2

Figure 6.2: 一元二次函數求微分圖1

如果我們把數字替換成變數，圖6.3顯示，函數\(y=f(x)\)從點\((\xi)\)增加到\((\xi+h)\)時，\(y\)從\(f(\xi)\)增加到\(f(\xi+h)\)。

Figure 6.3: 一元二次函數求微分圖2

函數的斜率是\({m=\frac{\mathit{f(\xi+h)-f(\xi)}}{\mathit{\xi+h-\xi}}} \quad h>0\)。
如果用微分方式，對於可微分的函數，求通過函數上任何一點的斜率：\(\mathrm{lim}_{\mathit{h}\rightarrow 0}\frac{\mathit{f(\xi+h)-f(\xi)}}{\mathit{h}}\)。
把 \(\xi+h\)以及\(\xi\)的差距代入到\(y=x^2-x\)函數，

\[\begin{align*} f(\xi+h)-f(\xi) & =(\xi+h)^2-(\xi+h)-(\xi^2-\xi) \\ & = \xi^2+2\xi\cdot h + h^2-\xi-h-\xi^2+\xi \\ & = 2\xi\cdot h + h^2 -h \tag{6.11} \end{align*}\]

再代入斜率的公式：

\[\begin{align*} m &=\frac{f(\xi+h)-f(\xi)}{\xi+h-\xi} \\ & = \frac{2\xi\cdot h + h^2 -h}{h} \\ & = 2\xi + h -1 \tag{6.12} \end{align*}\]

當\(h\)趨近於0，方程式 (6.12) 將接近於 \(2\xi-1\)，也就是一條斜率為2的直線。我們把\(m\)改為\(f^{`}\)，並且定義\(f^{`}\)：

\[f^{`}=\frac{\partial y}{\partial x}\]

當\(y=x^{a}\)，以下關係成立：

\[f^{`}(x)=ax^{a-1}\]

圖6.3中的藍色虛線代表切線，定義為通過\((\xi,f(\xi))\)這個點而且斜率為\(m\)的直線。
綠色的點則是當\(f{`}=0\)時的切線，此時\(f(\xi)\)是\(f(x)\)的最低點。心算可知當\(\xi\)等於0.5時，\(f^{`}=0\)。

因此\(\hat{\beta}_{1}\)也代表當X變動一個單位時Y變動程度，因為對以下方程式的X微分，表示該函數的斜率：

\[\begin{align*} Y & =\hat{\beta}_{0}+ \hat{\beta}_{1}X+\hat{u} \\ \frac{\partial Y(\hat{\beta}_{0}, \hat{\beta}_{1})}{\partial X} & =\hat{\beta}_{1} \end{align*}\]

結論，對一個函數微分，會得到特定點通過該函數的切線的斜率，也就是變化的速度。以\(f(x)=x^2-x\)函數而言，\(f^{\prime}(X)=2x-1\)，因此當X=1時，斜率等於1，以此類推。

6.3.1 範例

★ Granovetter and Soong (1988) 提出一個函數說明白人對於社區裡黑人居住的容忍程度的函數如下： \[ f(x)=R\left[1-\frac{x}{N_{w}}\right]x \]

其中\(x\)是白人的比例，\(R\)是白人對黑人的忍受度，\(N_{w}\)是黑人人數，\(f(x)\)代表黑人人數與白人搬走的關係。
如果\(N_{w}=100\)，\(R=5\)，請求出當白人的比例x是20%時，白人移出的比例變化是多少？
代入\(N_{w}=100\)，\(R=5\)，得到：

\[ f(x)=5x-\frac{x^2}{20} \]

使用斜率的公式\(\mathrm{lim}_{\mathit{h}\rightarrow 0}\frac{\mathit{f(\xi+h)-f(\xi)}}{\mathit{h}}\)：

\[\begin{align*} f'(x) = \mathrm{lim}_{\mathit{h}\rightarrow 0}\frac{[5(x+h)-\frac{1}{20}(x+h)^2]-[5x-\frac{1}{20}x^2]}{h} \\ & = \mathrm{lim}_{\mathit{h}\rightarrow 0}\frac{5h-\frac{1}{20}(x^2+2xh+h^2)+\frac{1}{20}x^2}{h}\\ & = \mathrm{lim}_{\mathit{h}\rightarrow 0}\frac{5h-\frac{2}{20}xh-\frac{1}{20}h^2}{h}\\ & = \mathrm{lim}_{\mathit{h}\rightarrow 0}\left(5-\frac{1}{10}x-\frac{1}{20}h \right)\\ & = 5-\frac{1}{10}x \end{align*}\]

代入x=20，可以得到\(f'(20)=3\)。如果代入x=40，可以得到\(f'(40)=1\)。如果代入x=50，可以得到\(f'(50)=0\)。這個結果顯示，不同的X會對應到不同的切線有不同的斜率。
\(f(x)=5x-\frac{x^2}{20}\)的函數圖形如圖6.4：

# Create a scatter plot
plot(1, type = 'n', xlim = c(-6, 80), ylim = c(0, 130), xlab = 'X', ylab = 'Y')

f <- function(x)(5*x - (1/20)*x^2)
curve(f, -10, 80, col = 'red', ylab = expression(italic(f(x))),
      add = T)
grid(NULL, NULL, col = 'gray50')

# Define the coordinates of the point through which the line passes
x_point <- 20
y_point <- 5*x_point - (1/20)*x_point^2  

# Add the point to the plot
points(x_point, y_point, col = 'red', pch = 16)

# Add a line with a slope of 3 passing through the point (20, y)

abline(a = y_point - 3 * x_point, b = 3, col = 'blue')

Figure 6.4: 一元二次函數與切線

當白人的比例x是20%時，白人移出的比例變化是3。

6.4 三項最小平方法迴歸的衍伸

最小平方法迴歸所衍伸的幾項特性：

自變項與殘差並無相關。表示為：

\[\sum {x}_{i}\hat{u}_{i}=0\]

證明： \[\begin{align*} \sum x_{i}\hat{u}_{i} & = \sum x_{i}(y_{i}-\hat{\beta}_{0}-\hat{\beta}_{1}x_{i})\\ & = \sum x_{i}y_{i}-\hat{\beta}_{0}\sum x_{i}-\hat{\beta}_{1}\sum x_{i}^2 \end{align*}\]
根據normal equations，也就是式(6.5)，\(\hat{\beta}_{0}\sum x_{i}-\hat{\beta}_{1}\sum x_{i}^2=\sum x_{i}y_{i}\)

迴歸線必定通過y的平均值，也會通過x的平均值。寫成：

\[\bar{\hat{y}}=\bar{y}\]

證明： \[\begin{align*} \bar{\hat{y}} & =\frac{1}{n}\sum \hat{y}\\ & = \frac{1}{n}\sum(\bar{y}-\beta_{1}\bar{x}+\beta_{1}x_{i})\\ & = \frac{1}{n}n\bar{y}+\frac{\beta_{1}}{n}\sum (x_{i}-\bar{x}) \\ & = \bar{y}+\frac{\beta_{1}}{n}(\sum x_{i}-n\frac{\sum x_{i}}{n}) \\ & = \bar{y} \end{align*}\]

殘差與預測值沒有相關

\[\sum \hat{y}_{i}\hat{u}_{i}=0\]

\[\begin{align*} \sum \hat{y}_{i}\hat{u}_{i} & = \sum (\hat{\beta}_{0}+\hat{\beta}_{1}x_{i})\hat{u}_{i}\\ & = \sum \hat{\beta}_{0}\hat{u}_{i}+\hat{\beta}_{1}\sum x_{i}\hat{u}_{1} \\ &= \hat{\beta}_{0}\sum \hat{u}_{i}+\hat{\beta}_{1}\sum x_{i}\hat{u}_{1} \end{align*}\]

因為：

\[\sum \hat{u}_{i}= 0\qquad\sum x_{i}\hat{u}_{1}=0\]

所以：

\[\sum \hat{\beta}_{0}\hat{u}_{i}+\hat{\beta}_{1}\sum x_{i}\hat{u}_{1}= 0\]

image3<-here::here('Fig','regybar0.jpg')
knitr::include_graphics(image3)

Figure 6.5: 有Y平均值的迴歸線圖

image5<-here::here('Fig','regybarxbar.jpg')
knitr::include_graphics(image5)

Figure 6.6: 有X與Y平均值的迴歸線圖

image6<-here::here('Fig','regybar1.jpg')
knitr::include_graphics(image6)

Figure 6.7: 有Y平均值以及迴歸線的迴歸圖線

7 迴歸係數的意義

要知道迴歸的參數\(\beta_{0}\)以及\(\beta_{1}\)，我們必須估計\(\hat{\beta_{0}}\)跟\(\hat{\beta_{1}}\)，方法是盡量降低誤差的總和。我們用之前提到的RSS(Residual Squared Sum)表示

\[\begin{align*} u &= y_{i}-\hat{y}_{i} \\ &= y_{i}- \tilde{\beta}_{0}-\tilde{\beta}_{1}x_{i} \end{align*}\]

OLS迴歸的目標是最小化殘差，殘差可以計算為\(y_{i}-\hat{y}_{i}\)，但是這會給所有的\({y}_{i}\)一樣的權重。改用殘差的平方和 (\(\it{SSR}\)), 才會給離迴歸線越遠的\({y}_{i}\)比較大的權重，這樣才能反映離差的大小。

\[S(\tilde{u})=\mathcal{\sum\limits_{i=1}^n(y_{i}-\tilde{u})^2}\]

OLS迴歸的目的是找到最小化以下模型的\(\{\hat{\beta}_{0},\hat{\beta}_{1}\}\)

\[ \{\tilde{\beta}_{0}, \tilde{\beta}_{1}\}= \rm{minimize} (\sum _{i=1}^n(y_{i}-\tilde{\beta}_{0}-x_{i}\tilde{\beta}_{1})^2)\]

樣本平均數就是最佳的殘差，也就是最小平方法的樣本統計(estimator)

\[\hat{u} \equiv \tilde{u} = \frac{1}{n}\mathcal {\sum\limits_{i=1}^n y_{i}} \]

經過對上述方程式微分之後，迴歸係數可計算為

\[ \mathcal{\hat{\beta}_{1}=\frac{\sum y_{i}-\hat{\beta}_{0}n}{\sum x_{i}}} \] \[ \mathcal{\hat{\beta}_{0}=\frac{\sum x_{i}y_{i}-\hat{\beta}_{1}\sum x_{i}^2}{\sum x_{i}}} \]

經過推導，\(\hat{\beta_{1}}\)可改計算為：

\[ \hat{\beta_{1}}=\mathcal{n\sum (x_{i}-\bar{x})(y_{i}-\bar{y})} \]

7.1 從線性函數推導

假設\(Y=f(X)\)

當X從\(a_{1}\)改變到\(a_{2}\)，\(f(a_{1})\)也變成\(f(a_{2})\)

因此，X的變動程度造成Y的變動程度表示為\(\mathcal{m=\frac{\mathit{f(a_{2})-f(a_{1})}}{\mathit{a_{2}-a_{1}}}}\)，也就是線性函數的斜率。也有人表示為\(\frac{\Delta Y}{\Delta X}\)

其中\(a_{2}>a_{1}\)

如果用微分方式，對於可微分的函數，求通過函數上任何一點的斜率：\(\mathcal{\mathrm{lim}_{\mathit{h}\rightarrow 0}\frac{\mathit{f(x+h)-f(x)}}{\mathit{h}}}\)。

因此，微分函數之後代表X變動一個單位，Y的變動

另一種迴歸係數的計算方式：

\[\begin{align*} \hat{\beta}_{1} & = \frac{\sum_{i=1}^n(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum_{i=1}^n(x_{i}-\bar{x})^2} \\ & = \frac{Sample\hspace{.5em} covariance\hspace{.5em} of\hspace{.5em} X\hspace{.5em} and\hspace{.5em} Y}{sample\hspace{.5em} variance\hspace{.5em} of\hspace{.5em} X} \end{align*}\]

另一方面， \[ \hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1}\bar{x} \]

根據以上的公式可以得出兩個結論：
X的變異數越大，\(\hat{\beta}_{1}\)越小
X, Y的共變量越大，\(\hat{\beta}_{1}\)越大

7.2 各項名詞的意義

\(y_{i}-\bar{y}\)等於是用y的平均值無法預測到\(y_{i}\)的部分

\(\hat{u}_{i}=y_{i}-\hat{y}_{i}\)是用迴歸線無法預測到\(y_{i}\)的部分

\(\hat{y}-\bar{y}\)是y的預測值與平均值的差異。

\(\sum _{i}^n(y_{i}-\bar{y})^2\)=SST=Var[y]

也稱為Total Sum of Squares

y的變異量(variability)可表示如下：

Residual Sum of Squares

\[\begin{align*} \sum _{i}^n(\hat{y}_{i}-y_{i})^2 & = \sum _{i}^n\hat{u}^2\\ & = SSR \\ & = \rm{Var}[\hat{u}] \end{align*}\]

y的預測值與平均值的差異則是：

Explained Sum of Squares

\[\begin{align*} \sum _{i}^n(\hat{y}-\bar{y})^2 & = SSE\\ & = Var[\hat{y}] \end{align*}\]

這三者的關係寫成：SST=SSE+SSR

8 \(R^2\)的意義

因為SST=SSE+SSR，所以有以下的關係：

\[ \mathcal{\frac{SST}{SST}=\frac{SSE}{SST}+\frac{SSR}{SST}} \] 以及
\[ \mathcal{\Longleftrightarrow \frac{SSE}{SST}=1-\frac{SSR}{SST}\equiv R^2} \] \[R^2=\frac{\sum (\hat{y}-\bar{y})^2}{\sum (y-\bar{y})^2} \]

\(R^{2}\)可以詮釋成Y與X相關，\(0\leq R^{2} \leq 1\)。

x <- seq(1,10,length.out = 100)
y <- 2 + 1.2*x + rnorm(100,0,sd = 2)
cat("Y~X, R-squared:", summary(lm(y ~ x))$r.squared, "\n")
## Y~X, R-squared: 0.762
cat("X~Y, R-squared:", summary(lm(x ~ y))$r.squared, "\n")
## X~Y, R-squared: 0.762

如果\(R^{2}=0\), \(\hat{y}=\bar{y}\)，X完全無法解釋Y

如果\(R^{2}=1\), 每一個點\((x_{i},y_{i})\)落在迴歸線上

不同的資料有可能得到相近的\(R^2\)

8.1 \(R^2\)的限制

\(R^2\)指的是「目前」的模型解釋依變數的程度，我們當然可以改良模型以提高解釋的程度，但是前提是我們正確地測量依變數以及自變數。如果我們沒辦法正確地測量，即使\(R^2\)很高，也可能是偶然發生。

圖8.1的左邊圖顯示，很多點都聚集在迴歸線上，右因此\(R^2\)高達0.89。左圖顯示觀察值比較分散，\(R^2\)達0.68。但是\(\beta_{1}\)幾乎一樣，分別為1.6與1.7。因此，。

set.seed(116)
obs = rnorm(50, 10, 1)
dm<-data.frame(obs=rep(obs,2),
          Y = c(rep('y1',50), rep('y2', 50)),     
          value = c(1 + 1.5*obs + runif(50, 0, 3),
                1 + 1.5*obs + rchisq(50, 2)))
ggpubr::ggscatter(dm, x='obs', y='value',add='reg.line',
                  color='Y', palette = "jco") +
     facet_wrap(~ Y) +
    stat_cor(label.y = 30) +
  stat_regline_equation(label.y = 27)

Figure 8.1: 不同的\(R^2\)迴歸線

圖8.2則顯示改變起始的亂數得到相似但是不一樣的迴歸估計結果。也就是說，我們永遠不確定目前的模型是不是是真實的模型，因為每一次抽樣得到的樣本可能不同，所以不需要太強調模型適合度。

set.seed(21)
obs = rnorm(50, 10, 1)
dm<-data.frame(obs=rep(obs,2),
          Y = c(rep('y1',50), rep('y2', 50)),     
          value = c(1 + 1.5*obs + runif(50, 0, 3),
                1 + 1.5*obs + rchisq(50, 2)))
ggpubr::ggscatter(dm, x='obs', y='value',add='reg.line',
                  color='Y', palette = "jco") +
     facet_wrap(~ Y) +
    stat_cor(label.y = 30) +
  stat_regline_equation(label.y = 27)

Figure 8.2: 不同的亂數的迴歸線

不過，如果\(R^2\)很小，我們就不能用這個模型產生正確的預測值。

8.2 \(R^2\)與自變數

當依變數的變異程度變大，而其他變數不變，\(R^2\)有可能會變小。我們首先自訂一個參數為Y的變異數的函數。

set.seed(3939889)
r2.0 <- function(sig){
  x <- seq(1,10, length.out = 300)        # our predictor
  y <- 2 + 1.2*x + rnorm(300, 0, sd = sig) # our response; a function of x plus some random noise
  summary(lm(y ~ x))$r.squared   # print the R-squared value
}

然後產生20個連續的Y的變異數，以常態分佈模擬Y之後，放入迴歸模型模擬得到20個\(R^2\)，圖8.3，顯示Y的變異程度越大，\(R^2\)也越小。

sigmas <- seq(1, 10, length.out = 20)
rout <- sapply(sigmas, r2.0)  # apply our function to a series of sigma values
dt <- data.frame(sigmas, rout)

library(ggplot2)
ggplot(data=dt, aes(x = sigmas, y = rout)) +
   geom_line(col='#FF6600', size=1.5)  +
    geom_point() +
    labs(y=expression(R^2), x=expression(sigma^2)) +
    theme_classic()

Figure 8.3: 變異數與解釋變異量圖1

由上面的圖8.3可以看出，隨著Y的變異數\(\sigma\)變大，\(R^2\)由大變小，但是如果檢視迴歸係數，大概在1.5上下，見圖8.4：

set.seed(3939889)
beta <- function(sig){
  x <- seq(1,10,length.out = 300)        # our predictor
  y <- 2 + 1.2*x + rnorm(300,0, sd = sig) # our response; a function of x plus some random noise
  summary(lm(y ~ x))$coefficients[1,2]           # print the R-squared value
}

sigmas <- seq(1, 10,length.out = 20)
res <- sapply(sigmas, beta)     
dt <- data.frame(x=sigmas, y=res)
library(ggplot2)
ggplot(data=dt, aes(sigmas, res)) +
   geom_line(col='#FF6600', size=1.5)  +
    geom_point() +
  labs(y=expression(beta[1]), x=expression(sigma^2)) +
    theme_classic()

Figure 8.4: 變異數與解釋變異量圖2

因此，我們真正關心的是迴歸係數，而不是\(R^2\)，因為後者有可能隨著\(\sigma^2\)變大而趨近於0，但是迴歸係數穩定增加。這說明\(R^2\)並不代表「模型適合度」。

為什麼依變數的變異程度增加，而其他變數不變，\(R^2\)會降低？因為\(Y\)的變異數\(\sum (Y-\bar{Y})^2\)變大的同時，依變數被解釋的程度有可能下降，也就是SSE(\(\sum(\hat{Y}-\bar{Y}\))變小，所以SSR(\(\sum(Y-\hat{Y}\))會變大，同時\(R^2\)變小。借用上面的語法為例，見圖8.5：

showtext::showtext_auto()

set.seed(02138)

SSRSST <- function(sig){
  x <- seq(1,10,length.out = 300)         # our predictor
  y <- 2 + 1.2*x + rnorm(300, 0, sd = sig) # our response; a function of x plus some random noise
  m1 <- lm(y ~ x)
  return(sum((m1$fitted.values-mean(y))^2))
}

sigmas <- seq(1, 10, length.out = 20)
tmp<- sapply(sigmas, SSRSST)    # apply our function return SSR

plot(sigmas, tmp, main='SSE與Y的變異數散佈圖',
      ylab = expression(SSE),
     xlab = expression(sigma^2), col='#ee1122')

Figure 8.5: SST與變異數

當\(Y\)的變異數增加，理論上SST變大，\(\frac{SSE}{SST}\)變小，造成\(R^2\)變小。

所以\(Y\)的離散程度不必然會提高被\(X\)解釋的程度。

8.3 調整後\(R^2\)

調整後的\(R^2\)考慮樣本數\(n\)以及自變數的數目\(k\)： \[ R^2_{Adj}=1-\frac{(1-R^2)(n-1)}{n-k-1} \]

8.4 幾種資料的\(R^2\)

以下幾種資料的迴歸模型的\(R^2\)很大，但是散佈圖顯示X與Y的線性關係不一定成立。

第一種是多數觀察值集中在左下角，少數集中在右上方，如圖8.6：

DT<-data.frame(Y=c(1,1,1,1,1,2,2,2,2,2,3,3,3,10),
X=c(1,2,1,1,2,1,2,1,3,1,2,3,1,11))
m1 <- lm(Y ~ X, data=DT)
library(ggpubr)
ggscatter(DT, x = "X", y = "Y", add = "reg.line", color = '#aecc02') +
  stat_cor(label.x = 3, label.y = 32) +
  stat_regline_equation(label.x = 3, label.y = 35)

Figure 8.6: R2散佈圖1

第二種是觀察值呈現非線性關係，如圖8.7：

library(tidyverse)
DT<-tibble(
X=seq(1, 10, by=0.5),
Y=X^2)
library(ggpubr)
ggscatter(DT, x = "X", y = "Y", add = "reg.line", color = '#c12c02') +
  stat_cor(label.x = 4, label.y = 75) +
  stat_regline_equation(label.x = 4, label.y = 85)

Figure 8.7: R2散佈圖2

第三種是多數的觀察值呈現線性關係，但是許多觀察值集中在右邊，例如圖8.8：

DTnew<-tibble(X=c(7, 7.5, 8, 8.9, 8.7, 1, 2), 
              Y=c(19, 19.5, 20, 30 ,  23,   1, 1.3))
library(ggpubr)
ggscatter(DTnew, x = "X", y = "Y", add = "reg.line", color = '#21baec') +
  stat_cor(label.x = 1, label.y = 32) +
  stat_regline_equation(label.x = 1, label.y = 35)

Figure 8.8: R2散佈圖3

由圖8.6, 8.7, 8.8可以看出，\(R^2\)只能詮釋為自變數跟依變數之間的相關，不能視為模型適合資料的程度，也不能做為判斷模型的標準。

9 作業

1. 某一家民調公司進行手機加上市話民調，其中市話完成740通，手機完成460通。市話樣本接觸到40歲以下民眾與40歲以上民眾分別是45%以及55%，手機樣本接觸到40歲以下民眾與40歲以上民眾分別是75%以及25%。在被問到接受過幾年的學校教育時，得到以下的表格：

	20-39	40-90
Landline	13	11
Cellular	17	14

請問根據以上的資料，受訪者的平均受教育時間多長？

2. 使用DSS的資料country.csv，有170個國家的GDP(gdp, prior_gdp)的變數以及光源(light, prior_light)的變數，請問光源的變動比例與GDP的變動比例有正相關嗎？迴歸係數是多少？

3. 根據以上的模型，請問如果未來三年光源變動率是62%、57%以及52%，請預測GDP變動率是多少？（提示：用predict這個語法）

4. 如果模擬兩個常態分佈的變數如下，請問迴歸模型的係數是多少？

Y <- rnorm(1000, 0, 1)
X <- rnorm(1000, 0, 1)

5. 如果用迴圈重複100次迴歸模型估計，請問100個z值（\(\beta_{1}\)的係數除以標準誤）裏面，大於2有幾個？

6. 請針對mtcars這筆資料的hp以及mpg繪製散佈圖與迴歸線圖，並且計算迴歸係數。

7. 請使用UsingR的babies這筆資料，估計年齡(age)與身高(ht)對於嬰兒重量(wt)的影響。注意這些變數有一定的範圍。

8.一公司欲知各經銷站每月所賣出的雜誌數(Y)與業務員人數(X1)、廣告費(X2)（10萬元）之間的關係，觀察得下列之資料，該公司想估計X1影響Y的作用，以及X2影響Y的作用。請寫出最小平方法的估計結果。

X1<-c(1,2,3,2,5,5,2,4,3,3)
X2<-c(4,6,7,6,7,7,5,7,6,5)
Y<-c(30,38,44,38,50,52,33,50,45,40)
mag <- rbind(X1,X2,Y)
knitr::kable(mag, format = 'simple')

X1	1	2	3	2	5	5	2	4	3	3
X2	4	6	7	6	7	7	5	7	6	5
Y	30	38	44	38	50	52	33	50	45	40

9. 承上一題，請用R估計X1影響Y的作用，以及X2影響Y的作用。

10. 請問在nym.2002資料中，性別與年齡會影響完成馬拉松的時間嗎？

接續上一題，請加入性別與年齡的交互作用項，並且討論結果。

10 更新日期

## 最後更新日期 04/25/2025

社會科學統計方法

迴歸原理

蔡佳泓

5/1/2025