R語言與機(jī)器學(xué)習(xí)logistic回歸

1、寫在前面的廢話2014,又到了新的一年，首先祝大家新年快樂，也感謝那些關(guān)注我的博客的人現(xiàn)在想想數(shù)據(jù)挖掘課程都是去年的事了，一直預(yù)告著，盤算著年內(nèi)完工的分類算法也拖了一年了。本來打算去年就完成分類算法，如果有人看的話也順帶提提關(guān)聯(lián)分析，聚類神馬的，借著新年新氣象的借口來補(bǔ)完這一系列的文章,可是，這明明就是在發(fā)盡管這個是預(yù)告里的最后一篇，但是我也沒打算把這個分類算法就這么完結(jié)。盡管每一篇都很淺顯，每個算法都是淺嘗輒止的，在deeplearning那么火的今天，掌握這些東西算起來展絲得不能再展絲了?？紤]到一致性與完備性， 最后補(bǔ)上兩篇一樣naive的： 組合方法提高分類效率、 幾種分類方法的績效討論

2、。希望讀到的人喜歡。算法六：logistic回歸由于我們在前面已經(jīng)討論過了神經(jīng)網(wǎng)絡(luò)的分類問題（參見R語言與機(jī)器學(xué)習(xí)學(xué)習(xí)筆記（分類算法）（5）,如今再從最優(yōu)化的角度來討論logistic回歸就顯得有些不合適了。Logistic回歸問題的最優(yōu)化問題可以表述為：尋找一個非線性函數(shù)sigmoid的最佳擬合參數(shù)，求解過程可使用最優(yōu)化算法完成。它可以看做是用sigmoid函數(shù)作為二閾值分類器的感知器問題。今天我們將從統(tǒng)計的角度來重新考慮logistic回歸問題。一、logistic回歸及其MLE當(dāng)我們考慮解釋變量為分類變量如考慮一個企業(yè)是否會被并購，一個企業(yè)是否會上市，你的能否考上研究生這些問題時，考慮線

3、性概率模型P（yi=1）=00+B1xi顯然是不合適的，它至少有兩個致命的缺陷：1、概率估計值可能超過1,使得模型失去了意義；（要解決這個問題并不麻煩，我們將預(yù)測超過1的部分記為1,低于0的部分記為0,就可以解決。這個解決辦法就是計量里有一定歷史的tobit模型）2、邊際效應(yīng)假定為不變，通常來說不合經(jīng)濟(jì)學(xué)常識?？紤]一個邊際效應(yīng)遞減的模型（假定真實值為藍(lán)線），可以看到線性模型表現(xiàn)很差。但是sigmoid函數(shù)去擬合藍(lán)線確實十分合適的。于是我們可以考慮logistic回歸模型：f白KP-0+凡士)一J)-1+門即R+再，)假定有N個觀測樣本Y1,Y2,YN,設(shè)P(Yi=1|Xi)=兀(Xi)為給定條

4、件Xi下得到結(jié)果Yi=1的條件概率；而在同樣條件下得到結(jié)果Yi=0的條件概率為P(Yi=0|Xi)=1-兀(Xi),于是得到一個觀測值的概率P(Yi)=兀(Xi)Yi1-兀(Xi)1-Yi假設(shè)各觀測獨(dú)立，對logistic回歸模型來說，其對數(shù)似然函數(shù)為：二皿(/?)KSI匕加帝(/)4(1-I-T(黑i I于是便可求解出logistic模型的MLE二、logit還是probit?雖說sigmoid函數(shù)對邊際遞減的模型擬合良好，但是我們也要知道S型函數(shù)并非僅sigmoid函數(shù)一個，絕大多數(shù)的累積分布函數(shù)都是S型的。于是考慮F-1(P)(F為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù))也不失為一個很好的選擇。像這

5、樣的，對概率P做一點(diǎn)變換，讓變換后的取值范圍變得合理，且變換后我們能夠有辦法進(jìn)行參數(shù)估計的，就涉及到廣義線性模型理論中的連接函數(shù)。在廣義線性模型中我們把10g(P/(1-P)稱為logit,F-1(P)(F為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù))稱為probit。那么這里就涉及到一個選擇的問題：連接函數(shù)選logit還是probit?logistic回歸認(rèn)為二分類變量服從伯努利分布，應(yīng)當(dāng)選擇logit,而且從解釋的角度說，p/(1-p)就是我們常說的oddsratio,也就是軟件報告中出現(xiàn)的ORfi。但是probit也有它合理的一面，首先，中心極限定理告訴我們，伯努利分布在樣本夠多的時候就是近似正態(tài)分布的

6、；其次，從不確定性的角度考慮，probit認(rèn)為我們的線性概率模型服從正態(tài)分布，這也是更為合理的。我們來看一下經(jīng)過變換后，自變量和P的關(guān)系是什么樣子的：如果你確實想知道到底你的數(shù)據(jù)用哪一個方法好，也不是沒有辦法，你可以看下你的殘差到底是符合logit函數(shù)呢還是符合probit函數(shù)，當(dāng)然,憑肉眼肯定是看不出來的，因為這兩個函數(shù)本來就很接近，你可以通probit過函數(shù)的假定，用擬合優(yōu)度檢驗一下。但通常，估計不會有人非要這么較真，因為沒有必要。但是有一點(diǎn)是要注意的，logit模型較probit模型而言具有厚尾的特征，這也是為什么經(jīng)濟(jì)學(xué)論文愛用logit的原因。表1 1FrFr MtMt根型和1 1唯i

7、lil模用概率值累積止態(tài)分布函敷累和加gisticgistic概率分布函數(shù)yi1 1Pi 3 30.W130.W13-1.0-1.0.(I22S.(I22S-1.5-1.5U.0MU.0M強(qiáng)-LD-LD0.15M70.15M7禮2 2礎(chǔ)9 9-0.5-0.50*37750*3775也D DQ/000Q/0000 0 3 3徜町仇他1515加5 5140.0. 4134130.73110.73111.5口用出OI7AOI7A2J)2J)D.M772D.M772“LXLX戕削3.(13.(1色9 9第7 7我們以鶯尾花數(shù)據(jù)中的virginica,versicolor法分類有無差別。兩類數(shù)據(jù)分類為例

8、，看看兩種辦 9 9i1i14tt*上*務(wù)*n* *一一一一 T76-54376-543L-t* itk 二二工工 314仁式電中仁式電中( (% %+ +即)與logistic回歸的似然估計類似，我們可以很容易寫出多項logit的對數(shù)似然函數(shù)：L(fj=n(L(p)=22尸11j-1其中可取其中可取 1 當(dāng)且僅豈方案當(dāng)且僅豈方案 j 襁選中襁選中. .否貝否貝 取取 0.1 的概率可由的概率可由 ?=f)=f)=1 解出解出,+,+多項Logit模型雖然好用，但從上面的敘述可以看出，多項Logit模型最大的限制在于各個類別必須是對等的，因此在可供選擇的類別中，不可有主要類別和次要類別混雜在一

9、起的情形。例如在研究旅游交通工具的選擇時，可將交通工具的類別粗分為航空、火車、公用汽車、自用汽車四大類，但若將航空類別再依三家航空公司細(xì)分出三類而得到總共六個類別，則多項Logit模型就不適用，因為航空、火車、公用汽車、自用汽車均屬同一等級的主要類別，而航空公司的區(qū)別則很明顯的是較次要的類別，不應(yīng)該混雜在一起。在這個例子中，主要類別和次要類別很容易分辨，但在其他的研究中可能就不是那么容易，若不慎將不同層級的類別混在一起,則由多項Logit模型所得到的實證結(jié)果就會有誤差。對于分類模型，我們還會遇到被解釋變量中有分類變量的情形。對于連續(xù)變量解釋離散變量，且被解釋的離散變量是有順序的（這個是和多項l

10、ogit最大的區(qū)別）的情形，我們就需要考慮到orderlogit模型。其數(shù)學(xué)模型敘述如下：P P（i i2 2I I國，力=Fg ffTi=ff| |天5 5A A）= =1-1-百（左丁氏河其中，F(xiàn)（.）表示累積分布函數(shù)，當(dāng)F表示正態(tài)分布的分布函數(shù)時，對應(yīng)的是orderprobit;F表示logistic分布時，變對應(yīng)orderlogit。與logistic分布類似，我們可以很容易寫出其對數(shù)似然函數(shù)：N N冊冊m駐駐 外外= =ZZ/ /。耳。耳（砥為（砥為= =）區(qū)）區(qū), ,A力力 。必。必= =/ /）其中D D司）是指示函數(shù),當(dāng)j j=12=12才m m為真時，如釬/）=L L當(dāng)J J

11、閆,為不真時,22弓尸（k k對于樣本中某個M M來說，機(jī)個可能值，同時只有一個為真口出駐 3方=工圖=l)-Hl*og(7V=E)以陽)1=1四、dummyvariable在logistic回歸中，經(jīng)常會遇到解釋變量為分類變量的情形，比如收入：高、中、低；地域：北京、上海、廣州等。這里對分類變量而言就涉及一個問題：要不要將分類變量設(shè)置dummyvariable?這個問題的答案在線性模型中很顯然，必須要這么做!如果我們不設(shè)置啞變量,而是單純地賦值：北京=1,上海=2,廣州=3,即我們將自變量視作連續(xù)性的數(shù)值變量，但這僅僅是個代碼而己，并不意味著地域問存在大小次序的關(guān)系，即并非代表被解釋變量(響

12、應(yīng)變量)會按此順序線性增加或減少。即使是有序多分類變量，如家庭收入分為高、中、低三檔，各類別間的差距也是無法準(zhǔn)確衡量的，按編碼數(shù)值來分析實際上就是強(qiáng)行規(guī)定為等距，這顯然可能引起更大的誤差。但是在logistic回歸中， 由于logit(p)變化的特殊性， 在解釋定序變量時， 為了減少自由度(即解釋變量個數(shù))，我們常常將定序變量(如家庭收入分為高、中、低)視為連續(xù)的數(shù)值變量，而且經(jīng)濟(jì)解釋可以是XX每提高一個檔次，相應(yīng)的概率會提高expression(delta(XX)(expression的表達(dá)式還是很復(fù)雜的，不打了)。當(dāng)然減少變量個數(shù)是以犧牲預(yù)測精度為代價的。畢竟數(shù)據(jù)處理是一門藝術(shù)而非一門技術(shù)

13、，如何取舍還得具體問題具體分析。當(dāng)然，非定序的分類變量是萬萬不可將其視為數(shù)值變量的。五、廣義線性模型的R實現(xiàn)R語言提供了廣義線性模型的擬合函數(shù)glm(),其調(diào)用格式如下：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,.)參數(shù)說明：Formula:回歸形式，與lm()函數(shù)的formula參數(shù)用法一致Family：設(shè)置

14、廣義線性模型連接函數(shù)的典則分布族，glm()提供正態(tài)、指數(shù)、gamma逆高斯、Poisson、二項分布。我們的logistic回歸使用的是二項分布族binomial0Binomial族默認(rèn)連接函數(shù)為logit,可設(shè)置為probit。Data：數(shù)據(jù)集鶯尾花例子使用的R代碼：logit.fit-glm(SpeciesPetal.Width+Petal.Length,family=binomial(link=logit),data=iris51:150,)logit.predictions0,virginica,versicolor)table(iris51:150,5,logit.predicti

15、ons)probit.fit-glm(SpeciesPetal.Width+Petal.Length,family=quasibinomial(link=probit),data=iris51:150,)probit.predictions0,virginica,versicolor)table(iris51:150,5,probit.predictions)程序包mlogit提供了多項logit的模型擬合函數(shù):mlogit(formula,data,subset,weights,na.action,start=NULL,alt.subset=NULL,reflevel=NULL,nests=

16、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.ind

17、ex=FALSE,ranked=FALSE,.)參數(shù)說明：formula:mlogit提供了條件logit,多項logit,混合logit多種模型，對于多項logit的估計模型應(yīng)寫為：因變量0|自變量，如：mode0|incomedata:使用mlogit.data函數(shù)使得數(shù)據(jù)結(jié)構(gòu)符合mlogit函數(shù)要求。Choice:確定分類變量是什么Shape：如果每一行是一個觀測，我們選擇wide,如果每一行是表示一個選擇，那么就應(yīng)該選擇longoalt.var:對于shape為10ng的數(shù)據(jù)，需要標(biāo)明所有的選擇名稱選才wwide的數(shù)據(jù)示例：price,pierprice.boatprice,chart

18、ercat157.930157.930182.93015.11410.53434.534161.87424.33459.33415.13455.93084.930106.93041.51471.014192.47428,93463.9341 .chartercharterl.pierpier2 .beachbeachhead(Fishnng)modeprice.beach1charter157.9302charter15.1143boat161.8744pier15.1345boat106.9306charter192.474catch.piercatch.boatcatch.charteri

19、ncone0.05030.26010.04510.15740.45220.24130.07890.16430.539170833320,46711250.0001.02663750.0000.53912083.3330.05030.10820.45220.16650.32404583.3320.39754583.332FTshhead(Fish)modeprice-beach1.beachFALSE157.9301.boatFALSE157.9301,charterTRUE157.930l.pierFALSE157.9302 .beachFALSE15.1142.boatFALSE15.114

20、price,pierprice,boatprice,chart157.930157.930182.9157,930157.930182.9157.930157.930182*9157,930157.930182.915.11410.53434.5catch.beach1.beach1.boat1 .charter1,pier2 .beich0.06780.06780.06780.06780.1049catch.pler0.0503O.O5O30.05030.05030.0451catch.boatcatch-charter10-26010.26010.26010.26010.15740.539

21、17080.53917080.53917080.53917080.46711251.beach1,boat0.1049altbeachbeat0.04510.15740.4671125選才?long的數(shù)據(jù)示例:以巾shing數(shù)據(jù)為例，來說明如何使用mlogitlibrary(mlogit)data(Fishing,package=mlogit)Fish TMv-mlogit.data(Trave1ModeTchoi十a(chǎn)lt.var=“ode”data(TrvelMode,package=AER)head(TravelMode)h史ad(TM)toto1 1OOOO8484 s s7773687

22、77368- -o o. .c c*_*_9-10270849-1027084一eo71865eo71865V13413V13413incomesize351351351351302302,shape=1individualmodechoicewait1. air1 .trainl.bus1.carairFAL5E69trainFALSE34busFALSE35carTRUE02.ai2.trainairFALSE64trainFALSE44t t915915O O8181S S532153532153O OU Utravel10037241718068354gcostincorr703713

23、703303683843包的multinom作介紹了，可以參見根據(jù)EconometricsinR一書，將回歸方法總結(jié)一下一文。程序包MASS1供polr()函數(shù)可以進(jìn)行orderedlogit或probit回歸。用法如下：polr(formula,data,weights,start,.,subset,na.action,contrasts=NULL,Hess=FALSE,model=TRUE,method=c(logistic,probit,cloglog,cauchit)參數(shù)說明：Formula:回歸形式，與lm()函數(shù)的formula參數(shù)用法一致Data：數(shù)據(jù)集Method:默認(rèn)為ord

24、erlogit,選擇probit時變?yōu)閛rderprobit模型。以housing數(shù)據(jù)為例說明函數(shù)用法：house.plr-polr(SatInfl+Type+Cont,weights=Freq,data=housing)house.plrsummary(house.plr,digits=3)這些結(jié)果十分直觀，易于解讀，所以我們在這里省略所有的輸出結(jié)果。再看手寫數(shù)字案例：最后，我們回到最開始的那個手寫數(shù)字的案例，我們試著利用多項logit重做這個案例。(這個案例的描述與數(shù)據(jù)參見kNN算法一章)特征的選擇可參見神經(jīng)網(wǎng)絡(luò)一章。由于手寫數(shù)字的特征選取很容易導(dǎo)致回歸系數(shù)矩陣是降秩的，所以我們使用nne

25、t包的multinom()函數(shù)代替mlogit()。運(yùn)行下列代碼：setwd(D:/R/data/digits/trainingDigits)names-list.files(D:/R/data/digits/trainingDigits)data-paste(train,1:1934,sep=)for(iin1:length(names)assign(datai,as.matrix(read.fwf(namesi,widths=rep(1,32)label-factor(rep(0:9,c(189,198,195,199,186,187,195,201,180,204)feature-mat

26、rix(rep(0,length(names)*25),length(names),25)for(iin1:length(names)featurei,1-sum(get(datai),16)featurei,2-sum(get(datai),8)featurei,3-sum(get(datai),24)featurei,4-sum(get(datai)16,)featurei,5-sum(get(datai)11,)featurei,6-sum(get(datai)21,)featurei,7-sum(diag(get(datai)featurei,8-sum(diag(get(datai)

27、,32:1)featurei,9-sum(get(datai)17:32,17:32)featurei,10-sum(get(datai)1:8,1:8)featurei,11-sum(get(datai)9:16,1:8)featurei,12-sum(get(datai)17:24,1:8)featurei,13-sum(get(datai)25:32,1:8)featurei,14-sum(get(datai)1:8,9:16)featurei,15-sum(get(datai)9:16,9:16)featurei,16-sum(get(datai)17:24,9:16)featurei

28、,17-sum(get(datai)25:32,9:16)featurei,18-sum(get(datai)1:8,17:24)featurei,20-sum(get(datai)17:24,17:24)featurei,21-sum(get(datai)25:32,17:24)featurei,22-sum(get(datai)1:8,25:32)featurei,23-sum(get(datai)9:16,25:32)featurei,24-sum(get(datai)17:24,25:32)featurei,25-sum(get(datai)25:32,25:32)data1-data

29、.frame(feature,label)#降秩時mlogit不可用#data10-mlogit.data(data1,shape=wide,choice=label)#m1-mlogit(label0|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(tabl

30、e(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(iin1:length(name)assign(data1i,as.matrix(read.fwf(namei,widths=rep(1,32)feature-matrix(rep(0,length(name)*25),length(name),25)for(iin1:length(name)featurei,1-sum

31、(get(data1i),16)featurei,2-sum(get(data1i),8)featurei,3-sum(get(data1i),24)featurei,4-sum(get(data1i)16,)featurei,5-sum(get(data1i)11,)featurei,6-sum(get(data1i)21,)featurei,7-sum(diag(get(data1i)featurei,8-sum(diag(get(data1i),32:1)featurei,9-sum(get(data1i)17:32,17:32)featurei,10-sum(get(data1i)1:8,1:8)featurei,11-sum(get(data1i)9:16,1:8)featurei,12-sum(get(data1i)17:24,1:8)featurei,13-sum(get(data1i)25:32,1:8)featurei,14-sum(get(data1i)1:8,9:16)featurei,15-sum(get(data1