社會科學統計方法

迴歸原理

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可看出，年齡越大，完成比賽所需的時間越長。

2 迴歸的意義

迴歸的目的是找到依變數或是反應變數(outcome)Y與一個或是一個以上的預測變數(predictor)X的函數關係。

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

母體迴歸函數(Population regressio function)一般表示為\(E(Y|X)=\beta_{0}+X\beta_{1}+\epsilon\)，或者\(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\)。

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

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

1. 預測：假設誤差項 \(\epsilon\) 的平均值為0。\(\hat{Y}=\hat{f}(X)\)，也就是\(f(X)\)可以預測\(Y\)，因為\(\hat{f}\)只是\(f\)的預測，所以兩者之間有誤差。即使\(\hat{f}\)等於\(f\)，\(\hat{Y}=f(X)+\epsilon\)，因此讓\(E(\epsilon^2)=Var(\epsilon)\)越小，預測就越準確。不過，我們無法減少\(\epsilon\)，因此預測的正確性(accuracy)有上限，實際上無法得知上限有多少。
2. 推論：我們關心的是\(X\)與\(Y\)之間的關係，\(Y\)如何隨著\(X\)而變動。例如，某家公司收集了許多消費者資料，他們想要知道那些背景的消費者對於他們的產品有興趣，也就是從許多變數中找到與\(Y\)相關的變數。此外，\(X\)與\(Y\)有可能是非線性的關係。線性模型的推論比較容易詮釋，但是不一定能預測所有觀察值。非線性模型的預測能力比較好，但是不容易詮釋。
3. 預測與推論兩者並重：一方面估計哪些房屋的特徵會影響房價，一方面預測特定房屋的價格被低估或者高估。

迴歸名詞解釋

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

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

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

2.1 有母數與無母數迴歸

有母數模式(parametric)意為函數是由一組母數所規範；具有某種固定的型態，例如線性、對數等等。無母數(non-parametric)的模式則為未知的函數型態，簡單地講，也就是沒有母數的限制。

母數迴歸模式的最大優點是可以推論迴歸結果，但是其缺點就是迴歸形式固定而呆板。無母數迴歸的優點是符合資料程度高，但是難以推論。

無母數迴歸不假設變數之間有線性關係，目的是選擇合適的平滑參數，讓通過樣本的迴歸線能夠最好地逼近真實的迴歸曲線，也就是讓均方誤差為最小。

無母數迴歸可以作為測試最小平方法迴歸模型是否線性的基礎。

無母數迴歸及散佈圖可以作為瞭解資料分佈的基礎，有可能是非線性分佈。

無母數迴歸係數無法像線性迴歸係數可以進行\(t\)檢定或是推論母體的參數。它只能讓我們觀察最接近每一個點的線。

無母數迴歸使用的樣本統計(estimator)有很大的彈性。帶寬越小，越有彈性，預測值的變異數也越大。

帶寬越大，越不具彈性，與真實模型之間的誤差越大。

2.2 無母數模型

• 核密度估計(Kernel Density Estimation, kde)裡面的\(K\)是一種單峰(unimodal)、對稱的函數，像是標準常態分佈，也稱為高斯(Gaussian)，是核函數的一種，也就是：

\[ \frac{1}{\sqrt{2\pi}}\text{exp}(-\frac{x^2}{2})\]

• 已知核函數，核密度可估計如下：

\[\begin{align*} \hat{f}(x;h)=\frac{1}{nh}\sum_{i=1}^{n}K\frac{x-X_{i}}{h} \end{align*}\]

• K是核函數，類似權重，\(h\)用來控制密度分布的平滑程度，先設定為1。\(n\)則是觀察值的數目。
• 細看這個公式，先假設\(h=1\)，可以想像我們把特定的數與每一個觀察值相減，然後代入核函數，加總之後，再除以全部樣本數\(n\)，就會得到特定的數的核密度。\(h\)越大，得到的核密度分佈越平滑。
• 如何決定\(h\)的大小？這一篇Andrey Akinshin寫的文章提到\(h\approx 1.06\times \hat{\sigma}n^{-1/5}\)適用於核函數為高斯分佈。kedd這個套件專門顯示不同的核函數以及參數的圖形，這是相關的說明。
• 我們用\(\textbf{UsingR::too.young}\) 這筆資料中的男性約會年齡為例，可以畫圖 2.1如下，請注意我們指定bw=5，這是一個特定的帶寬：

lattice::densityplot(UsingR::too.young$Male, bw=5)

Figure 2.1: 核密度圖1

• 計算約會年齡為20歲的核密度如下：

k=function(x){
  (1/sqrt(2*pi))*exp(-(x^2)/2)
}
xs <- UsingR::too.young$Male
x1 = 20
n = length(xs); h = 1
pdf = sum(k(x1-xs))/(n*h)
print(pdf)
## [1] 0.03701

• 結果是0.037左右，大概就是圖 2.1的結果。圖 2.2顯示bw=6，結果x=20的機率下降到0.03：

lattice::densityplot(UsingR::too.young$Male, bw=6)

Figure 2.2: 核密度圖2

• 如果把h改成1.2，得到的核密度大約是0.03，圖 2.2得到近似的估計：

n = length(xs); h = 1.2
pdf = sum(k(x1-xs))/(n*h)
print(pdf)
## [1] 0.03084

2.2.1 區域平均值(Simple Local Average)

• 估計\(E[Y|X]\)的方法之一是「可移動的區域平均」(moving local average)。假設\([X_{0}-h,X_{0}+h]\)這個區間內有許多y，計算這些y的平均值便是「可移動的區域平均」。
• 其中，\(h\)稱為帶寬(bandwidth)。
• 如果在區間內的每一個y都視為同樣的，有相同的權重，稱為「單一核密度」(Uniform Kernel)。又稱為Nadaraya–Watson kernel regression estimate.
• 根據「可移動的區域平均」公式，依照不同的帶寬，可計算得到不同的區域平均值。 將這些區域平均值連起來就是無母數迴歸線，我們用carData::Chirot資料為例，如圖2.3。

im1<-here::here('Fig','nonpreg1.jpg')
knitr::include_graphics(im1)

Figure 2.3: 無母數迴歸線1

2.2.2 區域多項迴歸(Local weighted polynomial regression)}

• 前一個無母數迴歸中，因為每一區間內的觀察值得到同樣的權重，所以形成線條狀的迴歸線。
• 在資料中我們取若干觀察值來計算每一個觀察值到特定觀察值(focal point)之間的距離，不同的區間大小決定該區間會用到多少比例的觀察值，也決定曲線的平滑程度。區間越大，平滑程度越大。。
• 進一步，我們用多元的函數來決定權重。根據權重可以計算出特定觀察值的模型值。模型可表示為：

\(Y_{i}=\alpha_{0}+\beta_{1}(x_{i}-x_{0})+\cdots +\beta_{p}(x_{i}-x_{0})^{p}+E_{i}\)

• 在最小化誤差的平方之後，把估計點連起來平滑曲線，稱為local weighted polynomial function或是Bivariate local polynomial regression.
• 這個只用鄰近資料點的方法，可以避免某些特別大或是特別小觀察值的影響，如圖2.4。

im2<-here::here('Fig','nonpreg2.jpg')
knitr::include_graphics(im2)

Figure 2.4: 無母數迴歸線2

• 可用以下指令來畫無母數迴歸線：

lwl<-loess(intensity ~ inequality, data=carData::Chirot)
plot(intensity ~ inequality, data=carData::Chirot)
j <- order(carData::Chirot$inequality)
lines(carData::Chirot$inequality[j], lwl$fitted[j], col='red', lwd=2, lty=2)

2.2.3 合併

• 合併線性迴歸與無母數迴歸在一張圖，可以看出資料點的散佈並非線性。
• 有許多診斷反應變數與解釋變數是否有線性關係的方式，例如檢視預測值與殘差的散佈圖，如果散佈圖呈現特定形狀，有可能兩者之間非線性關係，如圖2.5。

im3<-here::here('Fig','nonpreg3.jpg')
knitr::include_graphics(im3)

Figure 2.5: 無母數迴歸線3

從以上可知：

• 從點估計的均方誤差和可以想像，變異數與誤差之間有衝突，一個太大，另一個就會太小。

\(MSE_{pred.-y_{i}}=Var_{pred.}+Bias_{\mathit{f_{i}}-y_{i}}^2\)

• 無母數迴歸使用的樣本統計(estimator)有很大的彈性。
• 帶寬越小，越有彈性，預測值的變異數也越大。帶寬越大，越不具彈性，與真實模型之間的誤差越大。

im4<-here::here('Fig','nonpa_4.jpg')
knitr::include_graphics(im4)

Figure 2.6: 無母數迴歸線4

im5<-here::here('Fig','nonpa_8.jpg')
knitr::include_graphics(im5)

Figure 2.7: 無母數迴歸線5

im6<-here::here('Fig','loess_r.jpg')
knitr::include_graphics(im6)

Figure 2.8: 無母數迴歸線6

---

2.2.4 無母數迴歸的重要性

• 無母數迴歸可以作為測試最小平方法迴歸模型是否線性的基礎。
• 無母數迴歸及散佈圖可以作為瞭解資料分佈的基礎，有可能是非線性分佈。
• 無母數迴歸係數無法像線性迴歸係數，可以進行\(t\)檢定或是推論母體的參數。它只能讓我們觀察最接近每一個點的線。
• 無母數統計不容易詮釋，因此我們以有母數統計的線性迴歸推論變數之間的關係。

3 有母數統計：線性迴歸

線性迴歸假設條件平均值呈現線性。

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

只有兩個係數。\(\beta_{0}\)稱為截距或常數，\(\beta_{1}\)稱為斜率係數。當X=0，\(E[Y|X=0]\)等於Y的期望值，也就是\(\beta_{0}\)

線性假設的其中之一是截距為固定。 使用觀察資料，或者訓練資料估計模型之後，可以得到模型的預測值，表示為： \(\hat{y}=\hat{\beta}_{0}+\hat{\beta}_{1}x\)

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

定義殘差的平方和為：

\(RSS=e_{1}^2+e_{2}^2+\ldots +e_{n}^2\)

3.1 什麼是E(Y|X)？

E(Y|X)表示固定X時，Y的期望值。E(Y|X)稱為conditional expectation function。例如，\(E(Y|X=1)\neq E(Y|X=2)\)

當X是隨機變數，E(Y|X=x)也是隨機變數，E(Y|X)隨著X變化而改變，可以寫成h(Y)。

因此，E(Y|X=x)就像Y一樣，有期望值也有變異數。X的變異數為：

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

• 或者

\[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)}\]

• 例如（見Rachel Fewster）：

\[\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} \end{equation}\]

• 而且，

\[\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} \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} \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.1.3 Law of Iterated Expectations

條件期待函數有一個特性是條件期待值的期待值剛好等於該自變數的期待值。這是因為對於所以自變數的值，都會有一個發生的機率，所以：

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

• 可以推導出：

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

稱為Law of Iterated Expectations(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

3.2 最小平方法

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

image1<-here::here('Fig','reg1.jpg')
knitr::include_graphics(image1)

Figure 3.1: 迴歸線圖

3.3 迴歸係數的意義

迴歸模型\(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})\)

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

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

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

• 用 \(\tt{lm()}\)估計迴歸方程式如3.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 3.1: 電視廣告與銷售量之間的關係

	Dependent variable:
	sales
TV	0.048***
	(0.003)
Constant	7.033***
	(0.458)
Observations	200
R2	0.612
Adjusted R2	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與迴歸線預測的之間會有差異。
• 畫圖3.2如下：

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 3.2: 電視廣告與銷售量之間的關係

• 加上預測值以及觀察值與預測值之間的殘差（參考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()

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

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

假設有\(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{4.1} \end{align}\]

• 我們稱\(L(y)\)為loss function。

4.1 推導過程

對(4.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。把這些平均值連起來就是迴歸線。

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

• 圖4.1顯示，當\(\tilde{u}=y-\hat{y}=\bar{y}\)=E(Y|X)時，y與其他觀察值之間的誤差平方和會最小。

set.seed(666)
Y <- rnorm(1000, mean = 10, sd = 2)
mean(Y)
## [1] 9.961
obs <- rnorm(50, 10, 3)

# generate empty list
squared_residuals <- rep(NA, length(obs))

for (i in 1:length(obs)) {
  squared_residuals[i] = sum((Y - obs[i])^2)
}

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=mean(Y) , 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 4.1: 誤差平方和函數圖

4.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*}\]

4.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{4.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{4.3} \end{align*}\]

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

4.1.3 normal equations

• 根據(4.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{4.4} \end{align*}\]

• 根據(4.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{4.5} \end{align*}\]

• 這兩個等式稱為normal equations

• 從normal equations的第1式((4.4)可以得到：

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

• 根據式(4.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{4.7} \end{align*}\]

• 把方程式(4.7)代入(4.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{4.8} \end{align*}\]

• 在(4.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{4.9} \end{align*}\]

• 在(4.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{4.10} \end{align*}\] （因為\(\sum \bar{x}^2=\bar{x}^2+\bar{x}^2+\cdots +\bar{x}^2=n\bar{x}^2\)）

4.2 OLS與共變數

根據(4.9)以及(4.10)，式(4.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-{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\)把以上的數據代入(4.8)，

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

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

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

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

• 圖4.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\)來表示。

p <- poly_from_zeros((1):0)
p
## -x + x^2
plot(p, xlab = expression(italic(x)), ylab = expression(italic(P(x))),
  main = parse(text = paste("italic(P(x) ==",
                             as.character(p, decreasing = TRUE),")")),
  xlim=c(-1,4))
x0 <- solve(deriv(p))        ## stationary points
lines(tangent(p, x0), col = "dark green", lty = "solid",
      limits = cbind(x0-1/4, x0+1/4))
points(x0, p(x0), col = "dark green")

lines(tangent(p, 2), lty=2)
x<-c(2,3,3) 
y<-c(2,2,6)
polygon(x, y, lty=2, lwd=2)
segments(-1.5, 2, 2, 2, lty=2, col='red')
segments(-1.5, 6, 3, 6, lty=2, col='red')

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

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

p <- poly_from_zeros((1):0)
plot(p, xlab = expression(italic(x)), ylab = expression(italic(f(x))),
  main = parse(text = paste("italic(f(x) ==",
              as.character(p, decreasing = TRUE),")")),
  xlim=c(-1,4), xaxt='n', yaxt='n')
axis(2, at=c(2, 6), c(expression(paste(f(xi))),expression(paste(f(xi+h)))))
axis(1, at=c(2, 3), c(expression(paste(xi)),expression(paste(xi+h))))
x0 <- solve(deriv(p))        ## stationary points
lines(tangent(p, x0), col = "dark green", lty = "solid",
      limits = cbind(x0-1/4, x0+1/4))
points(x0, p(x0), col = "dark green")

lines(tangent(p, 2), lty=2, col='blue', lwd=1.6)
x<-c(2,3,3) 
y<-c(2,2,6)
polygon(x, y, lty=2, lwd=2)
segments(-1.5, 2, 2, 2, lty=2, col='red')
segments(-1.5, 6, 3, 6, lty=2, col='red')

Figure 4.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{4.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{4.12} \end{align*}\]

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

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

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

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

• 圖4.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*}\]

4.4 四項最小平方法迴歸的衍伸

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

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

\[\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，也就是式(4.5)，\(\hat{\beta}_{0}\sum x_{i}-\hat{\beta}_{1}\sum x_{i}^2=\sum x_{i}y_{i}\)

2. 迴歸線必定通過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*}\]

3. 殘差與預測值沒有相關

\[\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 4.4: 有Y平均值的迴歸線圖

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

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

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

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

---

5 迴歸係數的意義

要知道迴歸的參數\(\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})} \]

5.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}\)越大

5.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

6 \(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.6985
cat("X~Y, R-squared:", summary(lm(x ~ y))$r.squared, "\n")
## X~Y, R-squared: 0.6985

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

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

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

6.1 \(R^2\)的限制

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

圖6.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 6.1: 不同的\(R^2\)迴歸線

• 圖6.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 6.2: 不同的亂數的迴歸線

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

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

當依變數的變異程度變大，而其他變數不變，\(R^2\)會變小（為什麼？），但是\(X\)與\(Y\)之間的關係並沒有改變。我們首先自訂一個參數為Y的變異數的函數。

set.seed(3939889)
r2.0 <- function(sig){
  x <- seq(1,10,length.out = 100)        # our predictor
  y <- 2 + 1.2*x + rnorm(100,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\)，圖6.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 6.3: 變異數與解釋變異量圖1

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

set.seed(3939889)
beta <- function(sig){
  x <- seq(1,10,length.out = 100)        # our predictor
  y <- 2 + 1.2*x + rnorm(100,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)
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 6.4: 變異數與解釋變異量圖2

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

為什麼依變數的變異程度變大，而其他變數不變，\(R^2\)會變小？因為\(Y\)的變異數變大的同時，SSE雖然變大，但是SST也變大。借用上面的語法舉其中兩個\(\sigma\)為例，見表6.1：

set.seed(5525252)
Y.pre <- function(sig){
  x <- seq(1,10,length.out = 100)        # our predictor
  y <- 2 + 1.2*x + rnorm(100,0, sd = sig) # our response; a function of x plus some random noise
  m1<-lm(y ~ x)
  return(m1$fitted.values)          # fitted values
}
Y.mean <- function(sig){
  x <- seq(1,10,length.out = 100)
  y <- 2 + 1.2*x + rnorm(100,0, sd = sig) 
  return(mean(y))          # fitted values
}
Y.values <- function(sig){
  x <- seq(1,10,length.out = 100)
  y <- 2 + 1.2*x + rnorm(100,0, sd = sig) # our response
  return(y)
}
sigmas <- c(1, 10)
tmp1<- sapply(sigmas, Y.pre)    # apply our function 
tmp2<- sapply(sigmas, Y.mean)
tmp3<- sapply(sigmas, Y.values)
library(dplyr); library(data.table)
tmp <- data.table(Y.pre1=tmp1[,1], Y.mean1=tmp2[1], Y.values1=tmp3[,1], 
    Y.pre2=tmp1[,2], Y.mean2=tmp2[2], Y.values2=tmp3[,2])
  
tmp %>% summarise(SSR1=sum(Y.pre1-Y.mean1)^2, SST1=sum(Y.values1-Y.mean1)^2, 
      R.squared1=SSR1/SST1, 
    SSR2=sum(Y.pre2-Y.mean2)^2, SST2=sum(Y.values2-Y.mean2)^2,
    R.squared2=SSR2/SST2) 
##    SSR1 SST1 R.squared1  SSR2 SST2 R.squared2
## 1 208.3 11.6      17.95 42504 5031      8.448
kable(head(tmp, n=10), caption='SST與變異數') %>%
  kable_styling(full_width = F)

Table 6.1: SST與變異數

Y.pre1	Y.mean1	Y.values1	Y.pre2	Y.mean2	Y.values2
3.375	8.619	2.796	4.158	7.655	-2.276
3.483	8.619	2.727	4.270	7.655	3.049
3.592	8.619	2.268	4.382	7.655	17.000
3.701	8.619	3.778	4.495	7.655	17.140
3.810	8.619	3.584	4.607	7.655	4.024
3.919	8.619	3.414	4.719	7.655	-10.736
4.028	8.619	4.375	4.832	7.655	5.456
4.137	8.619	1.992	4.944	7.655	-1.539
4.246	8.619	5.026	5.056	7.655	12.123
4.354	8.619	6.581	5.168	7.655	-10.050

因此當\(Y\)的變異數增加，雖然預測值的離散程度因此增加，也就是預測值不會集中在某一個點，但是佔所有變異量的比例可能下降，\(R^2\)變小。

6.3 調整後\(R^2\)

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

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

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

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

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 6.5: R2散佈圖1

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

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 6.6: R2散佈圖2

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

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 6.7: R2散佈圖3

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

---

7 作業

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

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

3.一公司欲知各經銷站每月所賣出的雜誌數(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

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

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

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

8 更新日期

## 最後更新日期 05/13/2022