社會科學統計方法

迴歸原理

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的平均值。所以模型的對象是\(E[Y|X]\)，也就是X變化時，Y的平均值變化大小。或者是說固定X情況下，Y的平均值大小。

母體迴歸函數(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\)的方法。統計學習的目的是：

1. 預測：假設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)有上限，實際上無法得知上限有多少。

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)用來估計這個未知的機率畢度函數(PDF)。
• 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這個套件專門顯示不同的核函數以及參數的圖形，這是相關的說明。

• 假設有三個觀察值3, 4, 7，他們來自什麼樣的分佈？用KDE估計如下：

 # Define data and bandwidth
x <- c(3, 4, 7)
h <- 1

# Define a sequence of points where we want to evaluate the density
x_grid <- seq(min(x) - 2*h, max(x) + 2*h, length=100)

# Function to calculate the normal kernel density for a single point
kernel <- function(x_i, mu, h) {
  dnorm(x_i, mu, h)  # dnorm is the R function for normal density function
}

# Calculate kernel density for each data point
densities <- sapply(x, function(xi) kernel(x_grid, xi, h))

# Sum the kernel densities from each data point
density_estimate <- apply(densities, 1, sum)

# Normalize by the number of data points
density_estimate <- density_estimate / length(x)

# Plot the kernel density estimate
plot(x_grid, density_estimate, type = "l", 
     main = "Kernel Density Estimation (Normal Kernel)",
     xlab = "Value", ylab = "Density")

- 我們用\(\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

2.2.4 誤差與變異數的關係

• 從無母數迴歸與線性迴歸的比較中，發現無母數迴歸線的點比線性迴歸線的點更集中在線旁邊，表示誤差比較小，但是無母數迴歸線上面的點之間的離散程度比較高，也就是比較有彈性但是也比較有不確定性。

• 用計算均方誤差和的方式，來觀察變異數與誤差之間有衝突，一個太大，另一個就會太小。

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

• 定義mean squared error (MSE)=\(E[(y - f(x))^2]\)，其中

• \(y=f(x)+\epsilon\)

• \(y\)是觀察值

• \(f(x)\)是真實的模型

• \(\hat{f}(x)\)是樣本的模型

• \(\epsilon\)為y與\(f(x)\)的誤差

• 計算MSE如下：

\[ E[(y - \hat{f}(x))^2] = E[y^2] - 2E[y\times \hat{f}(x)] + E[\hat{f}(x)^2] \]

• 經過整理之後，並且假設\(\epsilon\)的期待值等於0，得到：

\[\begin{align*} E[\epsilon^2] + \Big(E[f(x)] - E[\hat{f}(x))]\Big)^2 + Var(f(x)-\hat{f}(x)) \\ & = \text{irreducible error}+\text{bias}+\text{variance of prediction} \end{align*}\]

• 在同樣的MSE之下，當\(\Big(E[f(x)] - E[\hat{f}(x))]\Big)^2\)越大，樣本模型的預測平均值與真實預測值之間的差距越大，會造成預測值之間的變異程度（\(Var(f(x)-\hat{f}(x))\)）變小，也就是模型受到特定觀察值而改變預測值的機會越小（越沒有彈性）。

• 無母數迴歸使用的樣本統計(estimator)有很大的彈性。帶寬越小，越有彈性，預測值的變異數也越大。帶寬越大，越不具彈性，與真實模型之間的誤差越大。例如圖2.7與圖2.8顯示，不同自變數X的區間有不同的帶寬，區間越多，誤差越小，但是變異數越大。

• 原始資料的散佈圖：

with(carData::Chirot, plot(intensity ~ inequality))

Figure 2.6: 散佈圖

• 劃分X為4個區間，每個區間取一個平均數：

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

Figure 2.7: 無母數迴歸線4

• 劃分X為8個區間，每個區間取一個平均數：

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

Figure 2.8: 無母數迴歸線5

• 圖2.9同時比較線性迴歸線與非線性迴歸線：

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

Figure 2.9: 無母數迴歸線6

---

2.2.5 無母數迴歸的重要性

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

3 有母數統計：線性迴歸

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

\[ 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}\)。 - 因為平均數是一個常數，常數的期望值等於該常數，也就是：

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

• 當自變數減去自己平均值後，依變數的期望值等於截距。也就是說，當我們還沒有用最小平方法估計\(\beta_{1}\)，我們至少知道依變數在與Ｘ無關的情況下的平均值等於\(\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=x時，Y的期望值（平均值）。\(E[Y|X]\)稱為conditional expectation function(CEF)。例如，\(E(Y|X=1)=\)_{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.1)：

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

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

\[\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.2} \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.3} \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.3)，計算\(\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.4} \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

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

\[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

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

\[\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.5} \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