R語言glm函數學習：
【轉載時請注明來源】：http://www.cnblogs.com/runner-ljt/
Ljt
作為一個初學者，水平有限，歡迎交流指正。
glm函數介紹：
glm(formula, family=family.generator, data,control = list(...))
family：每一種響應分布（指數分布族）允許各種關聯函數將均值和線性預測器關聯起來。
常用的family：
binomal(link='logit') ----響應變量服從二項分布，連接函數為logit，即logistic回歸
binomal(link='probit') ----響應變量服從二項分布，連接函數為probit
poisson(link='identity') ----響應變量服從泊松分布，即泊松回歸
control:控制算法誤差和最大迭代次數
glm.control(epsilon = 1e-8, maxit = 25, trace = FALSE)
-----maxit:算法最大迭代次數，改變最大迭代次數：control=list(maxit=100)
glm函數使用：
>
> data<-iris[1:100,]
> samp<-sample(100,80)
> names(data)<-c('sl','sw','pl','pw','species')
> testdata<-data[samp,]
> traindata<-data[-samp,]
>
> lgst<-glm(testdata$species~pl,binomial(link='logit'),data=testdata)
Warning messages:
1: glm.fit:算法沒有聚合
2: glm.fit:擬合機率算出來是數值零或一
> summary(lgst)
Call:
glm(formula = testdata$species ~ pl, family = binomial(link = "logit"),
data = testdata)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.836e-05 -2.110e-08 -2.110e-08 2.110e-08 1.915e-05
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -83.47 88795.25 -0.001 0.999
pl 32.09 32635.99 0.001 0.999
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1.1085e+02 on 79 degrees of freedom
Residual deviance: 1.4102e-09 on 78 degrees of freedom
AIC: 4
Number of Fisher Scoring iterations: 25
>
注意在使用glm函數就行logistic回歸時，出現警告：
Warning messages:
1: glm.fit:算法沒有聚合
2: glm.fit:擬合機率算出來是數值零或一
同時也可以發現兩個系數的P值都為0.999，說明回歸系數不顯著。
第一個警告：算法不收斂。
由於在進行logistic回歸時，依照極大似然估計原則進行迭代求解回歸系數，glm函數默認的最大迭代次數 maxit=25，當數據不太好時，經過25次迭代可能算法 還不收斂，所以可以通過增大迭代次數嘗試解決算法不收斂的問題。但是當增大迭代次數后算法仍然不收斂，此時數據就是真的不好了，需要對數據進行奇異值檢驗等進一步的處理。
>
> lgst<-glm(testdata$species~pl,binomial(link='logit'),data=testdata,control=list(maxit=100))
Warning message:
glm.fit:擬合機率算出來是數值零或一
> summary(lgst)
Call:
glm(formula = testdata$species ~ pl, family = binomial(link = "logit"),
data = testdata, control = list(maxit = 100))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.114e-05 -2.110e-08 -2.110e-08 2.110e-08 1.162e-05
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -87.18 146399.32 -0.001 1
pl 33.52 53808.49 0.001 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1.1085e+02 on 79 degrees of freedom
Residual deviance: 5.1817e-10 on 78 degrees of freedom
AIC: 4
Number of Fisher Scoring iterations: 26
>
如上，通過增加迭代次數，解決了第一個警告，此時算法收斂。
但是第二個警告仍然存在，且回歸系數P=1，仍然不顯著。
第二個警告：擬合概率算出來的概率為0或1
首先，這個警告是什么意思？
我們先來看看訓練樣本的logist回歸結果，擬合出的每個樣本屬於'setosa'類的概率為多少？
> >lgst<-glm(testdata$species~pl,binomial(link='logit'),data=testdata,control=list(maxit=100)) >p<-predict(lgst,type='response') >plot(seq(-2,2,length=80),sort(p),col='blue') >
可以看出訓練樣本為'setosa'類的概率不是幾乎為0，就是幾乎為1，並不是我們預想中的logistic模型的S型曲線，這就是第二個警告的意思。
那么問題來了，為什么會出現這種情況？
（以下內容只是本人參考一些解釋的個人理解）
這種情況的出現可以理解為一種過擬合，由於數據的原因，在回歸系數的優化搜索過程中，使得分類的種類屬於某一種類（y=1)的線性擬合值趨於大，分類種類為另一 類(y=0)的線性擬合值趨於小。
由於在求解回歸系數時，使用的是極大似然估計的原理，即回歸系數在搜索過程中使得似然函數極大化：
所以在搜索過程中偏向於使得y=1的h(x)趨向於大，而使得y=0的h(x)趨向於小。
即系數Θ使得 Y=1類的 -ΘTX 趨向於大，使得Y=0類的 -ΘTX 趨向於小。而這樣的結果就會導致P(y=1|x;Θ)-->1 ； P(y=0|x;Θ)-->0 .
那么問題又來了，什么樣的數據會導致這樣的過擬合產生呢？
先來看看上述logistic回歸中種類為setosa和versicolor的樣本pl值的情況。（橫軸代表pl值，為了避免樣本pl數據點疊加在一起，增加了一個無關的y值使樣本點展開）
可以看出兩類數據明顯的完全線性可分。
故在回歸系數搜索過程中只要使得一元線性函數h(x)的斜率的絕對值偏大，就可以實現y=1類的h(x)趨向大，y=0類的h(x)趨向小。
所以當樣本數據完全可分時，logistic回歸往往會導致過擬合的問題，即出現第二個警告：擬合概率算出來的概率為0或1。
出現了第二個警告后的logistic模型進行預測時往往是不適用的，對於這種線性可分的樣本數據，其實直接使用規則判斷的方法則簡單且適用（如當pl<2.5時則直接判斷為setosa類，pl>2.5時判斷為versicolor類）。
以下，對於不完全可分的二維訓練數據展示logistic回歸過程。
>
> data<-iris[51:150,]
> samp<-sample(100,80)
> names(data)<-c('sl','sw','pl','pw','species')
> testdata<-data[samp,]
> traindata<-data[-samp,]
>
> lgst<-glm(testdata$species~sw+pw,binomial(link='logit'),data=testdata)
> summary(lgst)
Call:
glm(formula = testdata$species ~ sw + pw, family = binomial(link = "logit"),
data = testdata)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.82733 -0.16423 0.00429 0.11512 2.12846
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -12.915 5.021 -2.572 0.0101 *
sw -3.796 1.760 -2.156 0.0310 *
pw 14.735 3.642 4.046 5.21e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 110.85 on 79 degrees of freedom
Residual deviance: 24.40 on 77 degrees of freedom
AIC: 30.4
Number of Fisher Scoring iterations: 7
>#畫擬合概率曲線圖
> p<-predict(lgst,type='response')
> plot(seq(-2,2,length=80),sort(p),col='blue')
>
>#畫訓練樣本數據散點圖
>a<-testdata$species=='versicolor'
> x1<-testdata[a,'sw']
> y1<-testdata[a,'pw']
> x2<-testdata[!a,'sw']
> y2<-testdata[!a,'pw']
> summary(testdata$sw)
Min. 1st Qu. Median Mean 3rd Qu. Max.
2.000 2.700 2.900 2.881 3.100 3.800
> summary(testdata$pw)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 1.300 1.600 1.672 2.000 2.500
>
> plot(x1,y1,xlim=c(1.5,4),ylim=c(.05,3),xlab='sw',ylab='pw',col='blue')
> points(x2,y2,col='red')
>
> #畫分類邊界圖，即畫h(x)=0.5的圖像
> x3<-seq(1.5,4,length=100)
> y3<-(3.796/14.735)*x3+13.415/14.735
> lines(x3,y3)
擬合概率曲線圖：
（基本上符合logistic模型的S型曲線）
訓練樣本散點圖及分類邊界：
（畫logistic回歸的分類邊界即畫曲線h(x)=0.5）
