R實戰案例：利用演算法識別糖尿病患者（R語言實現）

01-28

Part One：作業要求

數據源：PimaIndiansDiabetes2.csv

數據介紹：印第安人糖尿病資料庫，包括以下信息：

作業需求：我們需要利用"pregnant","glucose", "pressure", "triceps", "insulin","mass",

"pedigree","age"等變數來預測"diabetes"的值。因為數據有缺失值，需要對數據先進行探索，完成缺失值的插補後，再利用分類演算法建立預測模型，識別糖尿病患者。

項目要求：

1、 數據探索：利用課上講到的列表顯示缺失值和圖形探究缺失數據兩種方式對缺失值模式進行探究。

2、 數據完善：要求不能直接刪除缺失數據，至少需要利用兩種方式對缺失值進行插補。

3、 數據分區：需要按照變數diabetes來進行等比例抽樣，其中80%作為訓練集train數據，20%作為測試集test數據。

4、 建立模型及評估：利用分類演算法對train數據集建立分類預測模型，並對test數據集進行預測，利用混淆矩陣查看模型評估效果。

Part Two：個人探索過程

1. 數據探索：

a. 列表顯示缺失值

library(mice)

md.pattern(indian)

b. 圖形探究缺失數據(aggr()函數)

library(VIM)

aggr(indian,prop=T,numbers=T)

aggr(indian,prop=F,numbers=T)

aggr(indian,prop=F,numbers=F)

小結：根據列表和圖形結果，共有五列數據存在缺失值："glucose", "mass","pressure", "triceps", "insulin",分別缺失5,11,35,227,374個值。

2. 數據完善：

a. 隨機森林插補

library(missForest) # 通過隨機森林進行插補

z=missForest(indian)

indian.mf=z$ximp

md.pattern(indian.mf) # 列表驗證

aggr(indian.mf,prop=T,numbers=T) # 圖形驗證

b. 回歸模型插補

# 回歸模型插補

indian.lm <- indian

# 插補glucose列缺失值

ind<-which(is.na(indian.lm[,2])==T) # 返回地2列為NA的行號

data_NL <- indian.lm[-ind,] #獲取第2列所有不為NA的數據

data_NA <- indian.lm[ind,] #獲取第2列為NA的數據

# 構建回歸模型

lm = lm(glucose~preg.nt+pedigree+age,data=data_NL)

summary(lm) #查看模型的顯著性：pedigree+age極其顯著，preg.nt不顯著

lm = lm(glucose~pedigree+age,data=data_NL) # 重新構建回歸模型

indian.lm[ind,2] =round(predict(lm,data_NA)) #插補數據

依次通過以上方法對其它欄位進行插補：

# 插補pressure列缺失值

ind<-which(is.na(indian.lm[,3])==T) # 返回地3列為NA的行號

data_NL <- indian.lm[-ind,] #獲取第3列所有不為NA的數據

data_NA <- indian.lm[ind,] #獲取第3列為NA的數據

# 構建回歸模型

lm = lm(pressure~preg.nt+pedigree+age,data=data_NL)

summary(lm) #查看模型的顯著性：age極其顯著

lm = lm(pressure~age,data=data_NL) # 重新構建回歸模型

indian.lm[ind,3] =round(predict(lm,data_NA))

# 插補triceps列缺失值

ind<-which(is.na(indian.lm[,4])==T) # 返回地4列為NA的行號

data_NL <- indian.lm[-ind,] #獲取第4列所有不為NA的數據

data_NA <- indian.lm[ind,] #獲取第4列為NA的數據

# 構建回歸模型

lm = lm(triceps~preg.nt+pedigree+age,data=data_NL)

summary(lm) #查看模型的顯著性：pedigree+age顯著

lm = lm(triceps~pedigree+age,data=data_NL) # 重新構建回歸模型

indian.lm[ind,4] =round(predict(lm,data_NA))

# 插補insulin列缺失值

ind<-which(is.na(indian.lm[,5])==T) # 返回地5列為NA的行號

data_NL <- indian.lm[-ind,] #獲取第5列所有不為NA的數據

data_NA <- indian.lm[ind,] #獲取第5列為NA的數據

# 構建回歸模型

lm = lm(insulin~preg.nt+pedigree+age,data=data_NL)

summary(lm) #查看模型的顯著性：preg.nt+pedigree+age顯著

indian.lm[ind,5] =round(predict(lm,data_NA))

# 插補mass列缺失值

ind<-which(is.na(indian.lm[,6])==T) # 返回地6列為NA的行號

data_NL <- indian.lm[-ind,] #獲取第6列所有不為NA的數據

data_NA <- indian.lm[ind,] #獲取第6列為NA的數據

# 構建回歸模型

lm = lm(mass~preg.nt+pedigree+age,data=data_NL)

summary(lm) #查看模型的顯著性：pedigree顯著

lm = lm(mass~pedigree,data=data_NL) # 重新構建回歸模型

indian.lm[ind,6] =round(predict(lm,data_NA))

md.pattern(indian.lm) #列表驗證回歸插補結果

aggr(indian.lm,prop=T,numbers=T) # 圖形驗證

小結：通過兩種方式（隨機森林，回歸建模）對數據集進行插補得到兩組數據集indian.mf、indian.lm。

3. 數據分區與建模

### 數據分區

library(caret)

# 構建result 存放結果

result <- data.frame(model=c("naiveBayes","C5.0","CART","ctree","bagging","boosting","隨機森林"),errTrain_lm=rep(0,7),errTest_lm=rep(0,7),errTrain_rf=rep(0,7),errTest_rf=rep(0,7))

# 構建result.10 存放結果(10折交叉)

result.10 <-data.frame(model=c("決策樹","隨機森林","人工神經網路"),errTrain_lm=rep(0,3),errTest_lm=rep(0,3),errTrain_rf=rep(0,3),errTest_rf=rep(0,3))

##構建10折交叉驗證

library(rpart)

control <- trainControl(method = "repeatedcv",number = 10,repeats = 3)

rpart.model.mf<- train(diabetes~.,data=indian.mf,method="rpart",trControl=control)

rpart.model.lm<- train(diabetes~.,data=indian.lm,method="rpart",trControl=control)

# 隨機森林

rf.model.mf <- train(diabetes~.,data=indian.mf,method="rf",trControl=control) #8

rf.model.lm <- train(diabetes~.,data=indian.lm,method="rf",trControl=control) #2

# 神經網路

nnet.model.mf <- train(diabetes~.,data=indian.mf,method="nnet",trControl=control)

nnet.model.lm <- train(diabetes~.,data=indian.lm,method="nnet",trControl=control)

# 查看結果

rpart.model.mf

rpart.model.lm

rf.model.mf

rf.model.lm

nnet.model.mf

nnet.model.lm

通過循環，把兩種組數據（隨機森林補差，回歸補差）進行普通建模，10折交叉驗證獲取最優參數建模，構建混淆矩陣，評估模型。

for(a in 1:2)

{

# 分別使用回歸模型補差數據集隨機森林補差數據集進行建模

ind <- createDataPartition(switch(a,indian.lm$diabetes,indian.mf$diabetes),times = 1,p=0.8,list = F)

train <- indian.lm[ind,] # 構建訓練集

test <- indian.lm[-ind,] # 構建測試集

prop.table(table(indian.lm$diabetes))

prop.table(table(train$diabetes))

prop.table(table(test$diabetes))

### 建模和評估

# 使用naiveBayes函數建立樸素貝葉斯分類器

library(e1071)

str(indian.lm)

naiveBayes.model <- naiveBayes(diabetes~.,data = train) # 構建模型

# 預測結果

train_predict <- predict(naiveBayes.model,newdata=train)

test_predict <- predict(naiveBayes.model,newdata=test)

# 構建混淆矩陣

tableTrain <- table(actual=train$diabetes,predict=train_predict)

tableTest <- table(actual=test$diabetes,predict=test_predict)

# 計算誤差率

errTrain <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")

errTest <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")

# 決策樹模型

library(C50)

C5.0.model <-C5.0(diabetes~.,data=train)

library(rpart)

rpart.model <-rpart(diabetes~.,data=train)

library(party)

ctree.model <- ctree(diabetes~.,data=train)

result[1,switch(a,2,4)] <-errTrain

result[1,switch(a,3,5)] <-errTest

for(i in 1:3){

# 預測結果

train_predict <-predict(switch(i,C5.0.model,rpart.model,ctree.model),newdata=train,type=switch(i,"class","class","response"))

test_predict <-predict(switch(i,C5.0.model,rpart.model,ctree.model),newdata=test,type=switch(i,"class","class","response"))

# 構建混淆矩陣

tableTrain <- table(actual=train$diabetes,predict=train_predict)

tableTest <- table(actual=test$diabetes,predict=test_predict)

# 計算誤差率

result[i+1,switch(a,2,4)] <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")

result[i+1,switch(a,3,5)] <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")

}

result

# 使用隨機森林建模

library(adabag)

bagging.model <- bagging(diabetes~.,data=train)

boosting.model <-boosting(diabetes~.,data=train)

library(randomForest)

randomForest.model <- randomForest(diabetes~.,data=train)

for(i in 1:3){

# 預測結果

train_predict <-predict(switch(i,bagging.model,boosting.model,randomForest.model),newdata = train)

test_predict <- predict(switch(i,bagging.model,boosting.model,randomForest.model),newdata = test)

# 構建混淆矩陣

tableTrain <- table(actual=train$diabetes,predict=switch(i,train_predict$class,train_predict$class,train_predict))

tableTest <- table(actual=test$diabetes,predict=switch(i,test_predict$class,test_predict$class,test_predict))

# 計算誤差率

result[i+4,switch(a,2,4)] <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")

result[i+4,switch(a,3,5)] <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")

}

# 通過10折交叉驗證選擇最優參數

library(caret)

# 查看結果

rpart.model1

rf.model

nnet.model

# 利用rpart函數構建分類樹

rpart.model1 <- rpart::rpart(diabetes~.,data = train,control = c(cp=switch(a,0.01741294, 0.02985075)))

# 利用randomForest函數構建隨機森林

rf.model <-randomForest::randomForest(diabetes~.,data = train,mtry=switch(a,2,8))

# 利用nnet函數建立人工神經網路

nnet.model <- nnet::nnet(diabetes~.,data = train,size=switch(a,5,3),decay=0.1)

for(i in 1:3){

# 預測結果

train_predict <- predict(switch(i,rpart.model1,rf.model,nnet.model),newdata = train,type="class")

test_predict <- predict(switch(i,rpart.model1,rf.model,nnet.model),newdata = test,type="class")

# 構建混淆矩陣

tableTrain <- table(actual=train$diabetes,predict=train_predict)

tableTest <- table(actual=test$diabetes,predict=test_predict)

# 計算誤差率

result.10[i,switch(a,2,4)] <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")

result.10[i,switch(a,3,5)] <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")

}

# 查看模型評估

result

# 10折交叉驗證結果

result.10

總結：通過對兩組數據建模，構建混淆矩陣的結果主要可以得出如下結論：

1. 整體而言，隨機森林的模型要優於其它模型

2. 隨機森林中使用randomForest和boosting函數建模的所有訓練集誤差率為0%，測試集數據的誤差率達到24~30%

3. 整體而言，當前的所有模型還沒有達到最優，需要進一步探索。

PimaIndiansDiabetes2.zip

作者：王小二不在家 R語言中文社區專欄作者
出處：R語言作業：利用演算法識別糖尿病患者
公眾號：R語言中文社區（r_shequ）加微信：tstianwang，可以到R語言中文社區群，跟各個作者交流互動。