分析并预测糖尿病患者 R

data：diabetes.csv

---

一、缺失值处理

# 读入data
library(readr)
data <- read_csv("diabetes.csv")# 基本分析data
summary(data)# 判断NA & 异常值
library(VIM)
aggr(data,prop=F,numbers=T) # 直方图表示缺失值
data[!complete.cases(data),] # 缺失行# 异常值改为NA
data <- data %>% mutate(across(c("Glucose","BloodPressure","SkinThickness","Insulin","BMI"),~ifelse(. == 0,NA,.)))# 随机森林进行缺失值补充 （实际为rf中的mf）
# 适用：1.混合类型变量 & 非线性关系2.NMAR&高缺失比例（传统方法误差大）
# 不适用：1.大数据集（如>10万样本） 改用快速插补 2.MCAR 简单插补library(missForest)
set.seed(42)  # 确保可重复性
z <- missForest(diabetes, verbose = TRUE) 
diabetes.mf <- z$ximp  # .mf：后缀表示该数据是经 missForest() 函数处理后的结果，z 是一个列表，其中 ximp 为插补后的完整数据集 。# 评测随机森林填充效果
library(ggplot2)
ggplot() + geom_density(aes(x = data$Glucose, fill = "原始数据"), alpha = 0.5) + geom_density(aes(x = diabetes.mf$Glucose, fill = "填补后"), alpha = 0.5) +    labs(title = "Glucose分布对比") # 多重插补法 相当于自动多次回归
library(mice)
imputed_data <- mice(diabetes_mi,m = 5,           # 生成5套插补数据集maxit = 10,      # 最大迭代次数method = 'pmm',  # 使用预测均值匹配（适用于数值变量）seed = 42        # 确保可重复性
)complete_data <- complete(imputed_data, "long")  # 合并所有插补集
final_data <- complete(imputed_data, 1)  # 使用第一套插补结果
library(ggplot2)
# 效果检验：
# 1.分布一致性：对比插补前后Glucose的分布，确保右偏态等特征保留
ggplot() +geom_density(aes(x=diabetes_mi$Glucose, color="Original")) +geom_density(aes(x=final_data$Glucose, color="MICE"))# 2.对比前后数据分布
summary(original_data$Glucose)  # 原始数据
summary(imputed_data$Glucose)   # 填补后数据
# 线性回归进行缺失值补充
# 注意如果目标变量依靠的变量中存在NA ，则结果依然为NA 
diabetes.lm <- diabetes  # 复制原始数据集
ind <- which(is.na(diabetes.lm[,2]) == T)  # 找出第2列（Glucose）为NA的行号
data_NL <- diabetes.lm[-ind,]  # 非缺失数据（完整数据），用于训练回归模型
data_NA <- diabetes.lm[ind,]   # 缺失数据（需插补部分），用于后续预测插补lm <- lm(Glucose ~ Pregnancies + DiabetesPedigreeFunction + Age, data = data_NL)
summary(lm) # 评估 p>0.05，需剔除以简化模型 。# 优化 Pregnancies
lm <- lm(Glucose ~ DiabetesPedigreeFunction + Age, data = data_NL)
diabetes.lm[ind,2] <- round(predict(lm, data_NA))  # 预测并四舍五入summary(lm) # 通过R^2值进行判断拟合是否合适

二、基本数据分析

# 基本统计分析
summary(data)# 详细统计
library(psych)
describe(diabetes[c(列名)])# 可视化
#直方图
par(mfrow = c(1, 2)) 
hist(表$列名, main = "描述题目", xlab = "单位")
# 箱线图
boxplot(``)# 输出离群点
bp <- diabetes$BloodPressure[!is.na(diabetes$BloodPressure)]  # 排除NA值q1 <- quantile(bp, 0.25)q3 <- quantile(bp, 0.75)iqr <- q3 - q1lower_bound <- q1 - 1.5 * iqrupper_bound <- q3 + 1.5 * iqroutliers <- bp[bp < lower_bound | bp > upper_bound]print(outliers)  # 输出所有离群点数值# 偏度 & 峰度
# 偏态：若结果>0，数据右偏；若<0则左偏。
# 峰态：若结果>3，尖峰分布；<3低峰态（数据较平坦）。
library(moments)
skewness_val <- skewness(bloodpressure)  # 偏态系数
kurtosis_val <- kurtosis(bloodpressure)  # 峰态系数#正态性检验 # p-value<0.05非正态
shapiro.test(bloodpressure)
# 翘起 非正太
qqnorm(bloodpressure, main = "BloodPressure Q-Q Plot")
qqline(bloodpressure, col = "red")# 缺失值处理
# 检验某列是否有NA值
diabetes[!complete.cases(diabetes$BloodPressure),]
# 1.偏度不大 avg
# 2.存在异常值 median
median_bp <- median(diabetes$BloodPressure[diabetes$BloodPressure != 0], na.rm = TRUE)
diabetes$BloodPressure[diabetes$BloodPressure == 0] <- median_bp
# 3.是否存在线性关系
# 散点图观察  通过符号分类
# 1. 加载数据并计算相关系数矩阵
diabetes <- read.csv("diabetes.csv")
cor_matrix <- cor(diabetes[, -9], use = "complete.obs")  # 排除分类变量Outcome# 2. 定义符号分类函数
symbolize_cor <- function(r) {ifelse(abs(r) >= 0.7, "★★★",ifelse(abs(r) >= 0.3, "★★",ifelse(abs(r) >= 0.1, "★", "○")))
}# 3. 生成符号化矩阵
symbol_matrix <- apply(cor_matrix, c(1,2), symbolize_cor)
diag(symbol_matrix) <- "-"  # 对角线设为"-"# 4. 打印符号化矩阵（控制台输出）
print(symbol_matrix)# 5. 增强可视化（热图+符号叠加）
library(ggplot2)
library(reshape2)
melted_cor <- melt(cor_matrix) # 融化为长格式，适合可视化（如热力图）或分组分析# 数据映射与基础图层ggplot
#  图层叠加 geom_tile
# 符号化标记 geom_text
# 颜色梯度与图例   scale_fill_gradient2
# 主题与细节调整 theme_minimal()：简洁主题，移除背景网格和边框。 coord_fixed()：固定纵横比为 1，确保色块为正方形。
ggplot(melted_cor, aes(Var1, Var2, fill = value)) + geom_tile(color = "white") +geom_text(aes(label = symbolize_cor(value)), size = 4) +   scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0, limit = c(-1,1), name = "Correlation") +theme_minimal() +theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +coord_fixed() 

三、建模

# 数据分区&建模
# 1. 重构结果存储框架（移除填充方法区分）
result <- data.frame(model = c("naiveBayes", "C5.0", "CART", "ctree", "bagging", "boosting", "随机森林"),errTrain = rep(0, 7),  # 统一训练误差errTest = rep(0, 7)    # 统一测试误差
)result.10 <- data.frame(model = c("决策树", "随机森林", "人工神经网络"),errTrain = rep(0, 3),errTest = rep(0, 3)
)# 2. 10折交叉验证设置（保留原始配置）
# 当data不大时，训练集和测试集的划分情况会影响预测的准确性以及预测模型参数的选择
# so运用交叉验证生成多个超参组合从而判断出max 准确率 & 模型的参数
# n折表示分成n种情况
library(rpart)
control <- trainControl(method = "repeatedcv", number = 10, repeats = 3) # repeatedcv重复交叉验证# 3. 仅使用final_data训练模型 ----
# 在进行训练前，确保 Outcome 是因子（分类问题）而非数值。文档1中 Outcome 为二分类（0/1），需转换
final_data$Outcome <- as.factor(final_data$Outcome)
# 决策树
rpart.model <- train(Outcome ~ ., data = final_data, method = "rpart", trControl = control)
# 随机森林
rf.model <- train(Outcome ~ ., data = final_data, method = "rf", trControl = control)
# 神经网络
nnet.model <- train(Outcome ~ ., data = final_data, method = "nnet", trControl = control,preProcess = c("center", "scale"))  # 添加标准化

四、模型比较

# 4. 误差计算逻辑
# 决策树
result.10$errTrain[1] <- 1 - mean(predict(rpart.model, final_data) == final_data$Outcome)
result.10$errTest[1] <- 1 - max(rpart.model$results$Accuracy)  # 取最佳超参数结果
# 随机森林
result.10$errTrain[2] <- 1 - mean(predict(rf.model, final_data, type = "raw") == final_data$Outcome)
# 若需预测类别标签（如 Outcome=0/1），使用 type = "raw";若需预测概率（如 Outcome=1 的概率），使用 type = "prob"
result.10$errTest[2] <- 1 - max(rf.model$results$Accuracy)  
# 神经网络
# pred_train <- predict(nnet.model, final_data, type = "prob") # 神经网络输出为概率值，需阈值化
# pred_class <- ifelse(pred_train > 0.5, 1, 0)  # 概率转类别
# result.10$errTrain[3]  <- 1 - mean(pred_class == final_data$Outcome)result.10$errTrain[3] <- 1 - mean(predict(nnet.model, final_data) == final_data$Outcome)
result.10$errTest[3] <- 1 - max(nnet.model$results$Accuracy)  # 循环实现
for (i in 1:3){result.10$errTrain[i] <- 1 - mean(predict(switch(i,rpart.model,rf.model,nnet.model),final_data) == final_data$Outcome)
}for (i in 1:3){result.10$errTest[i] <- 1 - max(switch(i,rpart.model,rf.model,nnet.model)$results$Accuracy)
}# 使用多重插补法的final_data 进行建模ind <- createDataPartition(final_data$Outcome,times = 1,p=0.8,list = F)
train <- final_data[ind,] # 构建训练集
test <- final_data[-ind,] # 构建测试集# 查看目标变量的不平衡程度
prop.table(table(final_data$Outcome))
prop.table(table(train$Outcome))
prop.table(table(test$Outcome))### 建模和评估
# 使用naiveBayes函数建立朴素贝叶斯分类器
library(e1071)
str(final_data)
naiveBayes.model <- naiveBayes(Outcome~.,data = train) # 构建模型# 预测结果
train_predict <- predict(naiveBayes.model,newdata=train)
test_predict <- predict(naiveBayes.model,newdata=test)# 构建混淆矩阵
tableTrain <- table(actual=train$Outcome,predict=train_predict)
tableTest <- table(actual=test$Outcome,predict=test_predict)# 计算误差率
result[1,2] <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")
result[1,3] <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")# 决策树模型
library(C50)
C5.0.model <-C5.0(Outcome~.,data=train)library(rpart) 
rpart.model <-rpart(Outcome~.,data=train)library(party)
ctree.model <- ctree(Outcome~.,data=train)# 预测结果
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$Outcome,predict=train_predict)tableTest <- table(actual=test$Outcome,predict=test_predict)# 计算误差率result[i+1,2] <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")result[i+1,3] <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")
}# 使用集成学习算法Bagging和AdaBoost建模  
library(adabag)
bagging.model <- bagging(Outcome~.,data=train)
boosting.model <-boosting(Outcome~.,data=train)# 随机森林
library(randomForest)
randomForest.model <- randomForest(Outcome~.,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$Outcome,predict=switch(i,train_predict$class,train_predict$class,train_predict))tableTest <- table(actual=test$Outcome,predict=switch(i,test_predict$class,test_predict$class,test_predict))# 计算误差率result[i+4,2] <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")result[i+4,3] <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")
}# 利用rpart函数构建分类树
rpart.model1 <- rpart::rpart(Outcome~.,data = train,cp=0.01395349)
# rpart::rpart()：调用rpart包的rpart函数，用于构建分类或回归决策树。
# cp剪枝复杂度参数，R默认 .01，cp越小，树的分裂越细致（可能过拟合）；cp越大，树的结构越简单（可能欠拟合）。# 交叉验证法（recommand）求best cpmodel <- rpart(Outcome ~ ., data=train, cp=0)  # 生成完整树printcp(model)  # 查看交叉验证结果optimal_cp <- model$cptable[which.min(model$cptable[,"xerror"]), "CP"] # 0.01395349# 选择xerror最小的cp值，或xerror落在最小xerror + xstd（1倍标准误差）范围内的最简cp（即“1-SE规则”）# 利用randomForest函数构建随机森林
rf.model <-randomForest::randomForest(Outcome~.,data = train,mtry=7)
# mtry=2：表示每棵决策树在节点分裂时随机选择的特征数量为2。
# 默认值：分类问题通常为sqrt(p)（p为总特征数），回归问题为p/3。
# 影响：较小的mtry会增加树之间的差异性（减少过拟合），但可能牺牲精度；#     较大的mtry可能提高单棵树性能，但增加相关性（可能过拟合# 确定最优 mtry 值：rfcv()交叉验证 or 利用OOB误差 min# rfcv()
train_x <- train[, -which(names(train) == "Outcome")]  # 自变量数据
train_y <- train$Outcome # 目标变量
set.seed(123) 
rf_cv <- rfcv(train_x, train_y, cv.fold=5, scale="step", step=-1)  # 测试所有可能的mtry
optimal_mtry <- rf_cv$n.var[which.min(rf_cv$error.cv)] # 7
plot(rf_cv$n.var, rf_cv$error.cv, type = "b", xlab = "Number of Variables (mtry)", ylab = "Cross-Validation Error")
optimal_mtry# 利用nnet函数建立人工神经网络
nnet.model <- nnet::nnet(Outcome~.,data = train,size=2,decay=0.4)# size 隐藏层的神经元数量 ，nnet仅支持单隐藏层# decay 权重衰减参数（正则化项），用于防止过拟合，decay值越大，正则化强度越强（默认值为0，即无正则化）。# 选择 best 参数组（size & decay）：使用交叉验证 caret包自动化调参
library(caret)
grid <- expand.grid(size = 1:10, decay = seq(0, 1, by = 0.1))
ctrl <- trainControl(method = "cv", number = 5)  # 5折交叉验证
model <- train(Outcome ~ ., data = train, method = "nnet", trControl = ctrl, tuneGrid = grid, trace = FALSE)
print(model$bestTune)  # 输出最优参数 2   0.4# 预测结果
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$Outcome,predict=train_predict)tableTest <- table(actual=test$Outcome,predict=test_predict)# 计算误差率result.10[i,2] <-paste0(round((sum(tableTrain)-sum(diag(tableTrain)))*100/sum(tableTrain),2),"%")result.10[i,3] <-paste0(round((sum(tableTest)-sum(diag(tableTest)))*100/sum(tableTest),2),"%")
}