1 課程目標

本週上課將介紹二元迴歸的基本原理，描述自變數與二元變數之間的關係。例如：

2 二元迴歸的基本原理

假設每一個\(Y\)代表從參數為\(P_{j}\)抽樣得到的伯努利實驗變數，\(Y_{j}=\sum^{N_{j}} Y_{ij}\)，也就是在\(j\)的群體中，累積\(Y=1\)的次數\(N_{j}\)等於\(Y_{j}\)。

什麼是伯努利實驗？只有兩種結果的單次實驗。例如：擲出硬幣得到正面或者反面？參加考試是否及格？明天是否會下雨？期望值\(E[Y|X]=p\)，標準差為\(\sqrt{p(1-p)}。\)

而從機器學習的觀點，如何根據訓練資料找到\(f(X)\)，未來可以預測行為是贊成或反對、結果是成功或是失敗，是估計二元迴歸模型的目的。

累積\(Y=1\)的次數\(N_{j}\)等於\(Y_{j}\)的模型表示為：

\[ \begin{eqnarray} \frac{Y_{j}}{N_{j}} & = & f_{j}\\ & = &\sum \beta_{k}X_{ik}+u_{i} \end{eqnarray} \]

假設\(j\)代表投給國民黨或是民進黨。如果有兩個以上，例如\(j=3\)。我們需要估計\(j-1\)個方程式。

以上的方程式\(\frac{Y_{j}}{N_{j}}\)等於機率\(P\)，而已知\(0\leq P\leq 1\)。根據最小平方法得到的\(\hat{Y}\)有可能大於1或是小於0，也就是得到錯誤的預測值。

此外，\(X\), \(Y\)之間可能不是線性關係；當\(X\)變動一單位，\(Y\)變動的單位可能因為\(X\)在不相同的值而不相等。因此，\(p\)必須經過轉換，才不會出現超出上下限的預測值，也才會有線性關係。

2.0.1 實例

\(\blacksquare\) 什麼樣的人會無法償還信用卡的債務？

library(ISLR); library(stargazer)
Default$default.n <- as.numeric(Default$default)
Default$default.n <- Default$default.n-1 #No=0, Yes=1
M1 <- lm(default.n ~ balance, data=Default)
Default$predicted<-predict(M1)
stargazer (Default, type='html')

根據上面的估計結果，可得出以下的表格，其中\(\hat{y}_{i}\)的極小值小於0。


Statistic	N	Mean	St. Dev.	Min	Max

balance	10,000	835.375	483.715	0.000	2,654.323
income	10,000	33,516.980	13,336.640	771.968	73,554.230
default.n	10,000	0.033	0.179	0	1
predicted	10,000	0.033	0.063	-0.075	0.270

左邊的圖表示部分的\(\hat{y}_{i}\)小於0，右邊的圖表示經過轉換，\(\hat{y}{i}\)介於0與1之間。

此外，用OLS估計二元的依變數，其模型誤差\(\epsilon\)的變異數會隨著自變數的值而改變，也就不是固定不動，造成係數的標準誤不是最小。

最後，OLS模型如果估計二元的依變數誤差也不是呈現常態分佈。

2.1 連結函數：依變數的轉換

為了轉換機率的依變數，我們介紹對數。

對數log的特性

\[exp^{log(a)} = a\]

指數exponential的特性

\[\mathrm{log}(exp^{a})=a \] \[exp^{a-b}=\frac{exp^{a}}{exp^{b}} \] \[exp^{a+b}=exp^{a}exp^{b} \] 我們實際計算看看：

a=100
exp(log(a))

## [1] 100

log(exp(a))

## [1] 100

進一步計算指數：

a=1
exp(a)

## [1] 2.718282

exp(-a)

## [1] 0.3678794

假設\(Y'=log(Y)\)，\(log(Y)\equiv Y'\equiv X\beta^{*}+\epsilon\)

當\(X\)變動一個單位，\(Y\)變動 \(exp^{\beta^{*}}\)單位，約等於\(\beta^{*}\%\)

此處\(log(Y)\)稱為Link Function，寫成\(F(\cdot)=F(Y)=Y'=log(Y) = X\beta+\epsilon\)

3 Logit 模型的基本原理

logistic function:

\[ \begin{eqnarray} \text{logit}(Y=1|X) & = & \text{log} \frac{P_{i}}{1-P_{i}}\\ & = & \frac{p(X)}{1-p(X)}\\ &=&\beta_{0}+X\beta_{1} \end{eqnarray} \] 因為\(\text{exp}(\text{log}(a))=a\)，所以 \[ \frac{p(X)}{1-p(X)} = exp^{\eta} \] 其中， \[ \eta \equiv \sum \beta_{k}X_{ik}\]

為何logit做為一個連結函數，會產生\(\{-\infty, \infty \}\)的對應值？這是因為logit的函數是一個勝算比的型態。

何謂勝算？可表示成\(\frac{p(x)}{1-p(x)}\)。例如：

\(\blacksquare\) \(20\%\)的打擊率，表示打出安打跟無法上壘的勝算為1:4，因為\(\frac{0.2}{1-0.2}=\frac{1}{4}\)，或者說無法上壘的機率是打安打的4倍。
\(\blacksquare\) \(90\%\)的合格率，表示合格是不合格機率的9倍，因為\(\frac{0.9}{1-0.9}=9\)。

3.0.1 實例

\(\blacksquare\) 什麼樣的人會無法償還信用卡的債務？

library(ISLR); library(stargazer)
M2 <- glm(default ~ balance, data=Default, family=binomial(logit))
stargazer (M2, type='html')


	Dependent variable:

	default

balance	0.005^***
	(0.0002)

Constant	-10.651^***
	(0.361)


Observations	10,000
Log Likelihood	-798.226
Akaike Inf. Crit.	1,600.452

Note: ^***:p<0.01

由以上可知，\(\hat{\beta}_{1}=0.0055\)，代表增加一個單位的balance，會增加default的指數勝算(log odds)為0.0055，或者是勝算增加\(\text{exp}^{0.0055}=1.005515\)。

勝算比(odds ratio)是兩個勝算的比值，例如某家連鎖超市舉辦抽獎活動，經過統計，有八成的參加者是男性，兩成是女性，也就是\(p=0.8\)，\(q=1-p=0.2\)。分別計算勝算值：

\[ \text{odds(mobile)}=\frac{0.8}{0.2}=4 \] \[ \text{odds(land-line)}=\frac{0.2}{0.8}=0.25 \]

根據上述，男性參加抽獎是女性的16倍，因為\(\frac{4}{0.25}=16\)

那麼logit模型跟勝算比有什麼關係？因為經過轉換之後，

\[ \begin{eqnarray} \beta_{1}&=&\text{log(odds}_{x=1})-\text{log(odds}_{x=0})\\ &=&\text{log}(\frac{\text{odds}_{x=1}}{\text {odds}_{x=0}}) \end{eqnarray} \]

所以\(\text{exp}^{\beta_{1}}\)等於當=1是x=0的y=1的倍數。

M3 <- glm(default ~ student, data=Default, family=binomial(logit))
stargazer (M3, type='html')


	Dependent variable:

	default

studentYes	0.405^***
	(0.115)

Constant	-3.504^***
	(0.071)


Observations	10,000
Log Likelihood	-1,454.342
Akaike Inf. Crit.	2,912.683

Note:	p<0.05; p<0.01; p<0.001

\(\hat{\beta}_{1}=0.404\)，代表學生比非學生有\(\text{exp}^{0.404}=1.497\)倍的勝算發生違約。

3.1 實例

假設有一個自變數的值為1到5，估計logit模型得到ln\(\frac{P_{i}}{1-P_{i}}=1+0.5X\)

當x=1，exp(x\(\beta\))=exp(1.5)=4.48，而1+exp(3.5)=5.48，Pr(y=1|x=1)=\(\frac{4.48}{5.48}\)=0.81

當x=2，Pr(y=1|x=2)=0.88

當x=3，Pr(y=1|x=3)=0.92

可以看出當x增加1個單位，機率增加的程度會因為x從哪一個值增加而有所不同。

但是exp(\(\beta\))= exp(0.5)=1.64，表示增加勝算1.64-1=64%。增加勝算的幅度是固定的，並不會因為x的大小而不同。

而exp(\(\beta\))= exp(0.5)=1.64，表示x每增加1個單位，y是1對上y是0的勝算比，或者是勝算的對數（log of odds）增加1.64。

3.2 比較OLS與Logit模型估計結果

library(ISLR); library(stargazer)
Default$default.n <- as.numeric(Default$default)
Default$default.n <- Default$default.n-1 #No=0, Yes=1
M1 <- lm(default.n ~ balance, data=Default)
logM <- glm(default.n ~ balance, data=Default,
            family=binomial(logit))
stargazer(M1, logM, type="html")


	Dependent variable:

	default.n
	OLS	logistic
	(1)	(2)

balance	0.0001^***	0.005^***
	(0.00000)	(0.0002)

Constant	-0.075^***	-10.651^***
	(0.003)	(0.361)


Observations	10,000	10,000
R²	0.123
Adjusted R²	0.122
Log Likelihood		-798.226
Akaike Inf. Crit.		1,600.452
Residual Std. Error	0.168 (df = 9998)
F Statistic	1,396.816^*** (df = 1; 9998)

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

4 Probit 模型的基本原理

由於依變數為二元變數\(Y_{i}=0\)或1，服從伯努利分佈，參數為次數\(N\)與機率\(P\)，表示成\(P_{i}(Y_{i}=1)\)

為了轉換依變數\(P_{i}\)從\(\{0, 1\}\)變成\(\{-\infty, \infty \}\)，我們考慮常態分佈曲線又稱probit，中文為二元機率單元模型，又稱為成長曲線迴歸（黃紀、王德育，2012:87）。

probit function: \[ \text{Pr(y=1|x)}=F(Z)=\int_{-\infty}^{\beta_{0}+X\beta_{1}}\frac{1}{2\pi}exp(-u^2/2)du\equiv \Phi (Z)\]

probit function 轉換\(Z\)或者是\(X\beta\)成為常態分佈的\(Z\)分數，所以\(X\beta\)越大，Pr(Y=1)越大。

Logit, Probit兩種模型得到的結果類似，但是前者的係數比較容易詮釋。

4.1 Probit 模型

常態分布\(N(0, 1)\)的累積機率密度可計算如下：

當x\(\beta\)=-2，Pr(y=1)=\(\Phi(-2)=0.022\)

當x\(\beta\)=-1，Pr(y=1)=\(\Phi(-1)=0.158\)

當x\(\beta\)=2，Pr(y=1)=\(\Phi(2)=0.977\)

可用R計算：

pnorm(-2, 0, 1)

## [1] 0.02275013

pnorm(-1, 0, 1)

## [1] 0.1586553

pnorm(2, 0, 1)

## [1] 0.9772499

計算全部X\(\beta\)從-3到3在常態分佈時\(N(0,1)\)的累積機率密度：

j<-c()
k<-seq(-3, 3, 0.1)
n <- length(k)

for (i in 1:n)
{
  j[i] <- pnorm (k[i], 0, 1)
  
}

plot(k, j, main="CDF", xlab=expression(paste('X',beta)), ylab="Probability",
     type="l", lwd=1.5)

由以上可知，\(0<\Phi(Z)<1\)，而且累積的機率形成S形的曲線。

當x增加一個單位，Pr(y=1) 增加\(\beta\)的Z值。

Probit的累積機率(cdf)可畫成如下：

4.2 Probit模型的詮釋

參考UCLA的介紹：

Admit<-read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
Admit$rank <- factor(Admit$rank)
fit1 <- glm(admit ~ gre + gpa + rank, family = binomial(link = "probit"), 
    data = Admit)
stargazer(fit1, type='html')


	Dependent variable:

	admit

gre	0.001^**
	(0.001)

gpa	0.478^**
	(0.197)

rank2	-0.415^**
	(0.195)

rank3	-0.812^***
	(0.208)

rank4	-0.936^***
	(0.245)

Constant	-2.387^***
	(0.674)


Observations	400
Log Likelihood	-229.207
Akaike Inf. Crit.	470.413

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

Gre增加一個單位，z分數增加0.001

畫圖顯示：

newdata <- data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100), 
    4 * 4), gpa = rep(c(2.5, 3, 3.5, 4), each = 100 * 4), rank = factor(rep(rep(1:4, 
    each = 100), 4)))
newdata[, c("p", "se")] <- predict(fit1, newdata, 
                                   type = "response", se.fit = TRUE)[-3]
library(ggplot2)
ggplot(newdata, aes(x = gre, y = p, colour = rank)) + geom_line() + facet_wrap(~gpa)

5 最大概似法

最大概似法(Maximum Likelihood Estimation, MLE)的原理是對每一組\((x_{i}, y_{i})\)，嘗試Pr(\(y_{i}=1\))=\(\Phi(\text{X}\beta)\)之中的\(\beta\)值，最大化出現Pr(\(y_{i}=1\))的機率。

例如我們嘗試\(\beta\)值得到Pr(\(y_{i}=1\))=0.8，而且實際上\(y_{i}=1\)，我們可以說得到的概似機率為0.8。如果實際上\(y_{i}=0\)，我們可以說得到的概似機率為0.2。

我們用\(\mathcal{L}(y_{i}|\beta)\)代表有\(\beta\)的條件下的\(y_{i}\)的概似機率，如果連乘所有的機率\(\mathcal{L}(y_{1})\mathcal{L}(y_{2})\cdot\ldots\cdot\mathcal{L}(y_{n})=\prod_{i=1}^n\mathcal{L}(y_{i})\)

\[ \begin{eqnarray} \text{max}\sum_{i=1}^{n}[\text{log}\mathcal{L}(y_{i}|\hat{\pi})]&=&\sum \text{log}[\hat{\pi}^{y_{i}}(1-\hat{\pi})^{1-y_{i}}] \nonumber \\ & = & \sum y_{i}\text{log}(\hat{\pi})+(1-y_{i})\text{log}(1-\hat{\pi}) \end{eqnarray} \]

對方程式的\(\hat{\pi}\)取微分，並且設為0求出\(\hat{\pi}\)

\[\pi=\frac{\sum_{i=1}^n y_{i}}{n} \]

5.1 MLE

MLE的過程，需要使用Talyor series，讓概似機率逼近真實的機率。

概似機率的估計式的抽樣分佈會接近常態分佈，N(0, H)，其中H代表Hessian matrix （海森矩陣）

\(\hat{\pi}=\frac{\sum_{i=1}^n y_{i}}{n}\)表示樣本平均數可以最大化概似函數。

將\(\hat{\pi}\)代入\(\pi=\frac{1}{1+exp{-X\beta}}\)，應用重複加權(iteratively reweighted least squares (IRLS)）的演算式，聚合以下算式：

\[\hat{\beta}=\mathcal{(X^{T}WX)^{-1}X^{T}Wz} \]

進一步計算\(\beta\)的標準誤，表示為：

\[\hat{\beta}\sim \mathcal{N(\beta, \phi(X^{T}WX)^{-1})} \]

5.2 MLE的應用

概似機率的估計式的抽樣分佈會接近常態分佈，N(0, H)，其中H代表Hessian matrix （海森矩陣）

MLE也適用在線性迴歸，見Gary King (1989) 。最小平方法是MLE的特例。

貝氏統計也需要MLE構成可觀察資料的模型，乘以先驗的資訊得到後驗的機率分佈。

6 信賴區間

當依變數\(y_{i}\hspace{.5em}\sim\hspace{.5em} \text{Bin}(1,\pi_{i})\)，logit 模型可寫成：

\[\text{ln}(\frac{\pi_{i}}{1-\pi_{i}})=X\beta \]

或者是 logistic 迴歸模型：

\[ \pi_{i}=\frac{exp(X\beta)}{1+exp(X\beta)} \]

應用重複加權(IRLS)的演算式，計算\(\widehat{SE}\)，得到上下區間：

\[ (\hat{\pi}-z_{\alpha/2}\widehat{SE},\hspace{1em}\hat{\pi}+z_{\alpha/2}\widehat{SE}) \] 而且 \[\frac{\hat{\beta_{j}} - \beta_{j}}{\widehat{SE}}\sim z \]

其中，\(\widehat{SE}=\sqrt{x^{T}(X{T}WX)^{-1}x)}\)

運用IRLS 可以得到與R一樣的估計

特別注意要去掉遺漏值

見Patrick Breheny的語法

6.1 假設檢定

假設檢定是要測量係數的正確程度。

檢定虛無假設\(\text{H}_{0}\):\(\pi=\pi_{0}\)，根據這項假設計算\(x\beta\)，也就是\(\hat{\beta}_{1}/\textit{SE}(\hat{\beta}_{1})\)

如果\(\hat{\beta}_{1}=0\)，等於\(p(X)=\frac{e^{\beta_{0}}}{1+e^{\beta_{0}}}\)表示機率與X無關。但是我們不關心\(\beta_{0}\)。

但是實際上我們直接計算區間估計，判斷是否顯著不等於0。

區間估計分為Wald 與 Likelihood ratio兩種

6.1.1 Wald test

Abraham Wald 依據MLE發展出Wald confidence intervals, Wald test statistics等等

Wald 的方法依賴對於概似機率的趨近，如果概似機率不佳，也就是\(\hat{\beta}\)與\(\beta\)相差很大，Wald的檢定結果就會不夠好

6.1.2 LR test

LR 檢定的原理是比較兩個模型，一個是具有所有資訊(以\(\theta\)表示)的模型，一個是具有部分資訊(以\(\hat{\theta}\)表示)的模型。這兩個模型的比例表示為

\[\lambda=\frac{\mathcal{L}(\theta)}{\mathcal{L}(\hat{\theta})} \]

當\(n\longrightarrow \infty\)，

\[-2\hspace{.1em}\text{log}\hspace{.1em}\lambda\hspace{.5em}\xrightarrow{d}\hspace{.5em}\chi^2\]

根據以上的定理，可導出

\[-2\text{log}\frac{L(\hat{\beta}|\beta^{1}=\beta^{1}_{0})}{L(\hat{\beta})}\leq \chi^2_{1-\alpha, q} \]

其中\(q\)代表自由度，\(\chi^2_{1-\alpha, q}\)代表\(\chi^2\)分布的(1-\(\alpha\))百分位

7 多變數二元迴歸模型

我們可以延伸上述的logit迴歸模型到多變數：

\[ \text{log}(\frac{p(X)}{1-p(X)})=\beta_{0}+\beta_{1}X_{1}+\ldots+\beta_{p}X_{p} \]

因為多於一個變數，所以在詮釋其中一個變數的影響時，需要強調其他的變數在同樣水準。

在計算預測值時，必須設定變數特定的值，除了我們所關心的自變數。例如平均值、眾數等等。

7.1 實例

\(\blacksquare\)哪些因素影響信用卡違約？

library(ISLR); library(stargazer)
M5 <-glm(default ~ balance+income+student, data=Default, family=binomial(logit))
stargazer(M5, type='html')


	Dependent variable:

	default

balance	0.006^***
	(0.0002)

income	0.00000
	(0.00001)

studentYes	-0.647^***
	(0.236)

Constant	-10.869^***
	(0.492)


Observations	10,000
Log Likelihood	-785.772
Akaike Inf. Crit.	1,579.545

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

同樣收入的學生比非學生來得不會違約。

7.1.1 exp

以下指令可建立\(\text{exp}^{\beta}\)：

M5.or <- exp(coef(M5))
stargazer(M5, coef=list(M5.or), type='html')


	Dependent variable:

	default

balance	1.006^***
	(0.0002)

income	1.000^***
	(0.00001)

studentYes	0.524^**
	(0.236)

Constant	0.00002
	(0.492)


Observations	10,000
Log Likelihood	-785.772
Akaike Inf. Crit.	1,579.545

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

7.1.2 預測值

先建立一個資料框，包含三個自變數的平均值、最大值或者眾數。

allcoef <-data.frame(balance=rep(mean(Default$balance),2),
                     income=rep(mean(Default$income),2),
                     student=as.factor(c("Yes", "No")))
allcoef

##    balance   income student
## 1 835.3749 33516.98     Yes
## 2 835.3749 33516.98      No

代入模型M5

allcoef$pred <-predict(M5, newdata=allcoef, type="response")
allcoef<-cbind(allcoef, allcoef$pred)
allcoef

##    balance   income student        pred allcoef$pred
## 1 835.3749 33516.98     Yes 0.001328976  0.001328976
## 2 835.3749 33516.98      No 0.002534451  0.002534451

7.1.3 Logit, Probit

library(ISLR); library(stargazer)
M5 <-glm(default ~ balance+income+student, data=Default, family=binomial(logit))
M6 <-glm(default ~ balance+income+student, data=Default, family=binomial(probit))
stargazer(M5, M6, type='html')


	Dependent variable:

	default
	logistic	probit
	(1)	(2)

balance	0.006^***	0.003^***
	(0.0002)	(0.0001)

income	0.00000	0.00000
	(0.00001)	(0.00000)

studentYes	-0.647^***	-0.296^**
	(0.236)	(0.119)

Constant	-10.869^***	-5.475^***
	(0.492)	(0.238)


Observations	10,000	10,000
Log Likelihood	-785.772	-791.609
Akaike Inf. Crit.	1,579.545	1,591.217

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

8 作業

1. 請用以下的資料估計性別與錄取的關係：

library(kableExtra)
DT <-data.frame(admission=c(rep(1,10),rep(0,10)),
                gender=c(rep(1,7), rep(0,3),rep(1,3),rep(0,7)))
DT %>% 
  kable('html') %>% 
  kable_styling()

admission	gender
1	1
1	1
1	1
1	1
1	1
1	1
1	1
1	0
1	0
1	0
0	1
0	1
0	1
0	0
0	0
0	0
0	0
0	0
0	0
0	0

2. 請估計以下資料中，暴露於某種物質(Ag)與是否存活(Surv)的關係

library(kableExtra)
DT <-data.frame(Surv=c(rep(1,11),rep(0,22)),
                Ag=c(rep(1,9), rep(0,2),rep(1,8),rep(0,14)))
DT %>% 
  kable('html') %>% 
  kable_styling()

Surv	Ag
1	1
1	1
1	1
1	1
1	1
1	1
1	1
1	1
1	1
1	0
1	0
0	1
0	1
0	1
0	1
0	1
0	1
0	1
0	1
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0

3. 接續上題，請比較probit與logit的模型估計結果

社會科學統計方法

二元迴歸(v1.5)