R语言实现逻辑回归(LR)以及绘制ROC曲线和混淆矩阵

逻辑回归（Logistic Regression）是一种在统计学中处理因变量多分类问题的回归模型，其在生物信息学领域中有广泛的应用。今天分享一下如何用R实现逻辑回归模型，并绘制ROC曲线和混淆矩阵。

逻辑回归(LR)用途

分类问题：逻辑回归可以用于预测生物样本的分类问题，如疾病分类、肿瘤分型、蛋白质功能预测等。通过使用已知的特征和相应的标签信息，逻辑回归模型可以根据样本的特征来预测其所属的类别。
生物标志物预测：逻辑回归可以用于预测和鉴定生物标志物（biomarker），例如通过基因表达数据预测特定疾病的生物标志物。逻辑回归模型可以帮助确定哪些基因或特征与特定生物过程或疾病相关，并从中选择最有价值的标志物。
基因表达数据分析：逻辑回归可用于分析基因表达数据中的关联和预测。通过使用基因表达水平作为特征，可以预测样本的某些特性，如治疗反应、细胞类型、生物过程等。
DNA序列分析：逻辑回归可以用于分析DNA序列数据，如DNA序列的功能预测、转录因子结合位点预测、DNA修复位点预测等。通过使用序列特征和相应的功能或结合信息，逻辑回归模型可以预测DNA序列的不同功能或相互作用。

R语言实现LR

测试数据和代码

# 导入必要的包,没有安装的可以先安装一些

library(dplyr) #数据处理使用

library(data.table) #数据读取使用

library(caTools) #LR模型使用

library(pROC) #绘图使用

library(ggplot2) #绘图使用

library(ggpubr) #绘图使用

library(ggprism) #绘图使用

# 读取数据

data <- fread("./ML_data.txt",data.table = F)  # 替换为你的数据文件名或路径

数据长这个样子，一共35727行，214列。每一行代表一个样本，第一列是样本标签”cluster1或cluster2″，后面213列是213个特征。我们想根据213个特征，使用LR模型训练出一个能够对样本进行精准分类。

数据情况

构建LR模型

# 将第一列数据转换成1或0，分别代表cluster1和cluster2

data <- data %>% mutate(cluster = ifelse(cluster=="cluster1",1,0))

# 分割数据为训练集和测试集

set.seed(123)  # 设置随机种子，保证结果可复现

split <- sample.split(data$cluster, SplitRatio = 0.8)  # 将数据按照指定比例分割

train_data <- subset(data, split == TRUE)  # 训练集

test_data <- subset(data, split == FALSE)  # 测试集



# 构建逻辑回归模型

model <- glm(cluster ~ ., data = train_data, family = binomial)

# 查看模型摘要

summary(model)

在逻辑回归模型中，可以通过查看模型摘要来获取每个特征的重要性指标。逻辑回归模型的摘要通常提供了特征的系数（coefficient）或权重（weight），用于衡量每个特征对预测的贡献程度。

模型摘要

预测

# 预测

predictions <- predict(model, newdata = test_data, cluster = "response")



# 将预测结果转化为二分类（0和1）

threshold <- 0.5  # 设置阈值

predicted_classes <- ifelse(predictions >= threshold, 1, 0)



# 模型评估

confusion_matrix <- table(test_data$cluster, predicted_classes)

accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)

precision <- confusion_matrix[2, 2] / sum(confusion_matrix[, 2])

recall <- confusion_matrix[2, 2] / sum(confusion_matrix[2, ])

f1_score <- 2 * precision * recall / (precision + recall)



# 打印评估指标

print(confusion_matrix)

cat("Accuracy: ", accuracy, "\n")

cat("Precision: ", precision, "\n")

cat("Recall: ", recall, "\n")

cat("F1 Score: ", f1_score, "\n")

混淆矩阵

模型评价指标

confusion_matrix：混淆矩阵是用于评估分类模型性能的重要工具。它是一个二维矩阵，用于表示实际类别与预测类别之间的对应关系。
accuracy：准确率（Accuracy）是分类模型最常用的评估指标之一，它表示模型预测正确的样本占总样本数的比例。
precision：精确率（Precision）是衡量模型在预测为正类别中的准确性的指标。它表示预测为正类别的样本中真实为正类别的比例。在这里，precision 的计算公式是真正例（真实为正类别且预测为正类别）的样本数除以所有预测为正类别的样本数。
recall：召回率（Recall），也称为灵敏度（Sensitivity）或真正例率（True Positive Rate），衡量模型在正类别样本中能够正确识别的比例。它表示真实为正类别的样本中被正确预测为正类别的比例。在这里，recall 的计算公式是真正例的样本数除以所有真实为正类别的样本数。
f1_score：F1分数是综合考虑了精确率和召回率的综合指标，它是精确率和召回率的调和平均值。F1分数在不同类别不平衡或有偏斜的情况下更具有稳定性。

从这些指标来看，我们构建的LR模型可以很好的对样本进行分类，下面我们来绘制一下ROC曲线和混淆矩阵图。

R语言绘制ROC曲线

ROC（Receiver Operating Characteristic）曲线和AUC（Area Under the Curve）是评估二分类模型性能常用的指标。

ROC曲线是以模型的真阳性率（True Positive Rate，也称为灵敏度或召回率）为纵坐标，以模型的假阳性率（False Positive Rate）为横坐标所绘制的曲线。ROC曲线能够综合考虑分类模型在不同阈值下的预测能力。
AUC是ROC曲线下的面积，代表模型的整体性能。AUC值范围在0.5到1之间，其中0.5表示模型的预测性能等同于随机预测，而1表示模型的预测性能完美。通常来说，AUC值越接近1，模型的性能越好。

# 计算ROC曲线指标

roc_obj <- roc(test_data$cluster, predictions)

roc_auc <- auc(roc_obj)



# 将ROC对象转换为数据框

roc_data <- data.frame(1 - roc_obj$specificities, roc_obj$sensitivities)



# 绘制ROC曲线

ggplot(roc_data, aes(x = 1 - roc_obj$specificities, y = roc_obj$sensitivities)) +

  geom_line(color = "#0073C2FF", size = 1.5) +

  geom_segment(aes(x = 0, y = 0, xend = 1, yend = 1), linetype = "dashed", color = "gray") +

  geom_text(aes(x = 0.8, y = 0.2, label = paste("AUC =", round(roc_auc, 2))), size = 4, color = "black") +

  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +

  theme_pubr() +

  labs(x = "1 - Specificity", y = "Sensitivity") +

  ggtitle("ROC Curve") +

  theme(plot.title = element_text(size = 14, face = "bold"))+

  theme_prism(border = T)

R语言绘制混淆矩阵

混淆矩阵提供了对分类模型性能的全面评估。它展示了实际类别和预测类别之间的对应关系，可以清晰地看到模型的预测结果中真正例、真反例、假正例和假反例的数量或比例。绘制混淆矩阵可以将复杂的分类结果以直观的方式展示出来，使得结果更易于理解和解释。

# 创建混淆矩阵数据

confusion_matrix <- table(test_data$cluster, predicted_classes)



# 将混淆矩阵转换为数据框

confusion_matrix_df <- as.data.frame.matrix(confusion_matrix)

colnames(confusion_matrix_df) <- c("cluster1","cluster2")

rownames(confusion_matrix_df) <- c("cluster1","cluster2")

## 这里计算的是比例，也可以使用具体数量，个人感觉比例能够直观的展示

draw_data <- round(confusion_matrix_df / rowSums(confusion_matrix_df),2)

draw_data$real <- rownames(draw_data)

draw_data <- melt(draw_data)

# 绘制混淆矩阵热图

ggplot(draw_data, aes(real,variable, fill = value)) +

  geom_tile() +

  geom_text(aes(label = scales::percent(value))) +

  scale_fill_gradient(low = "#F0F0F0", high = "#3575b5") +

  labs(x = "True", y = "Guess", title = "Confusion matrix") +

  theme_prism(border = T)+

  theme(panel.border = element_blank(),

        axis.ticks.y = element_blank(),

        axis.ticks.x = element_blank(),

        legend.position="none")

文章转自微信公众号@Bio小菜鸟