邏輯迴歸研究因變數Y為分類變數與多個自變數X之間的迴歸問題。隨機變數X的取值為實數，隨機變數Y的取值為1或0。常用於預測某隨機事件發生概率的大小。

Logistic迴歸問題的最優化問題可以表述為：

尋找一個非線性函式Sigmoid的最佳擬合引數，求解過程可使用最優化演算法完成。它可以看做是用Sigmoid函式作為二閾值分類器的感知器問題。

一、logistic迴歸及其MLE

當我們考慮解釋變數為分類變數，如考慮一個企業是否會被併購，一個企業是否會上市這些問題時，考慮線性概率模型 P(yi=1)=β0+β1xi 顯然是不合適的，它至少有兩個致命的缺陷：

1、概率估計值可能超過1，使得模型失去了意義;(要解決這個問題並不麻煩，我們將預測超過1的部分記為1，低於0的

2、邊際效應假定為不變，通常來說不合經濟學常識。考慮一個邊際效應遞減的模型(假定真實值為藍線)，可以看到線性模型表現很差。

但是sigmoid函式去擬合藍線確實十分合適的。Sigmoid函式形式為

圖形如下所示

於是我們可以考慮logistic迴歸模型：

對輸入X進行分類的線性函式，其值域為實數域，通過Logistic迴歸模型定義式將線性函式轉換為概率。線性函式的值越接近於正無窮，概率值越接近於1；線性函式越接近於負無窮，概率值越接近於0。

假定有N個觀測樣本Y1,Y2,…,YN，設P(Yi=1|Xi)=π(Xi)為給定條件Xi下得到結果Yi=1的條件概率;而在同

1-Yi假設各觀測獨立，對logistic迴歸模型來說，其對數似然函式為：

因為

所以

將問題變成了以對數似然函式為目標函式的最優化問題，可以採用梯度下降法或擬牛頓法求解出logistic模型的MLE。

二、logit還是probit?

雖說Sigmoid函式對邊際遞減的模型擬合良好，但是我們也要知道S型函式並非僅sigmoid函式一個，絕大多數的累積分佈函式都是S型的。

於是考慮F-1(P)(F為標準正態分佈的累積分佈函式)也不失為一個很好的選擇。像這樣的，對概率P做一點變換，讓變換後的取值範圍變得合理，且變換後我們能夠有辦法進

logistic迴歸認為二分類變數服從伯努利分佈，應當選擇logit，而且從解釋的角度說，p/(1-p)就是我們常說的機率（odds ratio），也就是軟體報告中出現的OR值。但是probit也有它合理的一面，首先，中心極限定理告訴我們，伯努利分佈在樣本夠多的時候就是近似正態分佈的;其次，從不確定性的角度考慮，probit認為我們的線性概率模型服從正態分佈，這也是更為合理的。

我們來看一下經過變換後，自變數和P的關係是什麼樣子的：

如果你確實想知道到底你的資料用哪一個方法好，也不是沒有辦法，你可以看一下你的殘差到底是符合logit函式呢還是符合probit函式，當然，憑肉眼肯定是看不出來的，因為這兩個函式本來就很接近，你可以通過函式的假定，用擬合優度檢驗一下。但通常，估計不會有人非要這麼較真，因為沒有必要。但是有一點是要注意的，logit模型較probit模型而言具有厚尾的特徵，這也是為什麼經濟學論文愛用logit的原因。

我們以鳶尾花資料中的virginica,versicolor兩類資料分類為例，看看兩種辦法分類有無差別。

probit.predictions

versicolor virginica

versicolor 47 3

virginica 3 47

logit.predictions

versicolor virginica

versicolor 47 3

virginica 3 47

從上圖與比較表格均可以看出，兩者差別不大。

三、多項式logit與order logit

對於二項分類模型的一個自然而然的推廣便是多項分類模型。

我們借鑑神經網路裡提到的異或模型，有：

按照上面這種方法，給定一個輸入向量x，獲得最大hθ(x)的類就是x所分到的類。

選擇最大的 hθ(x)十分好理解：在類別選擇問題中，不論要選的類別是什麼，每一個類別對做選擇的經濟個體來說都有或多或少的效用(沒有效用的類別當然不會被考慮) ，一個類別的脫穎而出必然是因為該類別能產生出最高的效用。

我們將多項logit模型的數學表述敘述如下：

假定對於第i個觀測，因變數Yi有M個取值：1,2，…，M，自變數為Xi，則多項logit模型為：

與logistic迴歸的似然估計類似，我們可以很容易寫出多項logit的對數似然函式：

多項 Logit模型雖然好用，但從上面的敘述可以看出，多項 Logit 模型最大的限制在於各個類別必須是對等的，因此在可供選擇的類別中，不可有主要類別和次要類別混雜在一起的情形。例如在研究旅遊交通工具的選擇時，可將交通工具的類別粗分為航空、火車、公用汽車、自用汽車四大類，但若將航空類別再依三家航空公司細分出三類而得到總共六個類別，則多項 Logit 模型就不適用，因為航空、火車、公用汽車、自用汽車均屬同一等級的主要類別，而航空公司的區別則很明顯的是較次要的類別，不應該混雜在一起。在這個例子中，主要類別和次要類別很容易分辨，但在其他的研究中可能就不是那麼容易，若不慎將不同層級的類別混在一起，則由多項 Logit 模型所得到的實證結果就會有誤差。

對於分類模型，我們還會遇到被解釋變數中有分類變數的情形。對於連續變數解釋離散變數，且被解釋的離散變數是有順序的(這個是和多項logit最大的區別)的情形，我們就需要考慮到order logit模型。

其數學模型敘述如下：

其中，F(.)表示累積分佈函式，當F表示正態分佈的分佈函式時，對應的是order probit;F表示logistic分佈時，變對應order logit。

與logistic分佈類似，我們可以很容易寫出其對數似然函式：

四、啞變數(dummy variable)

在logistic迴歸中，經常會遇到解釋變數為分類變數的情形，比如收入：高、中、低;地域：北京、上海、廣州等。這裡對分類變數而言就涉及一個問題：要不要將分類變數設定dummy variable?

這個問題的答案線上性模型中很顯然，必須要這麼做!!!如果我們不設定啞變數，而是單純地賦值：北京=1，上海=2，廣州=3，即我們將自變數視作連續性的數值變數，但這僅僅是一個程式碼而己，並不意味著地域間存在大小次序的關係，即並非代表被解釋變數(響應變數)會按此順序線性增加或減少。即使是有序多分類變數，如家庭收入分為高、中、低三檔，各類別間的差距也是無法準確衡量的，按編碼數值來分析實際上就是強行規定為等距，這顯然可能引起更大的誤差。

但是在logistic迴歸中，由於logit(p)變化的特殊性，在解釋定序變數時，為了減少自由度(即解釋變數個數)，我們常常將定序變數(如家庭收入分為高、中、低)視為連續的數值變數，而且經濟解釋可以是XX每提高一個檔次，相應的概率會提高expression(delta(XX))(expression的表示式還是很複雜的，不打了)。當然減少變數個數是以犧牲預測精度為代價的。畢竟資料處理是一門藝術而非一門技術，如何取捨還得具體問題具體分析。當然，非定序的分類變數是萬萬不可將其視為數值變數的。

五、廣義線性模型的R實現

R語言提供了廣義線性模型的擬合函式glm()，其呼叫格式如下：

glm(formula, family = gaussian, data,weights, subset,na.action, start = NULL, etastart, mustart, offset,control= list(...), model = TRUE, method = "glm.fit",x =FALSE, y = TRUE, contrasts = NULL, ...)

引數說明：

Formula：迴歸形式，與lm()函式的formula引數用法一致

Family：設定廣義線性模型連線函式的典則分佈族，glm()提供正態、指數、gamma、逆高斯、Poisson、二項分

布。我們的logistic迴歸使用的是二項分佈族binomial。Binomial族預設連線函式為logit，可設定為probit。

Data：資料集

鳶尾花例子使用的R程式碼：

logit.fit <- glm(Species~Petal.Width+Petal.Length,

family = binomial(link = 'logit'),

data = iris[51:150,])

logit.predictions <- ifelse(predict(logit.fit) > 0,'virginica', 'versicolor')

table(iris[51:150,5],logit.predictions)

probit.fit <- glm(Species~Petal.Width+Petal.Length,

family = quasibinomial(link = 'probit'),

data = iris[51:150,])

probit.predictions <- ifelse(predict(probit.fit) >0,'virginica', 'versicolor')

table(iris[51:150,5],probit.predictions)

程式包mlogit提供了多項logit的模型擬合函式：

mlogit(formula, data, subset, weights,na.action, start = NULL,

alt.subset = NULL, reflevel = NULL,

nests = NULL, un.nest.el = FALSE, unscaled = FALSE,

heterosc = FALSE, rpar = NULL, probit = FALSE,

R = 40, correlation = FALSE, halton = NULL,

random.nb = NULL, panel = FALSE, estimate = TRUE,

seed = 10, ...)

mlogit.data(data, choice, shape = c("wide","long"), varying = NULL,

sep=".",alt.var = NULL, chid.var = NULL, alt.levels = NULL,id.var = NULL, opposite = NULL, drop.index = FALSE,

ranked = FALSE, ...)

引數說明：

formula：mlogit提供了條件logit，多項logit，混合logit多種模型，對於多項logit的估計模型應寫為：因變數~0|自變數,如：mode ~ 0 | income

data：使用mlogit.data函式使得資料結構符合mlogit函式要求。

Choice：確定分類變數是什麼

Shape：如果每一行是一個觀測，我們選擇wide，如果每一行是表示一個選擇，那麼就應該選擇long。

alt.var：對於shape為long的資料，需要標明所有的選擇名稱

選擇wide的資料示例：

選擇long的資料示例：

以fishing資料為例，來說明如何使用mlogit。

library(mlogit)

data("Fishing", package = "mlogit")

Fish <- mlogit.data(Fishing,shape = "wide",choice = "mode")

summary(mlogit(mode ~ 0 | income, data = Fish))

這個輸出的結果與nnet包中的multinom()函式一致。由於mlogit包可以做的logit模型更多，所以這裡就不在對nnet

包的multinom作介紹了，可以參見《根據Econometrics in R一書，將回歸方法總結一下》一文。

程式包MASS提供polr()函式可以進行ordered logit或probit迴歸。用法如下：

polr(formula, data, weights, start, ..., subset, na.action,

contrasts = NULL, Hess = FALSE, model = TRUE,

method = c("logistic", "probit", "cloglog", "cauchit"))

引數說明：

Formula：迴歸形式，與lm()函式的formula引數用法一致

Data：資料集

Method：預設為order logit，選擇probit時變為order probit模型。

以housing資料為例說明函式用法：

house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)

house.plr

summary(house.plr, digits = 3)

這些結果十分直觀，易於解讀，所以我們在這裡省略所有的輸出結果。

再看手寫數字案例：

最後，我們回到最開始的那個手寫數字的案例，我們試著利用多項logit重做這個案例。(這個案例的描述與資料參見《kNN演算法》一章)

特徵的選擇可參見《神經網路》一章。

由於手寫數字的特徵選取很容易導致迴歸係數矩陣是降秩的，所以我們使用nnet包的multinom()函式代替mlogit()。

執行下列程式碼：

setwd("D:/R/data/digits/trainingDigits")

names<-list.files("D:/R/data/digits/trainingDigits")

data<-paste("train",1:1934,sep="")

for(i in 1:length(names))

assign(data[i],as.matrix(read.fwf(names[i],widths=rep(1,32))))

label<-factor(rep(0:9,c(189,198,195,199,186,187,195,201,180,204)))

feature<-matrix(rep(0,length(names)*25),length(names),25)

for(i in 1:length(names)){

feature[i,1]<-sum(get(data[i])[,16])

feature[i,2]<-sum(get(data[i])[,8])

feature[i,3]<-sum(get(data[i])[,24])

feature[i,4]<-sum(get(data[i])[16,])

feature[i,5]<-sum(get(data[i])[11,])

feature[i,6]<-sum(get(data[i])[21,])

feature[i,7]<-sum(diag(get(data[i])))

feature[i,8]<-sum(diag(get(data[i])[,32:1]))

feature[i,9]<-sum((get(data[i])[17:32,17:32]))

feature[i,10]<-sum((get(data[i])[1:8,1:8]))

feature[i,11]<-sum((get(data[i])[9:16,1:8]))

feature[i,12]<-sum((get(data[i])[17:24,1:8]))

feature[i,13]<-sum((get(data[i])[25:32,1:8]))

feature[i,14]<-sum((get(data[i])[1:8,9:16]))

feature[i,15]<-sum((get(data[i])[9:16,9:16]))

feature[i,16]<-sum((get(data[i])[17:24,9:16]))

feature[i,17]<-sum((get(data[i])[25:32,9:16]))

feature[i,18]<-sum((get(data[i])[1:8,17:24]))

feature[i,19]<-sum((get(data[i])[9:16,17:24]))

feature[i,20]<-sum((get(data[i])[17:24,17:24]))

feature[i,21]<-sum((get(data[i])[25:32,17:24]))

feature[i,22]<-sum((get(data[i])[1:8,25:32]))

feature[i,23]<-sum((get(data[i])[9:16,25:32]))

feature[i,24]<-sum((get(data[i])[17:24,25:32]))

feature[i,25]<-sum((get(data[i])[25:32,25:32]))

}

data1 <- data.frame(feature,label)

#降秩時mlogit不可用

#data10<- mlogit.data(data1,shape = "wide",choice = "label")

#m1<-mlogit(label~0|X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13+X14+X15+X16+X17+X18+X19+X20+X21+X22+X23+X24+X25,data=data10)

library(nnet)

m1<-multinom(label ~ ., data = data1)

pred<-predict(m1,data1)

table(pred,label)

sum(diag(table(pred,label)))/length(names)

setwd("D:/R/data/digits/testDigits")

name<-list.files("D:/R/data/digits/testDigits")

data1<-paste("train",1:1934,sep="")

for(i in 1:length(name))

assign(data1[i],as.matrix(read.fwf(name[i],widths=rep(1,32))))

feature<-matrix(rep(0,length(name)*25),length(name),25)

for(i in 1:length(name)){

feature[i,1]<-sum(get(data1[i])[,16])

feature[i,2]<-sum(get(data1[i])[,8])

feature[i,3]<-sum(get(data1[i])[,24])

feature[i,4]<-sum(get(data1[i])[16,])

feature[i,5]<-sum(get(data1[i])[11,])

feature[i,6]<-sum(get(data1[i])[21,])

feature[i,7]<-sum(diag(get(data1[i])))

feature[i,8]<-sum(diag(get(data1[i])[,32:1]))

feature[i,9]<-sum((get(data1[i])[17:32,17:32]))

feature[i,10]<-sum((get(data1[i])[1:8,1:8]))

feature[i,11]<-sum((get(data1[i])[9:16,1:8]))

feature[i,12]<-sum((get(data1[i])[17:24,1:8]))

feature[i,13]<-sum((get(data1[i])[25:32,1:8]))

feature[i,14]<-sum((get(data1[i])[1:8,9:16]))

feature[i,15]<-sum((get(data1[i])[9:16,9:16]))

feature[i,16]<-sum((get(data1[i])[17:24,9:16]))

feature[i,17]<-sum((get(data1[i])[25:32,9:16]))

feature[i,18]<-sum((get(data1[i])[1:8,17:24]))

feature[i,19]<-sum((get(data1[i])[9:16,17:24]))

feature[i,20]<-sum((get(data1[i])[17:24,17:24]))

feature[i,21]<-sum((get(data1[i])[25:32,17:24]))

feature[i,22]<-sum((get(data1[i])[1:8,25:32]))

feature[i,23]<-sum((get(data1[i])[9:16,25:32]))

feature[i,24]<-sum((get(data1[i])[17:24,25:32]))

feature[i,25]<-sum((get(data1[i])[25:32,25:32]))

}

labeltest<-factor(rep(0:9,c(87,97,92,85,114,108,87,96,91,89)))

data2<-data.frame(feature,labeltest)

pred1<-predict(m1,data2)

table(pred1,labeltest)

sum(diag(table(pred1,labeltest)))/length(name)

經整理，輸出結果如下：(左邊為訓練集，右邊為測試集)

Tips: oddsratio=p/1-p 相對風險指數貝努力模型中 P是發生A事件的概率，1-p是不發生A事件的概率所以p/1-p是 發生與不發生的相對風險。OR值等於1，表示該因素對A事件發生不起作用;OR值大於1，表示A事件發生的可能性大於不發生的可能性;OR值小於1，表示A事件不發生的可能性大於發生的可能性。

Further reading：

Yves Croissant：Estimation of multinomial logit models in R : The mlogit Packages