1 課程目標
2 二元變數迴歸模型的基本原理
- 2.1 實例
- 2.2 連結函數：依變數的轉換
3 Logit 模型的基本原理
4 最大概似法(MLE)原理
- 4.1 最大概似法求解
- 4.2 常態分佈變數的MLE
5 估計二元勝算對數迴歸模型係數
- 5.1 二元勝算對數迴歸模型實例
- 5.2 勝算比(odds ratio)
  - 5.2.1 實例
  - 5.2.2 雙變數二元勝算對數迴歸模型
6 Probit 模型的原理
7 信賴區間
- 7.1 MLE的應用
- 7.2 假設檢定
  - 7.2.1 Wald test
  - 7.2.2 Liklihood ratio (LR) test
8 多變數二元勝算對數迴歸模型
- 8.1 實例
9 比較不同logit或者Probit模型的係數
10 小結
11 作業
12 更新日期

1 課程目標

線性迴歸適用於依變數是連續變數，自變數可以是連續或者類別變數。以最小平方法計算樣本的迴歸係數，推論母體的迴歸係數。也就是：

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

社會科學的研究對象經常是類別變數，例如有去投票跟沒去投票，民主與非民主國家，有接受外援與沒有接受等等。當依變數是類別變數時，如果使用最小平方法線性迴歸模型估計迴歸係數，有可能違反誤差項的變異數均一的假設，所以我們需要轉換依變數，成為新的線性模型。
本週上課將介紹二元變數迴歸的基本原理，描述自變數與二元依變數之間的關係。例如信用卡是否違約與信用卡未償金額的關係。並且介紹如何用R運算最大概似法函數。

2 二元變數迴歸模型的基本原理

假設每一個隨機變數\(Y\)代表從參數為\(P_{j}\)抽樣得到的伯努利實驗變數，也就是在\(j\)的群體中，累積\(Y=1\)的次數\(N_{j}\)等於\(Y_{j}\)。
伯努利實驗是只有兩種結果的單次實驗。例如：擲出硬幣得到正面或者反面？參加考試是否及格？明天是否會下雨？從網路抽出100張照片，其中有寵物的比例？隨機變數的期望值\(E[Y]=p\)，標準差為\(\sqrt{p(1-p)}\)。
如果重複類似實驗\(n\)次，得到\(k\)成功的結果將形成二項分布。\(X \sim B(n,p)\)。
從機器學習的觀點，如何根據訓練資料找到\(f(X)\)，未來加入新的資料之後可以預測結果，例如病人有沒有罹患癌症、學生有沒有通過測驗等等。
累積\(Y=1\)的次數\(N_{j}\)等於\(Y_{j}\)的模型表示為方程式 (2.1)：

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

假設\(j\)代表投給國民黨或是民進黨。如果有兩個以上類別，例如\(j=3\)。我們需要估計\(j-1\)個方程式。
以上的方程式\(\frac{Y_{j}}{N_{j}}\)等於機率\(p\)，並且已知\(0\leq p\leq 1\)。根據最小平方法得到的\(\hat{Y}\)有可能大於1或是小於0，也就是預測值超出觀察值\(p\)的範圍。
此外，\(X\), \(Y\)之間可能不是線性關係；當\(X\)變動一單位，\(Y\)變動的單位可能因為\(X\)在不相同的值而不相等。因此，\(p\)必須經過轉換，才不會出現超出上下限的預測值，也才會有線性關係。

2.1 實例

\(\blacksquare\) 什麼樣的人會無法償還信用卡的債務？是財務不好的人嗎？表 2.1 呈現用最小平方法（OLS）估計的結果：

library(ISLR); data("Default")
Default$default <- as.numeric(ISLR::Default$default)-1
#Default$default <- Default$default #No=0, Yes=1
M1 <- lm(default ~ balance, data=Default)
stargazer::stargazer(M1,  title='信用卡債務與財務',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 2.1: **信用卡債務與財務**

	Dependent variable:

	default

balance	0.0001^***
	(0.00000)

Constant	-0.075^***
	(0.003)


Observations	10,000
R²	0.123
Adjusted R²	0.122
Residual Std. Error	0.168 (df = 9998)
F Statistic	1,397.000^*** (df = 1; 9998)

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

表 2.1 顯示，財務越好的人，也越會發生信用卡債務違約。
圖 2.1以及圖2.2 顯示，模型所產生的部分預測值小於0，也就是小於依變數的最小值，顯然不合理。

Default$predicted<-predict(M1)
plot(seq(1,10000), sort(Default$predicted), type='l', xlab='')

Figure 2.1: 排序後的預測值

r.b<-range(Default$balance)
r.b
## [1]    0 2654
xbalance <- seq(0, 2655, 0.1)
ybalance <- predict(M1, list(balance=xbalance), type='response')
plot(Default$default ~ Default$balance, pch=16,
     col='#33aaeeee', xlab="Balance", ylab="Default")
lines(xbalance, ybalance, col='#eebb00')

Figure 2.2: 自變數與預測值

我們改用\(\tt{glm()}\)估計自變數與依變數關係，畫預測值與自變數的散佈圖2.3：

M2 <- glm(default ~ balance, data=Default, family = 
            binomial('logit'))
xbalance <- seq(0, 2655, 0.1)
ybalance <- predict(M2, list(balance=xbalance), type='response')
plot(Default$default ~ Default$balance, pch=16,
     col='#33aaeeee', xlab="Balance", ylab="Default")
lines(xbalance, ybalance, col='#eebb00')

Figure 2.3: glm模型的自變數與預測值

圖2.3顯示經過轉換，\(\hat{y}{i}\)介於0與1之間，跟原始的\(y\)不同。換句話說，

\[0\leq\rm{Pr(Y=1|X)}=\beta_{0}+\beta_{1}X\leq 1\quad\forall\hspace{.2em}X\]

此外，用最小平方法（OLS）估計二元的依變數，其模型誤差\(\epsilon\)的變異數會隨著自變數的值而改變，也就不是固定不動，造成係數的標準誤不是最穩定。
最後，OLS模型如果估計二元的依變數，其誤差不是呈現常態分佈，違反迴歸假設，檢定的結果可能有誤。

基於以上幾點，我們不使用線性迴歸模型估計自變數與二元變數的關係，改採用Logit或是Probit模型。

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

對數log或ln幫助我們轉換機率的依變數為負無限大到無限大。

log有以下特性：

log(a)的指的是以\(e\)為底對a取自然對數。\(e\)是數學常數，約等於2.718282(\(e\approx 2.718282\))。
log(a)的正式定義是積分，公式為：

\[\rm{log}(a)=\int_{1}^{a}\frac{1}{x}dx\]

\(\bullet\) 對數具有以下特性：

\[\mathrm{log}1=0\] \[\mathrm{log}_{a}a=1\] \[\mathrm{log}(a)+\mathrm{log}(b)=\mathrm{log}(a*b)\] \[\mathrm{log}(a)-\mathrm{log}(b)=\mathrm{log}(a/b)\] \[\mathrm{log}(a) =\infty\hspace{.1cm} \rm{if}\hspace{.1cm} a\approx \infty\] \[\mathrm{log}0=-\mathrm{Inf}\] \[\mathrm{log}a^k=k*\mathrm{log}a\] \[\rm{exp}^{log(a)} = a\] \[\mathrm{log}(\rm{exp}(a))=a\] \[\mathrm{log}(e)= 1\]

我們實際計算看看：

a=10; b=5; k=3
log10(a)
## [1] 1
log(a)+log(b); log(a*b)
## [1] 3.912
## [1] 3.912
log(a)-log(b); log(a/b)
## [1] 0.6931
## [1] 0.6931
log(a^k); k*log(a)
## [1] 6.908
## [1] 6.908
exp(log(a))
## [1] 10
log(exp(b))
## [1] 5

\(\bullet\)我們也會使用指數exponential來轉換對數的函數。指數有以下特性：

\[\rm{exp}^{0}=1\] \[\rm{exp}^{1}=2.718282\] \[\rm{exp}^{a-b}=\frac{\rm{exp}^{a}}{\rm{exp}^{b}} \] \[\rm{exp}^{a+b}=\rm{exp}^{a}\rm{exp}^{b} \]

指數運算的例子：

a=10; b=1
exp(0)
## [1] 1
exp(a-b); exp(a)/exp(b)
## [1] 8103
## [1] 8103
exp(a+b); exp(a)*exp(b)
## [1] 59874
## [1] 59874
log(exp(1))
## [1] 1

假設\(Y=\sum X\beta+\epsilon\)。當我們知道\(Y\)只有0跟1，我們對\(Y\)取對數，寫成\(Y'=\rm{log}(Y)\)，因此\(\rm{log}(Y)\equiv Y'\equiv \sum X\beta^{*}+\epsilon\)，我們估計的是\(\beta^{*}\)而不是\(\beta\) 。
此處\(\rm{log}(Y)\)稱為Link Function，寫成\(F(\cdot)=F(Y)=Y'=\rm{log}(Y) = X\beta+\epsilon\)
但是這種轉換並不如Probit或者Logit，因為\(-\infty \leq \rm{log}(Y)\leq 0\)，我們希望的模型是:

Pr(Y=1|X) 隨著X增加而增加
\(0\leq \rm(Pr(Y=1|X)\leq 1\quad\forall\quad X\)

3 Logit 模型的基本原理

Logit 迴歸模型可寫成:

\[\begin{align*} \rm{logit(Pr(Y=1|X))} & = \rm{log} \frac{P_{i}}{1-P_{i}}\\ &=\beta_{0}+X\beta_{1} \tag{3.1} \end{align*}\]

在方程式 (3.1) 中，\(\frac{P_{i}}{1-P_{i}}\)稱為勝算(odds)，所謂勝算指的是發生某個事件對不發生某個事件的機率。 - 何謂勝算？可表示成\(\frac{p(x)}{1-p(x)}\)。例如：

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

\(\bullet\) \(90\%\)的合格率，表示合格是不合格機率的9倍，因為\(\frac{0.9}{1-0.9}=9\)。

所以我們稱logistic regreesion model為二元勝算對數模型。但是也有人稱為勝算比(odds ratio)(Gould and Hardin)，這是因為勝算就是一種「比」(ratio)的概念，而且，因為\(\text{exp}(\text{log}(a))=a\)，所以

\[\begin{align*} \frac{p_{i}}{1-p_{i}} = exp^{\eta} \tag{3.2} \end{align*}\]

其中，

\[ \eta \equiv \sum \beta_{k}X_{ik}\]

也可以寫成（為什麼？）： \[\begin{align*} p_{i} & =\rm{logit}^{-1}(\eta_{i}) \\ & =\frac{e^\eta_{i}}{1+e^\eta_{i}} \end{align*}\]
在\(\tt{Stata}\)，\(Y\)的機率可以用logistic或者Probit的累積機率函數來表示，以\(F\)做為累積機率的符號如下：

\[\begin{align*} \rm{Pr}(Y=1|X)&=F(\boldsymbol{X'}\beta)\\ &=\frac{1}{1+exp^{-(\boldsymbol{X'}\beta)}} \tag{3.3} \end{align*}\]

這裡的\(\boldsymbol{X'}\)指的是一連串的變數X的矩陣，轉置之後跟\(\beta\)相乘，等於\(\beta_{0}+\beta_{1}X_{1}+\ldots + \beta_{k}X_{k}\)

為何要用logit做為一個連結函數，這是因為logit是一個勝算比的型態，它會產生\(\{-\infty, \infty \}\)的對應值，不受到\(0\le p_{i}\le 1\)的限制。

我們可以用\(\tt{car::logit()}\)這個函式得到logit轉換後的值如圖3.1：

Y<-seq(0, 1, by = 0.01)
Yprime<-car::logit(Y)
plot(Yprime ~ Y, xlab='Y', ylab=expression(paste("Y"^"'")),
     type='l', col='#bb2211')
abline(v=0.5, lty=2); abline(h=0, lty=2)

Figure 3.1: Logit轉換後的值

比較圖3.1以及圖3.2，都是從\(-\infty\)到\(\infty\)的曲線：

plot(log(Y/(1-Y)) ~ Y, type='l', col='#013ae2')

Figure 3.2: 勝算值與Y的散佈圖

如果把\(F(\cdot)\)裡面的累積機率密度函數改為常態分佈，稱為Probit或者機率單元模型。我們可以用\(\tt{qnorm(p)}\)得到轉換\(Y\)之後的結果，如圖3.3:

Y<-seq(0,1, by=0.01)
plot(qnorm(Y) ~ Y, type='l', col='#1111aa')

Figure 3.3: 勝算值與Y的散佈圖

結合圖3.3跟圖3.2，圖3.4顯示兩種累積機率密度函數很接近，所以這兩種模型得到的結果也很相近。

plot(qnorm(Y) ~ Y, type='l', col='#1111aa')
lines(log(Y/(1-Y)) ~ Y, type='l', col='#a0013ae2')
leg.text=c('Probit','Logit')
legend(0, 2.3, leg.text, col=c('#1111aa','#a0013ae2'), lty=c(1,1), bty='n')

Figure 3.4: 勝算值與Y的散佈圖

Logit以及Probit都符合我們希望的模型，也就是:

Pr(Y=1|X) 隨著X增加而增加
\(0\leq \rm Pr(Y=1|X)\leq 1\quad\forall\quad X\)

4 最大概似法(MLE)原理

那麼我們要如何估計這兩種模型的係數？我們可以用電腦透過最大概似法找到最佳的解。

最大概似法(Maximum Likelihood Estimation, MLE)的原理是對每一組\((x_{i}, y_{i})\)，嘗試Pr(\(y_{i}=1\))=\(\Phi(\text{X}\beta)\)之中的\(\beta\)值，最大化出現Pr(\(y_{i}=1\))的機率。參考StatQuest 的說明。
MLE也可以想成找到最佳的方法逼近資料的分佈，例如常態分佈、指數分佈、Gamma分布等等，讓同樣分佈的資料可以重複驗證。
例如我們嘗試\(\beta\)值得到Pr(\(y_{i}=1\))=0.8，而且實際上\(y_{i}=1\)，我們可以說得到的概似機率為0.8。如果實際上\(y_{i}=0\)，我們可以說得到的概似機率為0.2。
例如有一個隨機變數分佈如圖 4.1。假設該變數成常態分佈。首先我們嘗試以下兩個常態分佈分佈：

curve(dnorm(x, 0.35, 0.1), from=-1.5, to=1, col='#3322EE', lwd=2)
curve(dnorm(x, -0.8, 0.2), from=-1.5, to=1, col='#FF22EE', lwd=2, add=T)
x<-c(-1, -0.85, -0.6,-0.50, -0.45, -0.30, -0.25, -0.20, -0.15,
     -0.09,  0.00,  0.10,  0.16,  0.30,  0.50)
y<-rep(0, 15)
points(x,y, pch=16, cex=2, col='gray30')

Figure 4.1: 兩個常態分佈

顯然這兩個分佈的中心點與離散程度與實際資料有點差距。離資料的中心點以及離散程度越遠，表示概似程度(likelihood)越小。
所以當我們說資料的MLE估計時，表示我們極大化觀察到的資料接近特定分布的程度。我們改用不同的常態分佈來接近隨機變數，圖 4.2 顯示點離散程度較小：

curve(dnorm(x, -0.22, 0.38), from=-1.5, to=1, col='#FF3333', 
       lwd=2, xlab='y')
y<-c(-1, -0.85, -0.6,-0.50, -0.45, -0.30, -0.25, -0.20, -0.15, 
     -0.09,  0.00,  0.10,  0.16,  0.30,  0.50)
y2<-rep(0, 15)
points(y, y2, pch=16, cex=2, col='gray30')
abline(v=-0.22, lty=2, lwd=1.5, add=T)

Figure 4.2: 第三個常態分佈

4.1 最大概似法求解

4.1.1 常態分佈的最大概似法

常態分佈的機率密度函數如下：

\[\rm{Pr}(y| \mu, \sigma)=\frac{1}{\sqrt{2 \pi\sigma^2}}e^{-(y-\mu)^2/2\sigma^2}\]

常態分佈的最大概似函數如下：

\[\begin{equation*} \mathcal{L}(\mu, \sigma|y_{i}) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-(y-\mu)^2/2\sigma^2} \end{equation*}\]

在既定的\(\sigma\)以及\(y\)，我們不斷地嘗試各種\(\mu\)，以極大化概似機率。這些\(\mu\)對應的最大概似機率將成為一個常態分佈。當機率極大化時，常態分佈曲線也達到最高，斜率等於0，對應的\(\mu\)，理論上等於資料的平均數。
假設只有一個觀察值20。假設\(\sigma=2\)，代入300個\(\hat{\mu}\)求出最大化函數的值，並畫成圖 4.3：

y<-20; sigma=2; mu<-seq(10, 30, length.out = 300)
L <- c()
n=length(mu)
for (i in 1: n){
  L[i] <- 1/(sqrt(2*pi*sigma^2))*exp(-((y-mu[i])^2)/2*sigma^2) 
}
plot(mu, L, type='l', lwd=2, col='#FF3322', xlab=expression(hat(mu)))
abline(h=max(L), lty=2, lwd=2, col='#998833')

Figure 4.3: 最大概似法求解圖

在求出\(\mu\)之後，再根據\(\mu\)以及\(y\)，我們不斷地嘗試各種\(\hat{\sigma}\)，以極大化概似機率。這些\(\hat{\sigma}\)對應的分佈當機率極大化時，曲線也達到最高，斜率等於0。
如果有\(n\)個\(y\)，那麼概似機率需要連乘： \[\mathcal{L}(\mu, \sigma|y_{1},\ldots,y_{n})=\mathcal{L}(\mu, \sigma|y_{1})\times\ldots,\times\mathcal{L}(\mu, \sigma|y_{n})\]
對兩邊取log，再對\(\mu\)以及\(\sigma\)取微分，

\[\begin{align*} \rm{ln}[\mathcal{L}(\mu, \sigma|y_{1},\ldots,y_{n})]&=\rm{ln}\mathcal{(}\frac{1}{\sqrt{2 \pi\sigma^2}}e^{-(y_{1}-\mu)^2/2\sigma^2}\times\ldots\times\frac{1}{\sqrt{2 \pi\sigma^2}}e^{-(y_{2}-\mu)^2/2\sigma^2}\mathcal{)}\\ & =\rm{ln}\mathcal{(}\frac{1}{\sqrt{2 \pi\sigma^2}}e^{-(y_{1}-\mu)^2/2\sigma^2}\mathcal{)}+\ldots+\rm{ln}\mathcal{(}\frac{1}{\sqrt{2 \pi\sigma^2}}e^{-(y_{n}-\mu)^2/2\sigma^2}\mathcal{)} \tag{4.1} \end{align*}\]

其中

\[\begin{align*} \rm{ln}\mathcal{(}\frac{1}{\sqrt{2 \pi\sigma^2}}e^{-(y_{1}-\mu)^2/2\sigma^2}\mathcal{)}&=\rm{ln}(\frac{1}{\sqrt{2 \pi\sigma^2}})+\rm{ln}(e^{-(y_{1}-\mu)^2/2\sigma^2}) \\ & = \rm{ln}[(2\pi\sigma^2)^{-1/2}]-\frac{(y_{1}-\mu)^2}{2\sigma^2}\rm{ln}(e) \\ & = -\frac{1}{2}\rm{ln}[(2\pi\sigma^2)-\frac{(y_{1}-\mu)^2}{2\sigma^2} \\ & = -\frac{1}{2}\rm{ln}(2\pi)-\frac{1}{2}\rm{ln}(\sigma^2)-\frac{(y_{1}-\mu)^2}{2\sigma^2} \\ &=-\frac{1}{2}\rm{ln}(2\pi)-\rm{ln}(\sigma)-\frac{(y_{1}-\mu)^2}{2\sigma^2} \tag{4.2} \end{align*}\]

\(\bullet\)回想：\(\rm{log}(a^k)=k*log(a)\)以及\(\rm{log(e)=1}\)

計算\(n\)個\(\rm{ln}\mathcal{(}\frac{1}{\sqrt{2 \pi\sigma^2}}e^{-(y_{i}-\mu)^2/2\sigma^2}\mathcal{)}\)的總和，方程式(4.1) 變成：

\[\begin{align*} \rm{ln}[\mathcal{L}(\mu, \sigma|y_{1},\ldots,y_{n})]&=-\frac{n}{2}\rm{ln}(2\pi)-\it{n}\rm{ln}(\sigma)-\frac{(y_{1}-\mu)^2}{2\sigma^2}-\ldots-\frac{(y_{n}-\mu)^2}{2\sigma^2} \tag{4.3} \end{align*}\]

對方程式 (4.3) 的\(\mu\)求微分，獲得\(\mu\)的取log對數的最大概似機率函數的頂點：

\[\begin{align*} \frac{\partial}{\partial\mu}\rm{ln}[\mathcal{L}(\mu,\sigma|y_{1},\ldots,y_{n})]&= 0-0+\frac{(y_{1}-\mu)}{\sigma^2}+\ldots+\frac{(y_{n}-\mu)}{\sigma^2} \nonumber\\ &=\frac{1}{\sigma^2}[(y_{1}+\ldots+y_{n})-n\mu] \tag{4.4} \end{align*}\]

方程式 (4.4) 使用到The Chain Rule使\((y_{1}-\mu)^2\)微分為\(2(y_{1}-\mu)\)
對方程式 (4.4) 的\(\sigma\)求微分，獲得\(\sigma\)取log對數的最大概似機率函數的頂點：

\[\begin{align*} \frac{\partial}{\partial\sigma}\rm{ln}[\mathcal{L}(\mu,\sigma|y_{1},\ldots,y_{n})]&= 0-\frac{n}{\sigma}+\frac{(y_{1}-\mu)}{\sigma^3}+\ldots+\frac{(y_{n}-\mu)}{\sigma^3} \tag{4.5} \end{align*}\]

接下來我們計算當方程式 (4.4) 等於0時的\(\mu\)：

\[\begin{align*} 0&=\frac{1}{\sigma^2}[(y_{1}+\ldots+y_{n})-n\mu] \nonumber\\ &= [(y_{1}+\ldots+y_{n})-n\mu] \tag{4.6} \end{align*}\]

方程式 (4.6) 顯示，當\(\mu=\frac{\sum_{i}^{n}y_{i}}{n}\)時，也就是樣本平均值時，極大化最大概似函數。
接下來我們計算當方程式 (4.5) 等於0時的\(\sigma\)：

\[\begin{align} 0&= -\frac{n}{\sigma}+\frac{(y_{1}-\mu)^2}{\sigma^3}+\ldots+\frac{(y_{n}-\mu)^2}{\sigma^3}\nonumber \\ &= -n+\frac{(y_{1}-\mu)^2}{\sigma^2}+\ldots+\frac{(y_{n}-\mu)^2}{\sigma^2}\nonumber \\ &= -n+\frac{1}{\sigma^2}[(y_{1}-\mu)^2+\ldots+(y_{n}-\mu)^2]\nonumber \\ n&= \frac{1}{\sigma^2}[(y_{1}-\mu)^2+\ldots+(y_{n}-\mu)^2]\nonumber \\ n\sigma^2& = (y_{1}-\mu)^2+\ldots+(y_{n}-\mu)^2\nonumber \\ \sigma^2& = \frac{\sum_{1}^{n} (y_{i}-\mu)^2}{n} \tag{4.7} \end{align}\]

方程式 (4.7) 顯示，當\(\sigma\)等於樣本標準差時，極大化最大概似函數。

假設我們有14個觀察值來自常態分佈：6, 5, 10, 3, 1, 8, 2, 9, 0, 2.5, 3.5, 4.2, 4.5, 4.9。計算平均值為4.5，未校正過的標準差（分母=\(n\)）約為2.82。畫圖表示分佈如圖4.4

obs<-c(6, 5, 10, 3, 1, 8, 2, 9, 0, 2.5, 3.5, 4.2, 4.5, 4.9)
dat <- tibble(x=obs, y=0) 
p <- ggplot(data = dat, aes(x = x, y = y)) + geom_point(color='#CB1D00')
p

Figure 4.4: 14個觀察值

仿照圖 4.1，我們模擬兩個常態分佈來觀察哪一個比較接近觀察值？圖 4.5 顯示，\(\mu=4\)、\(\sigma=2.8\)的分布比較接近。

x<-c(6, 5, 10, 3, 1, 8, 2, 9, 0, 2.5, 3.5, 4.2, 4.5, 4.9)
y<-rep(0, 14)
plot(x, y, ylim=c(0, 0.3), xlim=c(-1,13))
curve(dnorm(x, 4, 2.8), from=-1, to=11, col='#3322EE', lwd=2, add=T)
curve(dnorm(x, 2, 1.5), from=-1, to=11, col='#FF22EE', lwd=2, add=T)

points(x,y, pch=16, cex=2, col='gray30')
leg.text=c(expression(paste(mu,'=4,', ' ',sigma,'=2.8')),
           expression(paste(mu,'=2,', ' ',sigma,'=1.5')))
legend(8, 0.28, leg.text, col=c('#3322EE','#FF22EE'), lty=c(1,1), bty='n')

Figure 4.5: 14個觀察值常態分佈

參考Daniel Linares畫圖表示對14個觀察值的資料，嘗試10個\(\sigma\)以及10個\(\mu\)極大化函數如圖4.6與圖4.7。

library(data.table); library(ggplot2)
set.seed(5000)
obs<-c(6, 5, 10, 3, 1, 8, 2, 9, 0, 2.5, 3.5, 4.2, 4.5, 4.9)
logL<-function(p) sum(log(dnorm(obs,p[1],p[2]))) # p=c(p[1],p[2])=c(mu,sigma)
library('plyr')
dLogL<-expand.grid(museq=seq(2, 7, 0.1),sigmaseq=seq(2,4,0.1))
dLogL<-ddply(dLogL,.(museq, sigmaseq),
             transform,logLike=logL(c(museq,sigmaseq)))
#Graph
ggplot(data=dLogL, aes(x=museq,y=logLike,
                      color=factor(sigmaseq)))+
             geom_line()

Figure 4.6: 極大化函數一

圖4.6顯示，當\(\mu\approx 4.5\)、\(\sigma\approx 3\)時，極大化最大概似函數。
圖4.7顯示，當\(\sigma=2.8\)時，極大化最大概似函數。

DT<-data.frame(dLogL)

is.na(DT) <- sapply(DT, is.infinite)

for (j in 1:ncol(DT)) set(DT, which(is.infinite(DT[[j]])), j, NA)

min_value<-min(DT$logLike)
max_value<-max(DT$logLike)

DT[which(DT$logLike==max_value),]
##     museq sigmaseq logLike
## 534   4.5      2.8  -34.39
dt <- DT %>% group_by(sigmaseq) 
#mean(obs) =4.5
#fix mu=4.5
dt <- dt %>% filter(museq==4.5)

ggplot(data=dt, aes(x=sigmaseq, y=logLike))+
  geom_line(col='#339999', size=1)    +
  geom_hline(aes(yintercept=max(logLike, na.rm=T)),   
               color="red", linetype="dashed", size=1) +
  labs(x=expression(hat(sigma)),y='Log Likelihood')  +
   scale_x_continuous(labels=seq(2,4, by=0.1),
                    breaks=seq(2,4, by=0.1))

Figure 4.7: 極大化函數二


#Negative log-likelihood
negLogL<-function(p) -logL(p)
estPar<-optim(c(20,10),negLogL) 
MLEparameters<-estPar$par
MLEparameters
## [1] 4.543 2.821

\(\sigma^2\)的期待值\(E(\sigma^2)=\sigma^2-\frac{1}{n}\sigma^2\)（參考Daijiang Li）。如果\(n\rightarrow \infty\)，\(E(\sigma^2)=\sigma^2\)

4.1.2 二項分佈的最大概似法

接下來說明二元變數的最大概似法。我們用\(\mathcal{L}(y_{i}|\pi)\)代表有\(\pi\)參數的條件下的\(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{equation*} \mathcal{L}(y_{i}) =\begin{cases} \pi & \text{if}\hspace{.15cm} y_{i}=1 \\ 1-\pi & \text{if} \hspace{.15cm} y_{i}=0 \end{cases} \end{equation*}\] \[\begin{align} \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}) \tag{4.8} \end{align}\]

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

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

概似機率的估計式的抽樣分佈會接近常態分佈，N(0, H)，其中H代表Hessian matrix （海森矩陣）。
\(\hat{\pi}=\frac{\sum_{i=1}^n y_{i}}{n}\)表示樣本平均數可以最大化概似函數。
以10個伯努利變數值為例，最大化概似函數之後畫圖4.8，約在\(p=0.696\)時達到最高點，而這10個變數值的平均值等於0.7（參考White）

set.seed(1000)
p.parameter <- 0.6
sequence <- rbinom(10, 1, p.parameter)

log.likelihood.sum <- function(sequence, p) {
  log.likelihood <- sum(log(p)*(sequence==1)) +  sum(log(1-p)*(sequence==0))
}

possible.p <-seq(0,1, length.out=100)
#jpeg('Likelihood_Concavity.jpg')
library('ggplot2')
g1<-qplot(possible.p,
      sapply(possible.p, function (p) {log.likelihood.sum(sequence, p)}),
      geom = 'line',
      main = 'Likelihood as a Function of p',
      xlab = 'p',
      ylab = 'Likelihood')
g1 + scale_x_continuous(labels=seq(0,1, by=0.1),
                        breaks=seq(0,1, by=0.1))

Figure 4.8: 10個伯努利變數值最大化概似函數

#dev.off()
dt<-data.frame(p=possible.p,
      value=sapply(possible.p, function (p) {log.likelihood.sum(sequence, p)}))
dt<-dt[-c(1,100),]
m.value=max(dt$value)
dt[which(dt$value==m.value),]
##        p  value
## 70 0.697 -6.109

4.1.3 二元分佈變數最大概似法實例

用信用卡是否違約為例，說明R的optimize函數應用在二元變數的方式。首先寫下log函數。

def <- as.numeric(ISLR::Default$default)
y<-c()
y[def==2]<-1
y[def==1]<-0
N<-1
logL <- function(p) sum(log(dbinom(y, N, p)))

然後用\(\tt{optimize()}\) 函數求極大值：

logL(0.1)
## [1] -1785
optimize(logL, lower=0, upper=1, maximum=TRUE)
## $maximum
## [1] 0.0333
## 
## $objective
## [1] -1460

計算平均值驗證一樣得到0.0333：

mean(y)
## [1] 0.0333

4.2 常態分佈變數的MLE

根據式 (4.3)，常態分佈的log-likelihood函數可寫成：

\[-\frac{n}{2}\rm{ln}(2\pi)-\it{n}\rm{ln}(\sigma)-\frac{(y_{1}-\mu)^2}{2\sigma^2}-\ldots-\frac{(y_{n}-\mu)^2}{2\sigma^2}\]

因為\(\rm{log}(a^k)=k*\rm{log}(a)\)，所以式 (4.3)的前半部分可以改寫為：

\[\begin{align*} -\frac{n}{2}\rm{ln}(2\pi)-\it{n}\rm{ln}(\sigma)&= -\frac{n}{2}\rm{ln}(2\pi)-\frac{2n}{2}\rm{ln}(\sigma)\\ &=-\frac{n}{2}\rm{ln}(2\pi)-\frac{n}{2}\rm{ln}(\sigma^2)\\ &=-\frac{n}{2}\rm{ln}(2\pi*\sigma^2) \end{align*}\]

所以常態分佈的log-likelihood函數可寫成：

\[\mathcal{L}(X|\mu,\sigma^2)=\frac{-n}{2}\rm{ln}(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum (x_{i}-\mu)^2\]

參考Joel S Steele畫圖 4.9顯示常態分佈，確認機率密度函數無誤。

# specify the single value normal probability function
norm_lik = function(x, mu, s){
y = 1/sqrt(2*pi*s^2)*exp((-1/(2*s^2))*(x-mu)^2)
}

# and plot it just to make sure
xval=seq(-3,3, length.out = 600)
yval=sapply(seq(-3,3, length.out=600),FUN=norm_lik,mu=0,s=1)
plot(xval, yval, type='l', ylab='',xlab='', main='Normal Distribution',
     col='#FF22EE', lwd=1.5)

Figure 4.9: 常態分佈機率密度函數圖

另一種最大概似函數的寫法：

# Likelihood function
llik = function(x,par){
  m=par[1]
s=par[2]
n=length(x)
# log of the normal likelihood
# -n/2 * log(2*pi*s^2) + (-1/(2*s^2)) * sum((x-m)^2)
ll = -(n/2)*(log(2*pi*s^2)) + (-1/(2*s^2)) * sum((x-m)^2) 
# return the negative to maximize rather than minimize
return(-ll) }

用\(\tt{optim()}\) 極大化函數，結果如下：

set.seed(270511)
xvec<-rnorm(600, 1.5, sqrt(pi))
# call optim with the starting values 'par', # the function (here 'llik'),
# and the observations 'x'
res0 = optim(par=c(.5,.5), llik, x=xvec)
#Result
res=rbind('observation'=c('mean'=mean(xvec),
                              'sd'=sd(xvec)),
        'optim'=res0$par)

knitr::kable(res)%>%
  kable_styling(bootstrap_options = "striped", full_width = F)

	mean	sd
observation	1.54	1.768
optim	1.54	1.767

結果估計是平均值為1.53時、標準差=1.77時最大化，剛好就是我們的觀察值資料。
極大化的過程中，最大概似值越大表示資料越符合特定的分佈，或者是- log likelihood越小越好。

5 估計二元勝算對數迴歸模型係數

就像圖4.8顯示，要估計二元變數迴歸模型的係數需要找出函數的梯度(gradient)等於0，這要經過許多次的疊代(iteration)才能找到等於0的值。代數很難估計出非線性模型的參數，而是要用數值方法(numerical methods)，方法是從一個值\(\theta_{0}\)開始，加上\(\xi_{0}\)改善結果，變成\(\theta_{1}=\theta_{0}+\xi_{0}\)，依此類推，直到\(\theta_{n+1}=\theta_{n}+\xi_{n}\)，讓最大概似函數的梯度等於0（水平線），而\(\theta_{n+1}\)與\(\theta_{n+2}\)並沒有太大差異，也就是達到收斂(convergence)，此時\(\hat{\theta}\)就是我們的最大概似估計值。

\(\xi_{n}\)包括兩個部分，一個是最大概似函數的曲度，也就是變動X一單位，Y變動的程度。另一個是方向。有兩個比較常用的演算法：Newton-Raphson以及BHHH。

那麼我們如何根據(3.3)以及方程式(4.8)，應用R進行最大概似法求出二元變數迴歸模型的係數？換句話說，我們的最大概似函數是什麼？

首先複習\(a\cdot\rm{ln}(x)=\rm{ln}(x^a)\)，如果\(a=-1\)，那麼\(\rm{ln}(\frac{1}{x})=\rm{ln}(x^{-1})=-\rm{ln}(x)\)。
其次，根據(3.3)，以下的關係成立：

\[p(y;X,\beta)=(1+\rm{exp}(-X\beta))^{-1}\]

根據二項分配的機率密度函數，\(\beta\)的最大概似函數可寫成：（參考Steven Wu）

\[\begin{align*} L(\beta)&=\rm{ln}\prod\limits_{i=1}^{n}p^{y_{i}}(1-p)^{1-y_{i}}\\ &=\sum\limits_{i=1}^{n}(y_{i}\rm{ln}(p)+(1-y_{i})\rm{ln}(1-p))\\ &=\sum\limits_{i=1}^{n}(y_{i}\rm{ln}(\frac{1}{1+\rm{exp}(-X\beta)})+(1-y_{i})\rm{ln}(1-(\frac{1}{1+\rm{exp}(-X\beta)})))\\ &=-\sum\limits_{i=1}^{n}(y_{i}\rm{ln}(1+\rm{exp}(-X\beta))+(1-y_{i})\rm{ln}(1+\rm{exp}(X\beta)) \tag{5.1} \end{align*}\]

根據(5.1)，參考Christopher Adolph，估計二元勝算對數迴歸模型參數的語法如下：

set.seed(02139)
# Dep. V.
y<- rbinom(1000, 1, prob = 0.67)
# Ind. V.
x <- cbind(X1=rnorm(1000, 1, 0.2))
# Likelihood function for logit
llk.logit <- function(param,y,x) {
  os <- rep(1,length(x[,1]))
  x <- cbind(os,x)
  b <- param[ 1 : ncol(x) ]
  xb <- x%*%b
  sum( y*log(1+exp(-xb)) + (1-y)*log(1+exp(xb)));
               # optim is a minimizer, so min -ln L(param|y)
}
# Fit logit model using optim

ls.result <- lm(y~x)  # use ls estimates as starting values
stval <- ls.result$coefficients  # initial guesses
logit.result <- optim(stval,llk.logit,method="BFGS",hessian=T,y=y,x=x)
                   # call minimizer procedure
pe.1 <- logit.result$par   # point estimates
vc.1 <- solve(logit.result$hessian)  # var-cov matrix
se.1 <- sqrt(diag(vc.1))    # standard errors
ll.1 <- -logit.result$value  # likelihood at maximum
#Result
res=rbind('Estimates'=rev(pe.1),
        'Std. Error'=rev(se.1),
        'Log likelihood'=ll.1)

knitr::kable(res)%>%
  kable_styling(bootstrap_options = "striped", full_width = F)

	x	(Intercept)
Estimates	-0.1773	0.9556
Std. Error	0.3310	0.3410
Log likelihood	-622.8996	-622.8996

我們用\(\tt{glm()}\)函數估計同樣的模型，比較表5.1，結果相同。

m1<-glm(y~x, family = binomial('logit'))
stargazer::stargazer(m1, title='二元勝算對數迴歸模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 5.1: **二元勝算對數迴歸模型**

	Dependent variable:

	y

x	-0.177
	(0.331)

Constant	0.956^***
	(0.341)


Observations	1,000
Log Likelihood	-622.900
Akaike Inf. Crit.	1,250.000

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

#logLik(m1)

5.1 二元勝算對數迴歸模型實例

\(\blacksquare\) 個人財務會不會跟無法償還信用卡的債務有關？我們用\(\tt{glm()}\)函數估計二元勝算迴歸模型，注意要設定\(\tt{family=binomial(logit)}\)。

library(ISLR)
M2 <- glm(default ~ balance, data=Default, family=binomial(logit))
stargazer::stargazer(M2,  title='信用卡債務與財務',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
      label=knitr::opts_current$get("label"))

Table 5.2: **信用卡債務與財務**

	Dependent variable:

	default

balance	0.005^***
	(0.0002)

Constant	-10.650^***
	(0.361)


Observations	10,000
Log Likelihood	-798.200
Akaike Inf. Crit.	1,600.000

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

根據式(3.2)，用勝算比的概念，表 5.2 估計結果可寫成方程式為：

\[\begin{align*} \frac{p(X)}{1-p(X)} & = exp^{-10.651+0.0055\cdot \mathtt{balance}}\\ & = exp^{-10.651}\cdot exp^{0.0055\cdot \mathtt{balance}} \end{align*}\]

由5.2可知，\(\hat{\beta}_{1}=0.0055\)，代表增加一個單位的信用卡扣掉貸款的消費金額(balance)，會增加違約(default)的指數勝算(log odds)0.0055，或者是違約對不違約的勝算增加\(\text{exp}^{0.0055}=1.005515\)。
我們畫圖 5.1 表示\(\frac{p(X)}{1-p(X)}\)與自變數的關係：

a=-10.651; b=0.0055
y=exp(a)*exp(b*ISLR::Default$balance)
yi=sort(y)
xi=sort(ISLR::Default$balance)
dt <-data.frame(yi, xi)

ggplot(dt, aes(x=xi, y=yi)) +
  geom_point(col='#CC9933')

Figure 5.1: 自變數與勝算的關係

另外畫圖 5.2 表示\(p(X)\)與\(\frac{e^\eta}{1+e^\eta}\)的關係：

a=-10.651; b=0.0055
xi=sort(ISLR::Default$balance)
ETA=a + b*xi
pi=exp(ETA)/(1+exp(ETA))

dt <-data.frame(x=xi, y=pi)

ggplot(dt, aes(x=xi, y=pi)) +
  geom_point(col='#FF9933')

Figure 5.2: 機率與指數的關係

圖 5.2 顯示，\(p\)值隨著自變數的增加而上升，但是不會小於0或是大於1。中間有一段接近線性函數。當\(p\)接近0，\(\rm{log}\frac{p(X)}{1-p(X)}\approx -\infty\)，當\(p\)接近1，\(\rm{log}\frac{p(X)}{1-p(X)}\approx \infty\)。

5.2 勝算比(odds ratio)

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

\[ \text{odds(male)}=\frac{0.8}{0.2}=4 \] \[ \text{odds(female)}=\frac{0.2}{0.8}=0.25 \]

根據上述，男性參加抽獎是女性的16倍，因為\(\frac{4}{0.25}=16\)
那麼logit模型跟勝算比有什麼關係？因為經過轉換之後，\(\beta\)是勝算比的對數值。

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

所以\(\text{exp}^{\beta_{1}}\)等於當=1是x=0的y=1的倍數。這個解釋特別適用於自變數是二元類別的模型。
例如表 5.3 納入自變數是有無學生身份，預測是否發生信用卡違約的模型為：

library(ISLR)
#Default$default <- as.numeric(ISLR::Default$default)-1
M3 <- glm(default ~ student, data=Default, family=binomial(logit))
stargazer::stargazer(M3,  title='學生身份與信用卡違約',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 5.3: **學生身份與信用卡違約**

	Dependent variable:

	default

studentYes	0.405^***
	(0.115)

Constant	-3.504^***
	(0.071)


Observations	10,000
Log Likelihood	-1,454.000
Akaike Inf. Crit.	2,913.000

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

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

5.2.1 實例

假設有一個自變數的值為1到5，估計logit模型得到\(\rm{log}\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。

5.2.2 雙變數二元勝算對數迴歸模型

我們也可以同時估計兩個自變數與依變數的關係如表 5.4：

M4 <- glm(default ~ balance + student, data=Default, family=binomial(logit))

stargazer::stargazer(M4,  title='雙自變數的Logit模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 5.4: **雙自變數的Logit模型**

	Dependent variable:

	default

balance	0.006^***
	(0.0002)

studentYes	-0.715^***
	(0.148)

Constant	-10.750^***
	(0.369)


Observations	10,000
Log Likelihood	-785.800
Akaike Inf. Crit.	1,578.000

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

比較前後三個模型，balance的係數幾乎一樣，但是studentYes從正變成負，也就是對於同樣的信用卡消費金額扣掉貸款之後，學生反而不會違約。

6 Probit 模型的原理

常態分佈的隨機變數\(X\)有對應常態分佈的機率(PDF)，以及累積機率(CDF)，可寫成\(P(X < x)\)，也就是累積到\(x\)之前的機率。也可表示為\(F(X)\)或者是\(\Phi(\cdot)\)。

對於任何的常態分佈的隨機變數，\(F\)函數可以轉換為標準化的累積機率分佈。

標準化常態分佈的隨機變數\(X\)經常表示為\(Z\)，所以\(Z\)服從標準化常態分佈： \(Z \sim \mathcal{N}(\mu = 0, \hspace{0.1cm}\sigma = 1)\)

例如我們想計算X<0.25的標準化CDF，而且\(\mu=1\)、\(\sigma^2=2\)。計算方式為： \(P(X < 0.25) = P(\frac{X - \mu}{\sigma} < \frac{0.25-\mu}{\sigma}) = P(Z < \frac{0.25 - 1}{2}) = \Phi\left(\frac{0.25 - 1}{2}\right) = 0.3538\)

x=0.25; mu=1; sd=2
pnorm(x, mu, sd)
## [1] 0.3538

我們可以畫圖6.1 如下：

x<-seq(-5, 5, length.out=1000)
CDF <- c()
for (i in 1: 1000){
   CDF[i] = pnorm(x[i], 1, 2)
}
dt <- data.frame(x, CDF)
library(ggplot2)
ggplot(data=dt, aes(x=x, y=CDF)) +
      geom_line(col='#CC33FF', size=1.5)

Figure 6.1: 累積機率密度圖

圖 6.1 顯示，標準化常態分佈的累積機率密度介於0與1之間。可以寫成\(\Phi(Z)\in [0,1]\)。當\(x=0\)，剛好累積\(50\%\)。
由於依變數為二元變數\(Y_{i}=0\)或1，服從伯努利分佈，參數為次數\(N\)與機率\(P\)，表示成\(P_{i}(Y_{i}=1)\)

logit轉換依變數\(Y_{i}\)變成\(P_{i}\)，介於\(\{0, 1\}\)。而log-odds再轉換成\(\{-\infty, \infty \}\)，標準常態分佈的累積機率，又稱probit或者成長曲線迴歸（黃紀、王德育，2012:87），依變數仍然是\(\{0, 1\}\)。稱為二元機率單元模型。

標準常態分佈可寫成:

\[ \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 轉換\(X\beta\)成為常態分佈的\(Z\)分數，所以\(X\beta\)越大，Pr(Y=1)越大。-
Probit模型可寫成：

\[\rm{Pr}(y=1|X)=\Phi(X\beta)\]

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

我們模擬重複3000次的實驗，每次有1000個樣本，而成功的機率為0.1到0.99的一致分佈中隨機抽取一個數。如果1000個樣本中，至少有\(k*p*d\)的成功次數，其中d是常態分佈隨機抽出的亂數，便把這個實驗結果稱為1，也就是有3000個實驗結果。再以常態分佈產生自變數\(X\)，並且以probit函數估計自變數的作用。

set.seed(100)
n=3000; k=1000; p=runif(1,0.1,0.99); d=rnorm(1, 1, 0.01)
V=rbinom(n, k, p)
Vy<-c()
for (i in 1: n){
         if(V[i]>=k*p*d)
           Vy[i]=1
          else(Vy[i]=0)
}
table(Vy)
## Vy
##    0    1 
## 1289 1711

表 6.1呈現估計的結果。

X=rnorm(n, p*d, sqrt(p*(1-p)))
fit1 <- glm(Vy ~ X, family=binomial(probit) )

stargazer::stargazer(fit1,  title='模擬資料中的自變數作用',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 6.1: **模擬資料中的自變數作用**

	Dependent variable:

	Vy

X	-0.094^**
	(0.048)

Constant	0.211^***
	(0.029)


Observations	3,000
Log Likelihood	-2,048.000
Akaike Inf. Crit.	4,099.000

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

接著代入係數，求出\(Z(X\beta)\)，也就是\(Z\)分數，見圖 6.2。

set.seed(100)
n=3000; k=1000; p=runif(1,0.1,0.99); d=rnorm(1, 1, 0.01)
newCDF <- c()
xval=sort(X)
for (i in 1: n){
   newCDF[i] = pnorm(0.210-0.094*X[i],  p*d, sqrt(p*(1-p)))
}
Z <-sort(newCDF)
dt<-data.frame(Vy, X, xval, Z)

library(ggplot2)
G=ggplot(data=dt, aes(x=xval, y=Vy)) +
   geom_point(size=1.5) 
G+geom_line(data=dt, aes(x=xval, y=Z), 
            col='#33EEFF', size=1.2)

Figure 6.2: 模擬資料之Z值圖

6.1 Probit 累積機率

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

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

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

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

可用R計算以上結果：

pnorm(-2, 0, 1)
## [1] 0.02275
pnorm(-1, 0, 1)
## [1] 0.1587
pnorm(2, 0, 1)
## [1] 0.9772

計算全部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)

Figure 6.3: 標準化常態分佈累積機率密度

由以上可知，\(0<\Phi(Z)<1\)，而且累積的機率形成S形的曲線。
當x增加一個單位，Pr(y=1) 增加\(\beta\)的Z值。
Probit的累積機率(cdf)可畫成如下：

v1<-here::here('Fig','probitcdf.png')
knitr::include_graphics(v1)

Figure 6.4: Probit累積機率函數

6.2 估計二元單元機率模型

根據式(3.3)，我們嘗試估計二元單元機率模型參數如下：（參考Sol\(\grave{e}\) Prillaman）

# Declare the number of observations to create.
nobs=1000
x1 = rnorm(nobs)*.15^.5 
x2 = rnorm(nobs)*.35^.5 
u = rnorm(nobs)*(.5)^.5
inside = 2*x1+2*x2+u

p = pnorm(inside)
y = rbinom(nobs,1,p)
mydata = data.frame(y=y,x1=x1,x2=x2)
X <- as.matrix(mydata[,2:3])
ll.probit <- function(beta, y=y, X=X){
  if(sum(X[,1]) != nrow(X)) X <- cbind(1,X)
  phi <- pnorm(X%*%beta, log = TRUE)
  opp.phi <- pnorm(X%*%beta, log = TRUE, lower.tail = FALSE)
  logl <- sum(y*phi + (1-y)*opp.phi)
  return(logl)
}
## Estimate the MLE:
p.res <- optim(par = rep(0, ncol(X) + 1),
             fn = ll.probit,
             y = y,
             X = X,
             control = list(fnscale = -1),
             hessian = T,
             method = "BFGS")

names(p.res$par) <- c('Intercept', colnames(X))

# Extract the coefficients from optim
coefs <- p.res$par
# Calculate the variance
varcov <- -solve(p.res$hessian)
# Calculate the standard errors
ses <- sqrt(diag(varcov))
# Log lik
ll.1 <- p.res$value

#Result

res=rbind('Estimates' = coefs,
        'Std. Error'= ses,
        'Log likelihood'=ll.1)

knitr::kable(res)%>%
  kable_styling(bootstrap_options = "striped", full_width = F)

	Intercept	x1	x2
Estimates	0.0001	1.4449	1.7239
Std. Error	0.0475	0.1319	0.1116
Log likelihood	-462.3379	-462.3379	-462.3379

用\(\tt{glm()}\)函數估計同樣的模型，比較表6.2，結果相同。

pm1<-glm(y~x1+x2, family = binomial('probit'))
stargazer::stargazer(pm1, title='二元單元機率迴歸模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 6.2: **二元單元機率迴歸模型**

	Dependent variable:

	y

x1	1.445^***
	(0.136)

x2	1.724^***
	(0.111)

Constant	0.0002
	(0.047)


Observations	1,000
Log Likelihood	-462.300
Akaike Inf. Crit.	930.700

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

#logLik(m1)

6.3 Probit模型係數的詮釋

參考UCLA的申請入學例子，依變數與自變數如表 6.3：

#ucla.url <- "https://stats.idre.ucla.edu/stat/data"
#admitlink<- paste(ucla.url, "binary.csv", sep = "/")
#download.file(url=admitlink, destfile = "./data/uclaadmit.csv" )
addt <- here::here('data', 'binary.txt')
Admit <-read.table(addt, header = T, sep = ',')
Admit$rank.f<- factor(Admit$rank)
knitr::kable(head(Admit),format = 'simple',
      booktabs=T,caption='UCLA申請入學')

Table 6.3: UCLA申請入學
admit	gre	gpa	rank	rank.f
0	380	3.61	3	3
1	660	3.67	3	3
1	800	4.00	1	1
1	640	3.19	4	4
0	520	2.93	4	4
1	760	3.00	2	2

接下來用glm函數估計自變數的係數，估計是否被錄取的模型如表 6.4：

fitadmit <- glm(admit ~ gre + gpa + rank.f, data = Admit, family = binomial("probit"))
stargazer::stargazer(fitadmit,  title='被錄取與否的Probit模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 6.4: **被錄取與否的Probit模型**

	Dependent variable:

	admit

gre	0.001^**
	(0.001)

gpa	0.478^**
	(0.197)

rank.f2	-0.415^**
	(0.195)

rank.f3	-0.812^***
	(0.208)

rank.f4	-0.936^***
	(0.245)

Constant	-2.387^***
	(0.674)


Observations	400
Log Likelihood	-229.200
Akaike Inf. Crit.	470.400

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

雖然我們可以說在其他變數相同的情況下，GRE分數增加一個單位，\(Z\)分數增加0.001，但是計算累積機率密度更能看出自變數變化時，依變數變化的程度。
模型估計結果為：

\[F(-2.387+0.001*\rm{GRE}+0.478*\rm{GPA}-0.415*\rm{rank2}-0.812*\rm{rank3}-0.936*\rm{rank4})\]

假設GPA成績為0、學生屬於第一等級，當GRE分數為500，累積機率密度為：

\[F(-2.387+0.001*500)=0.0295\]

用R計算：

pnorm(-2.387+0.001*500)
## [1] 0.02958

假設GPA成績為0、學生屬於第一等級，當GRE分數為600，累積機率密度為：

\[F(-2.387+0.001*600)=0.0369\]

pnorm(-2.387+0.001*600)
## [1] 0.03697

假設GPA成績為0、學生屬於第一等級，當GRE分數為700，累積機率密度為：

\[F(-2.387+0.001*600)=0.0458\]

pnorm(-2.387+0.001*700)
## [1] 0.0458

以上計算方式可以代入GPA的平均值，或者其他的GRE分數，得到不同的累積機率密度。
畫圖 6.5 顯示不同的GPA情況下，當GRE分數增加時，不同等級的學生的錄取機率增加的趨勢：

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.f= factor(rep(rep(1:4, 
    each = 100), 4)))
newdata[, c("p", "se")] <- predict(fitadmit, newdata, type = "response", se.fit = TRUE)[-3]

ggplot(newdata, aes(x = gre, y = p, colour = rank.f)) +   geom_line() + 
  facet_wrap(~gpa)

Figure 6.5: 綠取機率與分數的關係

7 信賴區間

當依變數\(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)}\]

將\(\hat{\pi}\)代入\(\pi=\frac{1}{1+exp{-X\beta}}\)，應用疊代再加權最小平方(Iteratively Reweighted Least Squared, IRLS)的演算式，聚合以下算式：

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

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

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

再應用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一樣的估計。特別注意要去掉遺漏值。

7.1 MLE的應用

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

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

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

7.2 假設檢定

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

檢定虛無假設\(\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兩種。

7.2.1 Wald test

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

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

7.2.2 Liklihood ratio (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\))百分位。

8 多變數二元勝算對數迴歸模型

我們可以延伸上述的logit迴歸模型到多變數： \[ \text{log}(\frac{p(X)}{1-p(X)})=\beta_{0}+\beta_{1}X_{1}+\ldots+\beta_{p}X_{p} \]
因為多於一個變數，所以在詮釋其中一個變數的影響時，需要強調其他的變數在同樣水準。
在計算預測值時，必須設定變數特定的值，除了我們所關心的自變數。例如平均值、眾數等等。

8.1 實例

\(\blacksquare\)哪些因素影響信用卡違約？用二元勝算對數迴歸模型估計如表 8.1：


M5 <-glm(default ~ balance+income+student, data=Default, family=binomial(logit))

stargazer::stargazer(M5,  title='影響信用卡違約的模型',
  type=ifelse(knitr::is_latex_output(),"latex","html"),
label=knitr::opts_current$get("label"))

Table 8.1: **影響信用卡違約的模型**

	Dependent variable:

	default

balance	0.006^***
	(0.0002)

income	0.00000
	(0.00001)

studentYes	-0.647^***
	(0.236)

Constant	-10.870^***
	(0.492)


Observations	10,000
Log Likelihood	-785.800
Akaike Inf. Crit.	1,580.000

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

同樣收入、信用卡積欠金額的學生比非學生來得不會違約，統計上達到顯著水準。信用卡積欠金額同樣具有統計上顯著的影響。

8.1.1 指數模型

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

M5.or <- exp(coef(M5))
knitr::kable(M5.or, format = 'simple',
             caption='被錄取與否的模型')

Table 8.2: 被錄取與否的模型
	x
(Intercept)	0.0000
balance	1.0058
income	1.0000
studentYes	0.5237

#stargazer::stargazer(M5.or, coef=list(M5.or),  title='被錄取與否的模型',
#  type=ifelse(knitr::is_latex_output(),"latex","html"),
#label=knitr::opts_current$get("label"))

8.1.2 預測值

\(\tt{predict()}\)函數可計算模型的預測值。
先建立一個資料框，包含三個自變數的平均值、最大值或者眾數。

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

balance	income	student
835.4	33517	Yes
835.4	33517	No

代入模型M5： \[\rm{Pr}(Y=1)=(e^{\eta})/(1+e^\eta)\] 其中 \[\eta=-10.869+0.006*\rm{balance}-0.647*\rm{student}\]

M5 <-glm(default ~ balance+income+student, data=Default, family=binomial(logit))
allcoef$pred <-predict(M5, newdata=allcoef, type="response")
allcoef<-cbind(allcoef, allcoef$pred)
allcoef
##   balance income student     pred allcoef$pred
## 1   835.4  33517     Yes 0.001329     0.001329
## 2   835.4  33517      No 0.002534     0.002534

8.1.3 比較Logit, Probit估計結果

表 8.3 顯示Logit, Probit模型的估計結果：

M5 <-glm(default ~ balance+income+student, data=Default, family=binomial(logit))
M6 <-glm(default ~ balance+income+student, data=Default, family=binomial(probit))
stargazer::stargazer(M5, M6, align=TRUE, title='Logit, Probit模型的估計比較',
  type=ifelse(knitr::is_latex_output(),"latex","html"), label=knitr::opts_current$get("label"))

Table 8.3: **Logit, Probit模型的估計比較**

	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.870^***	-5.475^***
	(0.492)	(0.238)


Observations	10,000	10,000
Log Likelihood	-785.800	-791.600
Akaike Inf. Crit.	1,580.000	1,591.000

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

Logistic迴歸模型的係數以及標準誤大約是Probit模型的兩倍。

9 比較不同logit或者Probit模型的係數

Allison (1999) 提出如何比較不同logit或者Probit模型的係數的方法，見Allison, Paul D. (1999). Comparing Logit and Probit Coefficients Across Groups. Sociological Methods & Research, 28(2), 186–208. https://doi.org/10.1177/0049124199028002003

10 小結

轉換二元變數後，極大化最大概似法函數，估計出Logit, Probit模型的參數。

以LR test 檢定迴歸模型。

以指數詮釋二元勝算對數迴歸模型的係數。

11 作業

1. 請隨機產生200個伯努利分佈的亂數，而且母體的機率是0.45。請問200個樣本的平均數是多少？請用最大概似法求出平均值。

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

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)))
knitr::kable(DT, "latex", caption = "Admission and Gender")

3. 請估計以下資料中，暴露於某種物質(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<-DT %>% 
  kable("latex", caption = "Survive and Age")

4. 接續上題，請比較probit與logistic的模型估計結果。

12 更新日期

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

社會科學統計方法

二元變數迴歸模型

蔡佳泓

5/27/2022