R語言邏輯迴歸、ROC曲線和十折交叉驗證

Tiaaaaa發表於2017-02-27
R語言邏輯迴歸

自己整理編寫的邏輯迴歸模板，作為學習筆記記錄分享。資料集用的是14個自變數Xi，一個因變數Y的australian資料集。


1. 測試集和訓練集3、7分組

australian <- read.csv("australian.csv",as.is = T,sep=",",header=TRUE)
#讀取行數
N = length(australian$Y)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
#ind=1的是0.7概率出現的行，ind=2是0.3概率出現的行
ind=sample(2,N,replace=TRUE,prob=c(0.7,0.3))
#生成訓練集(這裡訓練集和測試集隨機設定為原資料集的70%，30%)
aus_train <- australian[ind==1,]
#生成測試集
aus_test <- australian[ind==2,]

2.生成模型，結果匯出 

#生成logis模型，用glm函式
#用訓練集資料生成logis模型，用glm函式
#family：每一種響應分佈（指數分佈族）允許各種關聯函式將均值和線性預測器關聯起來。常用的family：binomal(link='logit')--響應變數服從二項分佈，連線函式為logit，即logistic迴歸
pre <- glm(Y ~.,family=binomial(link = "logit"),data = aus_train)
summary(pre)

#測試集的真實值
real <- aus_test$Y
#predict函式可以獲得模型的預測值。這裡預測所需的模型物件為pre，預測物件newdata為測試集,預測所需型別type選擇response,對響應變數的區間進行調整
predict. <- predict.glm(pre,type='response',newdata=aus_test)
#按照預測值為1的概率，>0.5的返回1，其餘返回0
predict =ifelse(predict.>0.5,1,0)
#資料中加入預測值一列
aus_test$predict = predict
#匯出結果為csv格式
#write.csv(aus_test,"aus_test.csv")

3.模型檢驗 

##模型檢驗
res <- data.frame(real,predict)
#訓練資料的行數，也就是樣本數量
n = nrow(aus_train)      
#計算Cox-Snell擬合優度
R2 <- 1-exp((pre$deviance-pre$null.deviance)/n)    
cat("Cox-Snell R2=",R2,"\n")
#計算Nagelkerke擬合優度，我們在最後輸出這個擬合優度值
R2<-R2/(1-exp((-pre$null.deviance)/n))  
cat("Nagelkerke R2=",R2,"\n")
##模型的其他指標
#residuals(pre)     #殘差
#coefficients(pre)  #係數，線性模型的截距項和每個自變數的斜率，由此得出線性方程表示式。或者寫為coef(pre)
#anova(pre)         #方差

4.準確率和精度 

true_value=aus_test[,15]
predict_value=aus_test[,16]
#計算模型精確度
error = predict_value-true_value
accuracy = (nrow(aus_test)-sum(abs(error)))/nrow(aus_test) #精確度--判斷正確的數量佔總數的比例
#計算Precision，Recall和F-measure
#一般來說，Precision就是檢索出來的條目（比如：文件、網頁等）有多少是準確的，Recall就是所有準確的條目有多少被檢索出來了
#和混淆矩陣結合，Precision計算的是所有被檢索到的item（TP+FP）中,"應該被檢索到的item（TP）”佔的比例；Recall計算的是所有檢索到的item（TP）佔所有"應該被檢索到的item（TP+FN）"的比例。
precision=sum(true_value & predict_value)/sum(predict_value)  #真實值預測值全為1 / 預測值全為1 --- 提取出的正確資訊條數/提取出的資訊條數
recall=sum(predict_value & true_value)/sum(true_value)  #真實值預測值全為1 / 真實值全為1 --- 提取出的正確資訊條數 /樣本中的資訊條數
#P和R指標有時候會出現的矛盾的情況，這樣就需要綜合考慮他們，最常見的方法就是F-Measure（又稱為F-Score）
F_measure=2*precision*recall/(precision+recall)    #F-Measure是Precision和Recall加權調和平均，是一個綜合評價指標
#輸出以上各結果
print(accuracy)
print(precision)
print(recall)
print(F_measure)
#混淆矩陣，顯示結果依次為TP、FN、FP、TN
table(true_value,predict_value)

5.ROC曲線的幾個方法 

#ROC曲線
# 方法1
#install.packages("ROCR")  
library(ROCR)     
pred <- prediction(predict.,true_value)   #預測值(0.5二分類之前的預測值)和真實值   
performance(pred,'auc')@y.values        #AUC值
perf <- performance(pred,'tpr','fpr')
plot(perf)
#方法2
#install.packages("pROC")
library(pROC)
modelroc <- roc(true_value,predict.)
plot(modelroc, print.auc=TRUE, auc.polygon=TRUE,legacy.axes=TRUE, grid=c(0.1, 0.2),
     grid.col=c("green", "red"), max.auc.polygon=TRUE,
     auc.polygon.col="skyblue", print.thres=TRUE)        #畫出ROC曲線，標出座標，並標出AUC的值
#方法3，按ROC定義
TPR=rep(0,1000)
FPR=rep(0,1000)
p=predict.
for(i in 1:1000)
  { 
  p0=i/1000;
  ypred<-1*(p>p0)  
  TPR[i]=sum(ypred*true_value)/sum(true_value)  
  FPR[i]=sum(ypred*(1-true_value))/sum(1-true_value)
  }
plot(FPR,TPR,type="l",col=2)
points(c(0,1),c(0,1),type="l",lty=2)

6.更換測試集和訓練集的選取方式，採用十折交叉驗證 

australian <- read.csv("australian.csv",as.is = T,sep=",",header=TRUE)
#將australian資料分成隨機十等分
#install.packages("caret")
#固定folds函式的分組
set.seed(7)
require(caret)
folds <- createFolds(y=australian$Y,k=10)

#構建for迴圈，得10次交叉驗證的測試集精確度、訓練集精確度

max=0
num=0

for(i in 1:10){
  
  fold_test <- australian[folds[[i]],]   #取folds[[i]]作為測試集
  fold_train <- australian[-folds[[i]],]   # 剩下的資料作為訓練集
  
  print("***組號***")
  
  fold_pre <- glm(Y ~.,family=binomial(link='logit'),data=fold_train)
  fold_predict <- predict(fold_pre,type='response',newdata=fold_test)
  fold_predict =ifelse(fold_predict>0.5,1,0)
  fold_test$predict = fold_predict
  fold_error = fold_test[,16]-fold_test[,15]
  fold_accuracy = (nrow(fold_test)-sum(abs(fold_error)))/nrow(fold_test) 
  print(i)
  print("***測試集精確度***")
  print(fold_accuracy)
  print("***訓練集精確度***")
  fold_predict2 <- predict(fold_pre,type='response',newdata=fold_train)
  fold_predict2 =ifelse(fold_predict2>0.5,1,0)
  fold_train$predict = fold_predict2
  fold_error2 = fold_train[,16]-fold_train[,15]
  fold_accuracy2 = (nrow(fold_train)-sum(abs(fold_error2)))/nrow(fold_train) 
  print(fold_accuracy2)
  
  
  if(fold_accuracy>max)
    {
    max=fold_accuracy  
    num=i
    }
  
}

print(max)
print(num)

##結果可以看到，精確度accuracy最大的一次為max,取folds[[num]]作為測試集，其餘作為訓練集。

7.得到十折交叉驗證的精確度，結果匯出 

#十折裡測試集最大精確度的結果
testi <- australian[folds[[num]],]
traini <- australian[-folds[[num]],]   # 剩下的folds作為訓練集
prei <- glm(Y ~.,family=binomial(link='logit'),data=traini)
predicti <- predict.glm(prei,type='response',newdata=testi)
predicti =ifelse(predicti>0.5,1,0)
testi$predict = predicti
#write.csv(testi,"ausfold_test.csv")
errori = testi[,16]-testi[,15]
accuracyi = (nrow(testi)-sum(abs(errori)))/nrow(testi) 

#十折裡訓練集的精確度
predicti2 <- predict.glm(prei,type='response',newdata=traini)
predicti2 =ifelse(predicti2>0.5,1,0)
traini$predict = predicti2
errori2 = traini[,16]-traini[,15]
accuracyi2 = (nrow(traini)-sum(abs(errori2)))/nrow(traini) 

#測試集精確度、取第i組、訓練集精確
accuracyi;num;accuracyi2
#write.csv(traini,"ausfold_train.csv")





相關文章