|
在日常學習或工作中經(jīng)常會使用線性回歸模型對某一事物進行預測,，例如預測房價、身高,、GDP,、學生成績等，發(fā)現(xiàn)這些被預測的變量都屬于連續(xù)型變量,。然而有些情況下,，被預測變量可能是二元變量，即成功或失敗,、流失或不流失,、漲或跌等,，對于這類問題，線性回歸將束手無策,。這個時候就需要另一種回歸方法進行預測,，即Logistic回歸。 在實際應用中,，Logistic模型主要有三大用途： 1）尋找危險因素,，找到某些影響因變量的"壞因素"，一般可以通過優(yōu)勢比發(fā)現(xiàn)危險因素,； 2）用于預測,，可以預測某種情況發(fā)生的概率或可能性大小,； 3）用于判別,，判斷某個新樣本所屬的類別。 Logistic模型實際上是一種回歸模型,，但這種模型又與普通的線性回歸模型又有一定的區(qū)別： 1）Logistic回歸模型的因變量為二分類變量,； 2）該模型的因變量和自變量之間不存在線性關系； 3）一般線性回歸模型中需要假設獨立同分布,、方差齊性等,，而Logistic回歸模型不需要； 4）Logistic回歸沒有關于自變量分布的假設條件,，可以是連續(xù)變量,、離散變量和虛擬變量； 5）由于因變量和自變量之間不存在線性關系,，所以參數(shù)(偏回歸系數(shù))使用最大似然估計法計算,。 logistic回歸模型概述 廣義線性回歸是探索“響應變量的期望”與“自變量”的關系，以實現(xiàn)對非線性關系的某種擬合,。這里面涉及到一個“連接函數(shù)”和一個“誤差函數(shù)”,，“響應變量的期望”經(jīng)過連接函數(shù)作用后，與“自變量”存在線性關系,。選取不同的“連接函數(shù)”與“誤差函數(shù)”可以構(gòu)造不同的廣義回歸模型,。當誤差函數(shù)取“二項分布”而連接函數(shù)取“l(fā)ogit函數(shù)”時，就是常見的“l(fā)ogistic回歸模型”,，在0-1響應的問題中得到了大量的應用,。 Logistic回歸主要通過構(gòu)造一個重要的指標：發(fā)生比來判定因變量的類別。在這里我們引入概率的概念,，把事件發(fā)生定義為Y=1,，事件未發(fā)生定義為Y=0，那么事件發(fā)生的概率為p,事件未發(fā)生的概率為1-p,，把p看成x的線性函數(shù),； 回歸中,，最常用的估計是最小二乘估計，因為使得p在[0,1]之間變換,，最小二乘估計不太合適,，有木有一種估計法能讓p在趨近與0和1的時候變換緩慢一些（不敏感），這種變換是我們想要的,，于是引入Logit變換,對p/(1-p)也就是發(fā)生與不發(fā)生的比值取對數(shù),，也稱對數(shù)差異比。經(jīng)過變換后,，p對x就不是線性關系了,。 logistic回歸的公式可以表示為： 其中P是響應變量取1的概率，在0-1變量的情形中,，這個概率就等于響應變量的期望,。 這個公式也可以寫成： 可以看出,，logistic回歸是對0-1響應變量的期望做logit變換,，然后與自變量做線性回歸。參數(shù)估計采用極大似然估計,，顯著性檢驗采用似然比檢驗,。 建立模型并根據(jù)AIC準則選擇模型后，可以對未知數(shù)據(jù)集進行預測,，從而實現(xiàn)分類,。模型預測的結(jié)果是得到每一個樣本的響應變量取1的概率，為了得到分類結(jié)果,，需要設定一個閾值p0——當p大于p0時,，認為該樣本的響應變量為1，否則為0,。閾值大小對模型的預測效果有較大影響,，需要進一步考慮。首先必須明確模型預測效果的評價指標,。 對于0-1變量的二分類問題,，分類的最終結(jié)果可以用表格表示為： 其中，d是“實際為1而預測為1”的樣本個數(shù),，c是“實際為1而預測為0”的樣本個數(shù),，其余依此類推。 顯然地,，主對角線所占的比重越大,，則預測效果越佳，這也是一個基本的評價指標——總體準確率(a+d)/(a+b+c+d),。 準確（分類）率=正確預測的正反例數(shù)/總數(shù) Accuracy=(a+d)/(a+b+c+d) 誤分類率=錯誤預測的正反例數(shù)/總數(shù) Error rate=(b+c)/(a+b+c+d)=1-Accuracy 正例的覆蓋率=正確預測到的正例數(shù)/實際正例總數(shù) Recall(True Positive Rate,，or Sensitivity)=d/(c+d) 正例的命中率=正確預測到的正例數(shù)/預測正例總數(shù) Precision(Positive Predicted Value,PV+)=d/(b+d) 負例的命中率=正確預測到的負例個數(shù)/預測負例總數(shù) Negative predicted value(PV-)=a/(a+c) 通常將上述矩陣稱為“分類矩陣”,。一般情況下，我們比較關注響應變量取1的情形,，將其稱為Positive（正例）,，而將響應變量取0的情形稱為Negative（負例）。常見的例子包括生物實驗的響應,、營銷推廣的響應以及信用評分中的違約等等,。針對不同的問題與目的，我們通常采用ROC曲線與lift曲線作為評價logistic回歸模型的指標,。 1）ROC曲線 設置了兩個相應的指標：TPR與FPR,。 TPR：True Positive Rate（正例覆蓋率），將實際的1正確地預測為1的概率,，d/(c+d),。 FPR：False Positive Rate，將實際的0錯誤地預測為1的概率,，b/(a+b),。 TPR也稱為Sensitivity（即生物統(tǒng)計學中的敏感度），也可以稱為“正例的覆蓋率”——將實際為1的樣本數(shù)找出來的概率,。覆蓋率是重要的指標,，例如若分類的目標是找出潛在的劣質(zhì)客戶（響應變量取值為1），則覆蓋率越大表示越多的劣質(zhì)客戶被找出,。 類似地,，1-FPR其實就是“負例的覆蓋率”，也就是把負例正確地識別為負例的概率,。 TPR與FPR相互影響,，而我們希望能夠使TPR盡量地大，而FPR盡量地小,。影響TPR與FPR的重要因素就是上文提到的“閾值”,。當閾值為0時，所有的樣本都被預測為正例,，因此TPR=1,，而FPR=1。此時的FPR過大,，無法實現(xiàn)分類的效果,。隨著閾值逐漸增大，被預測為正例的樣本數(shù)逐漸減少,，TPR和FPR各自減小,，當閾值增大至1時，沒有樣本被預測為正例，此時TPR=0,，F(xiàn)PR=0,。 由上述變化過程可以看出，TPR與FPR存在同方向變化的關系（這種關系一般是非線性的）,，即,，為了提升TPR（通過降低閾值），意味著FPR也將得到提升,，兩者之間存在類似相互制約的關系,。我們希望能夠在犧牲較少FPR的基礎上盡可能地提高TPR，由此畫出了ROC曲線,。 ROC曲線的全稱為“接受者操作特性曲線”（receiver operating characteristic）,，其基本形式為： 當預測效果較好時，ROC曲線凸向左上角的頂點,。平移圖中對角線,，與ROC曲線相切，可以得到TPR較大而FPR較小的點,。模型效果越好,，則ROC曲線越遠離對角線，極端的情形是ROC曲線經(jīng)過（0,，1）點,，即將正例全部預測為正例而將負例全部預測為負例,。ROC曲線下的面積可以定量地評價模型的效果,，記作AUC，AUC越大則模型效果越好,。 當我們分類的目標是將正例識別出來時（例如識別有違約傾向的信用卡客戶）,，我們關注TPR，此時ROC曲線是評價模型效果的準繩,。 2）lift曲線 在營銷推廣活動中,，我們的首要目標并不是盡可能多地找出那些潛在客戶，而是提高客戶的響應率,?？蛻繇憫适怯绊懲度氘a(chǎn)出比的重要因素。此時,，我們關注的不再是TPR（覆蓋率）,，而是另一個指標：命中率。 回顧前面介紹的分類矩陣,，正例的命中率是指預測為正例的樣本中的真實正例的比例,，即d/(b+d)，一般記作PV。 在不使用模型的情況下,，我們用先驗概率估計正例的比例,，即(c+d)/(a+b+c+d)，可以記為k,。 定義提升值lift=PV/k,。 lift揭示了logistic模型的效果。例如,，若經(jīng)驗告訴我們10000個消費者中有1000個是我們的潛在客戶,，則我們向這10000個消費者發(fā)放傳單的效率是10%（即客戶的響應率是10%），k=(c+d)/(a+b+c+d)=10%,。通過對這10000個消費者進行研究,，建立logistic回歸模型進行分類，我們得到有可能比較積極的1000個消費者,，b+d=1000,。如果此時這1000個消費者中有300個是我們的潛在客戶，d=300,，則命中率PV為30%,。此時，我們的提升值lift=30%/10%=3,，客戶的響應率提升至原先的三倍,，提高了投入產(chǎn)出比。 為了畫lift圖,，需要定義一個新的概念depth深度,，這是預測為正例的比例，(b+d)/(a+b+c+d),。 與ROC曲線中的TPR和FPR相同,，lift和depth也都受到閾值的影響。 當閾值為0時,，所有的樣本都被預測為正例,，因此depth=1，而PV=d/(b+d)=(0+d)/(0+b+0+d)=k,，于是lift=1,，模型未起提升作用。隨著閾值逐漸增大,，被預測為正例的樣本數(shù)逐漸減少,，depth減小，而較少的預測正例樣本中的真實正例比例逐漸增大,。當閾值增大至1時,，沒有樣本被預測為正例,，此時depth=0，而lift=0/0,。 由此可見,，lift與depth存在相反方向變化的關系。在此基礎上作出lift圖：
與ROC曲線不同,，lift曲線凸向（0,，1）點。我們希望在盡量大的depth下得到盡量大的lift（當然要大于1）,，也就是說這條曲線的右半部分應該盡量陡峭,。 至此，我們對ROC曲線和lift曲線進行了描述,。這兩個指標都能夠評價logistic回歸模型的效果,，只是分別適用于不同的問題： 如果是類似信用評分的問題，希望能夠盡可能完全地識別出那些有違約風險的客戶（不使一人漏網(wǎng)）,，我們需要考慮盡量增大TPR（覆蓋率）,，同時減小FPR（減少誤殺），因此選擇ROC曲線及相應的AUC作為指標,； 如果是做類似數(shù)據(jù)庫精確營銷的項目,，希望能夠通過對全體消費者的分類而得到具有較高響應率的客戶群，從而提高投入產(chǎn)出比,，我們需要考慮盡量提高lift（提升度）,，同時depth不能太小（如果只給一個消費者發(fā)放傳單,，雖然響應率較大,，卻無法得到足夠多的響應），因此選擇lift曲線作為指標,。 相關R應用包 普通二分類 logistic 回歸 用系統(tǒng)的 glm 因變量多分類 logistic 回歸 有序分類因變量：用 MASS 包里的 polrb 無序分類因變量：用 nnet 包里的 multinom 條件logistic回歸,，用 survival 包里的 clogit R運行邏輯回歸 R可以讓邏輯回歸建模過程變得很簡單。我們可以使用glm()函數(shù)進行邏輯回歸建模,，而且，它的擬合過程和我們之前所做的線性回歸沒有很大的區(qū)別,。在這篇博客中,，我們將介紹如何一步一步的進行邏輯回歸建模。 我們用r來看一下曲線的形狀： x <- seq(from = -10, to = 10, by = 0.01)
p在0和1附近變化不敏感,，ps:這條曲線幫了我一個大忙,，最近做RTB競價算法時，需要按照CTR設定一個CTR指導價,，然后以CPC做曲線的修正,，出于成本考慮，我們的出價不會高于某個固定值K，所以把這條曲線變換一下就派上大用場了,。（扯遠了） 這條曲線的主要特性是θ函數(shù)取值可以在-∞變到+∞,，p的取值在[0,1],且極值端變化不敏感。這就是我們想要的,。值得注意的是,，經(jīng)過logit變換，已經(jīng)不是線性模型,。 二,、相關應用例子：Binary Logistic（因變量只能取兩個值1和0虛擬因變量） 案例一：本文用例來自于John Maindonald所著的《Data Analysis and Graphics Using R》一書，其中所用的數(shù)據(jù)集是anesthetic,，數(shù)據(jù)集來自于一組醫(yī)學數(shù)據(jù),，其中變量conc表示麻醉劑的用量，move則表示手術(shù)病人是否有所移動,，而我們用nomove做為因變量,，因為研究的重點在于conc的增加是否會使nomove的概率增加。 首先載入數(shù)據(jù)集并讀取部分文件,，為了觀察兩個變量之間關系,，我們可以利cdplot函數(shù)來繪制條件密度圖 install.packages("DAAG") library(lattice) library(DAAG) head(anesthetic) move conc logconc nomove 1 0 1.0 0.0000000 1 2 1 1.2 0.1823216 0 3 0 1.4 0.3364722 1 4 1 1.4 0.3364722 0 5 1 1.2 0.1823216 0 6 0 2.5 0.9162907 1 cdplot(factor(nomove)~conc,data=anesthetic,main='條件密度圖',ylab='病人移動',xlab='麻醉劑量') 從圖中可見，隨著麻醉劑量加大,，手術(shù)病人傾向于靜止,。下面利用logistic回歸進行建模，得到intercept和conc的系數(shù)為-6.47和5.57,，由此可見麻醉劑量超過1.16(6.47/5.57)時,，病人靜止概率超過50%。 anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic) summary(anes1) 結(jié)果顯示： Call: glm(formula = nomove ~ conc, family = binomial(link = "logit"), data = anesthetic) Deviance Residuals: Min 1Q Median 3Q Max -1.76666 -0.74407 0.03413 0.68666 2.06900 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -6.469 2.418 -2.675 0.00748 ** conc 5.567 2.044 2.724 0.00645 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 41.455 on 29 degrees of freedom Residual deviance: 27.754 on 28 degrees of freedom AIC: 31.754 Number of Fisher Scoring iterations: 5 下面做出模型的ROC曲線 anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic) 對模型做出預測結(jié)果 pre=predict(anes1,type='response') 將預測概率pre和實際結(jié)果放在一個數(shù)據(jù)框中 data=data.frame(prob=pre,obs=anesthetic$nomove) 將預測概率按照從低到高排序 data=data[order(data$prob),] n=nrow(data) tpr=fpr=rep(0,n) 根據(jù)不同的臨界值threshold來計算TPR和FPR,，之后繪制成圖 for (i in 1:n){ threshold=data$prob[i] tp=sum(data$prob>threshold&data$obs==1) fp=sum(data$prob>threshold&data$obs==0) tn=sum(data$prob fn=sum(data$prob tpr[i]=tp/(tp+fn) #真正率 fpr[i]=fp/(tn+fp) #假正率 } plot(fpr,tpr,type='l') abline(a=0,b=1) R中也有專門繪制ROC曲線的包,，如常見的ROCR包，它不僅可以用來畫圖,，還能計算ROC曲線下面面積AUC,以評價分類器的綜合性能,，該數(shù)值取0-1之間，越大越好,。 library(ROCR) pred=prediction(pre,anesthetic$nomove) performance(pred,'auc')@y.values perf=performance(pred,'tpr','fpr') plot(perf) 還可以使用更加強大的pROC包,，它可以方便的比較兩個分類器，并且能自動標出最優(yōu)臨界點,，圖形看起來比較漂亮： install.packages("pROC") library(pROC) modelroc=roc(anesthetic$nomove,pre) plot(modelroc,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="blue",print.thres=TRUE) 上面的方法是使用原始的0-1數(shù)據(jù)進行建模,即每一行數(shù)據(jù)均表示一個個體,，另一種是使用匯總數(shù)據(jù)進行建模，先將原始數(shù)據(jù)按下面步驟進行匯總 anestot=aggregate(anesthetic[,c('move','nomove')],by=list(conc=anesthetic$conc),FUN=sum) 結(jié)果如下： conc move nomove 1 0.8 6 1 2 1.0 4 1 3 1.2 2 4 4 1.4 2 4 5 1.6 0 4 6 2.5 0 2 anestot$conc=as.numeric(as.character(anestot$conc)) anestot$total=apply(anestot[,c('move','nomove')],1,sum) anestot$total [1] 7 5 6 6 4 2 anestot$prop=anestot$nomove/anestot$total anestot$prop [1] 0.1428571 0.2000000 0.6666667 0.6666667 1.0000000 1.0000000 對于匯總數(shù)據(jù),，有兩種方法可以得到同樣的結(jié)果,，一種是將兩種結(jié)果的向量合并做為因變量,，如anes2模型。另一種是將比率做為因變量,，總量做為權(quán)重進行建模,，如anes3模型。這兩種建模結(jié)果是一樣的,。 anes2=glm(cbind(nomove,move)~conc,family=binomial(link='logit'),data=anestot) summary(anes2) 結(jié)果顯示如下： Call: glm(formula = cbind(nomove, move) ~ conc, family = binomial(link = "logit"), data = anestot) Deviance Residuals: 1 2 3 4 5 6 0.20147 -0.45367 0.56890 -0.70000 0.81838 0.04826 Coefficients: Estimate Std. Error z value Pr(>z) (Intercept) -6.469 2.419 -2.675 0.00748 ** conc 5.567 2.044 2.724 0.00645 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 15.4334 on 5 degrees of freedom Residual deviance: 1.7321 on 4 degrees of freedom AIC: 13.811 Number of Fisher Scoring iterations: 5 anes3=glm(prop~conc,family=binomial(link='logit'),weights=total,data=anestot) 結(jié)果和上面的一樣,。 根據(jù)logistic模型，我們可以使用predict函數(shù)來預測結(jié)果,，下面根據(jù)上述模型來繪圖 x=seq(from=0,to=3,length.out=30) y=predict(anes1,data.frame(conc=x),type='response') plot(prop~conc,pch=16,col='red',data=anestot,xlim=c(0.5,3),main='Logistic回歸曲線圖',ylab='病人靜止概率',xlab='麻醉劑量') lines(y~x,lty=2,col='blue') 案例二：利用iris數(shù)據(jù)集,，進行邏輯回歸二分類測試，該數(shù)據(jù)集是R語言自帶得數(shù)據(jù)集,，包括四個屬性,，和三個分類。邏輯回歸我們用glm函數(shù)實現(xiàn),，該函數(shù)提供了各種類型的回歸,，如：提供正態(tài)、指數(shù),、gamma,、逆高斯、Poisson,、二項,。我們用的logistic回歸使用的是二項分布族binomial。 index <- which(iris$Species == 'setosa') 將種類為setosa的數(shù)據(jù)排除出我們需要的數(shù)據(jù)集 ir <- iris[- index,] levels(ir$Species)[1] <- '' 生成訓練集 split <- sample(100,100*(2/3)) ir_train <- ir[split,] 生成測試集 通過訓練集建立模型 summary(model) 模型運行結(jié)果： Call: glm(formula = Species ~ ., family = binomial(link = "logit"),data = ir_train) Deviance Residuals: Min 1Q Median 3Q Max -1.339e-04 -2.100e-08 2.100e-08 2.100e-08 1.059e-04 Coefficients: Estimate Std. Error z value Pr(>z) (Intercept) -1502.72 363247.01 -0.004 0.997 Sepal.Length 12.45 66482.13 0.000 1.000 Sepal.Width -285.61 95437.92 -0.003 0.998 Petal.Length 154.76 115968.97 0.001 0.999 Petal.Width 869.60 204513.80 0.004 0.997 (Dispersion parameter for binomial family taken to be 1) Null deviance: 9.0949e+01 on 65 degrees of freedom Residual deviance: 4.0575e-08 on 61 degrees of freedom AIC: 10 Number of Fisher Scoring iterations: 25 通過anova()函數(shù) 對模型進行方差分析 anova(model, test="Chisq") 方差分析如下： Analysis of Deviance Table Model: binomial, link: logit Response: Species Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 65 90.949 Sepal.Length 1 18.934 64 72.015 1.353e-05 *** Sepal.Width 1 0.131 63 71.884 0.7176 Petal.Length 1 51.960 62 19.924 5.665e-13 *** Petal.Width 1 19.924 61 0.000 8.058e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 下面通過McFadden R2指標進一步對模型進行分析 install.packages("pscl") library(pscl) pR2(model) llh llhNull G2 McFadden r2ML r2CU -2.028752e-08 -4.547461e+01 9.094922e+01 1.000000e+00 7.479224e-01 1.000000e+00 為了得到分類結(jié)果,，需要設定一個閾值p0——當p大于p0時,，認為該樣本的響應變量為1，否則為0,。閾值大小對模型的預測效果有較大影響,，需要進一步考慮。首先必須明確模型預測效果的評價指標,。 求解訓練模型的最佳閥值 對模型做出預測結(jié)果 model <- glm(Species ~.,family=binomial(link='logit'),data=ir_train) pre1=predict(model,type='response') 將預測概率pre1和實際結(jié)果放在一個數(shù)據(jù)框中 data1=data.frame(prob=pre1,obs=ifelse(ir_train$Species=="virginica",1,0)) 將預測概率按照從低到高排序 data1=data1[order(data1$prob),] n=nrow(data1) tpr=fpr=rep(0,n) 根據(jù)不同的臨界值threshold來計算TPR和FPR,，之后繪制成圖 for (i in 1:n){ threshold=data1$prob[i] tp=sum(data1$prob>threshold&data1$obs==1) fp=sum(data1$prob>threshold&data1$obs==0) tn=sum(data$prob fn=sum(data$prob tpr[i]=tp/(tp+fn) #真正率 fpr[i]=fp/(tn+fp) #假正率 } plot(fpr,tpr,type='l') abline(a=0,b=1)
下面通過pROC包自動標出最優(yōu)臨界點(0.506) install.packages("pROC") library(pROC) modelroc1=roc(ifelse(ir_train$Species=="virginica",1,0),pre1) plot(modelroc1,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c("green","red"),max.auc.polygon=TRUE,auc.polygon.col="skyblue",print.thres=TRUE)
評估模型的預測效果 predict <- predict(model,type='response',newdata=ir_test) predict.results <- ifelse( predict> 0.506,"virginica","versicolor") misClasificError <- mean(predict.results != ir_test$Species) print(paste('Accuracy',1-misClasificError)) [1] "Accuracy 1" 最后一步，我們將通過畫ROC曲線,，并計算其AUC面積,，作為評估二類分類效果的一個典型測量 install.packages("ROCR") library(gplots) library(ROCR) p <- predict(model,type='response',newdata=ir_test) p.results <- ifelse( p> 0.5,1,0) pr <- prediction(p.results, ifelse(ir_test$Species=="virginica",1,0)) prf <- performance(pr, measure = "tpr", x.measure = "fpr") plot(prf) auc <- performance(pr, measure = "auc") auc <- [email protected][[1]] 0.9285714
auc real <- ir_test$Species data.frame(real,predict) res <- data.frame(real,predict =ifelse(predict>0.5,'virginca','versicorlor')) 查看模型效果 plot(res)
基于R的案例 以下這個例子主要參考《Introduction to statistical learning with R》，強烈推薦這本書,。尤其是回歸部分，講解的很透徹,。 數(shù)據(jù)以Smarket為例,，該數(shù)據(jù)包含2001-2005年標準普爾500指數(shù)1250個觀測值,，9個變量。9個變量分別為：年份,，前5個交易日的收益率,，前一個交易日和交易額（單位：billions），今天的回報率和走勢（上升或下降）,。 df=read.csv("house.csv",header=TRUE) str(df) data.frame': 1250 obs. of 9 variables: $ Year : num 2001 2001 2001 2001 2001 ... $ Lag1 : num 0.381 0.959 1.032 -0.623 0.614 ... $ Lag2 : num -0.192 0.381 0.959 1.032 -0.623 ... $ Lag3 : num -2.624 -0.192 0.381 0.959 1.032 ... $ Lag4 : num -1.055 -2.624 -0.192 0.381 0.959 ... $ Lag5 : num 5.01 -1.055 -2.624 -0.192 0.381 ... $ Volume : num 1.19 1.3 1.41 1.28 1.21 ... $ Today : num 0.959 1.032 -0.623 0.614 0.213 ... $ Direction: Factor w/ 2 levels "Down","Up": 2 2 1 2 2 2 1 2 2 2 ... #以前5個交易日的收益率和前一個工作日的交易額為自變量,，今天的走勢做因變量做Logistic回歸； glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume, data=Smarket,family="binomial") summary(glm.fit) Call: glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, family = "binomial", data = Smarket) Deviance Residuals: Min 1Q Median 3Q Max -1.446 -1.203 1.065 1.145 1.326 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.126000 0.240736 -0.523 0.601 Lag1 -0.073074 0.050167 -1.457 0.145 Lag2 -0.042301 0.050086 -0.845 0.398 Lag3 0.011085 0.049939 0.222 0.824 Lag4 0.009359 0.049974 0.187 0.851 Lag5 0.010313 0.049511 0.208 0.835 Volume 0.135441 0.158360 0.855 0.392 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1731.2 on 1249 degrees of freedom Residual deviance: 1727.6 on 1243 degrees of freedom AIC: 1741.6 Number of Fisher Scoring iterations: 3 有時候,，預測股市就是這么不靠譜,。很多自變量沒通過驗證，接下來我們基于AIC準則用逐步回歸做一下變量篩選,； 注：AIC一般公式為 AIC=2k-ln(L),其中k為參數(shù)的個數(shù),，L為估計模型最大化的似然函數(shù)值； step(glm.fit) Start: AIC=1741.58 Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume Df Deviance AIC - Lag4 1 1727.6 1739.6 - Lag5 1 1727.6 1739.6 - Lag3 1 1727.6 1739.6 - Lag2 1 1728.3 1740.3 - Volume 1 1728.3 1740.3 <none> 1727.6 1741.6 - Lag1 1 1729.7 1741.7 #此處略去一百字,； Direction ~ 1 Call: glm(formula = Direction ~ 1, family = "binomial", data = Smarket) Coefficients: (Intercept) 0.07363 Degrees of Freedom: 1249 Total (i.e. Null); 1249 Residual Null Deviance: 1731 Residual Deviance: 1731 AIC: 1733 這個結(jié)果足以讓你崩潰,，扔硬幣都比這靠譜。原來不用任何變量預測是最靠譜的,。 #模型的評估 glm.probs =predict(glm.fit,type ="response") glm.probs[1:10] #生成啞變量 contrasts(Smarket$Direction) glm.pred=rep ("Down " ,1250) glm.pred[glm .probs >.5]=" Up" table(glm .pred ,Direction ) mean(glm.pred== Direction ) [1] 0.5216 通過混淆矩陣計算分類的準確性僅有52%,，很不樂觀； 挖掘本質(zhì)上是挖掘有趣的的模式,，以上數(shù)據(jù)說明我們白忙活了一場,，但是我們通過實例學習了Logistic模型。 當然我們還可以通過數(shù)據(jù)變換或構(gòu)造新的變量來提升擬合降低AIC,，最終,，模型講究的還是泛化能力； 有時候：擬合很豐滿,，泛化很骨感,。 因變量多分類 logistic 回歸 1、概述：多元Logistic回歸模型被用來建立有多個輸出變量的模型,，且這些預測變量通過一個線性組合變成為一個最終的預測變量,。Multinomial Logistic 回歸模型中因變量可以取多個值.例如一個典型的例子就是分類一部電影，它是“有趣”,、“一般般”還是“無聊”,。 所需的應用包： library(foreign) library(nnet) library(ggplot2) library(reshape2) 數(shù)據(jù)集我們將以Titanic數(shù)據(jù)集作為實例進行分析。關于這個數(shù)據(jù)集,，在網(wǎng)上流傳了好幾個免費的版本,，而我則建議大家使用在Kaggle上的那個版本，因為它已經(jīng)被用過很多次了（如果你要下載這個數(shù)據(jù)集,，你需要在Kaggle注冊一個賬號）,。 數(shù)據(jù)清洗過程當我們拿到一個真實的數(shù)據(jù)集的時候,，我們需要確認這個數(shù)據(jù)集有哪些缺失值或者異常值。因此,，我們需要做這方面的準備工作以便我們后面對此進行數(shù)據(jù)分析,。第一步則是我們使用read.csv()來讀取數(shù)據(jù)集，可以在這里下載數(shù)據(jù)集,。
現(xiàn)在，我們要查閱一下缺失值,，并觀察一下每個變量存在多少異常值,，使用sapply()函數(shù)來展示數(shù)據(jù)框的每一列的值。
在這里對缺失值進行可視化會有幫助：Amelia包有一個特別的作圖函數(shù)missmap()來針對我們的數(shù)據(jù)集進行作圖,，并且對缺失值做標記,。
這個變量有非常多的缺失值，我們都不會使用他們,。我們也會把Passengerld去掉,，因為它只有索引和票。
觀察一下缺失值現(xiàn)在，我們需要觀察一下別的缺失值,。R可以通過在擬合函數(shù)內(nèi)設定一個參數(shù)集來產(chǎn)生擬合一個廣義線性模型,。不過，我個人認為,，我們應當盡可能的手動去除這些缺失值,。這里有很多種方法來說實現(xiàn)，而一種常用的方法就是以均值，或中位數(shù),，或者是一個眾數(shù)來代替缺失值,。這里,，我們會使用均值,。
為了能更容易的理解R是怎么處理分類變量的,，我們可以使用contrasts()函數(shù),。這個函數(shù)向我們展示了R是如何優(yōu)化變量的，并把它轉(zhuǎn)譯成模型,。
比如,，你可以看見在變量sex，female可被作為參考,。在乘船上的缺失值,，由于這里只有兩行，我們可以把這兩行刪掉(我們也可以用眾數(shù)取代缺失值,，然后保留這些數(shù)據(jù)點）,。
在繼續(xù)擬合過程之前，讓我們回顧一下數(shù)據(jù)清洗和重構(gòu)的重要性,。這樣的預處理通常對于獲得一個好而且預測能力強的模型是非常重要的,。 模型擬合我們把數(shù)據(jù)分成兩個部分，訓練和測試集,。訓練數(shù)據(jù)集可用于模型的擬合,，而測試數(shù)據(jù)集則對此進行測試。
現(xiàn)在,，讓我們擬合一下這個模型,。確認我們在glm()函數(shù)使用參數(shù)family=binomial。
現(xiàn)在,，使用summary()函數(shù)來獲得這個模型的結(jié)果,。
解讀邏輯回歸模型的結(jié)果現(xiàn)在，我們可以分析這個擬合的模型,，閱讀一下它的結(jié)果,。首先，我們可以看到SibSp,，F(xiàn)are和Embarked并不是重要的變量,。對于那些比較重要的變量，sex的p值最低,，這說明sex與乘客的生存幾率的關系是最強的,。這個負系數(shù)預測變量表明其它變量都一樣,，男乘客生存的機率更低。記住,，邏輯模型的因變量是對數(shù)機率：ln(odds) = ln(p/(1-p)) = ax1 + bx2 + … + z*xn,。male是一個優(yōu)化變量，男性的生還機率下載2.75個對數(shù)機率,，而年齡的增大則下降0.037個對數(shù)機率,。
空偏差和殘差的不同展示了我們的模型是空模型（只有截距）,。這個溝閱讀,，效果越好。分析這個圖表,，我們可以看到每增加一個變量,，同時也會降低其中的偏差。再一次,，添加Pclass,、Sex和Age能顯著的降低殘差。其它變量對于模型的提高沒什么用,，盡管SibSp的p值很低,。一個p值較大的值預示著一個沒有變量的模型多少解釋了模型的變化。最終,，你想要看到的就是偏差的降低和AIC,。
評估模型的預測能力在上述的步驟中,，我們簡單的評估一下模型的擬合效果。現(xiàn)在,，我們看一下使用心得數(shù)據(jù)集,，它的預測結(jié)果y是什么樣的。設定參數(shù)type=’response’,，R可以輸出形如P(Y=1|X）的概率,。我們的預測邊界就是0.5。如果P(Y=1|X)>0.5,，y=1或y=0,。記住，其它的決策邊界可能是更好的選擇,。
精度為0.84,，說明這個模型擬合效果不錯。然而，你要記住,，結(jié)果多少都依賴于手動改變數(shù)據(jù)集的數(shù)據(jù),。因此，如果你想要得到更精確的精度,，你最好執(zhí)行一些交叉檢驗,，例如k次交叉檢驗。
下面是ROC圖： |
|
|
