如何优雅地做logistic回归？

爱做科研的桶

于 2024-11-26 19:28:45 发布

阅读量867

点赞数 9

分类专栏： R语言数据分析文章标签：回归数据挖掘逻辑回归 r语言数据分析

本文链接：https://blog.youkuaiyun.com/llthxx/article/details/144054754

版权

R语言数据分析专栏收录该内容

8 篇文章

订阅专栏

什么是逻辑回归

逻辑回归（Logistic Regression）是一种广泛使用的分类算法，尤其在处理二分类问题时（例如，预测“是”或“否”结果时）非常有效。尽管名字中带有“回归”，它实际上是一种分类算法，用来预测离散的类别标签（通常是二分类问题），而不是连续的数值。它通过估计事件发生的概率，来进行预测。

1. 逻辑回归的基本概念

逻辑回归的核心思想是利用一个逻辑函数（或称sigmoid函数）来将线性回归的输出转换为概率值。通过这种方式，逻辑回归模型可以将预测值映射到 0 到 1 之间，这使得其非常适用于二分类问题。

2. 逻辑回归的作用

逻辑回归主要用于解决二分类问题，其作用如下：

二分类问题预测： 预测某个事件的概率，比如：一个邮件是否为垃圾邮件（是/否）、某个客户是否会购买某个产品（是/否）、某个病人是否患有某种疾病（是/否）等。
分类问题建模： 逻辑回归适用于任何二分类任务，特别是在处理存在多个自变量（特征）的情况下。
可解释性： 逻辑回归的模型系数 β\betaβ 具有很好的可解释性，可以了解每个特征对分类结果的影响。系数的符号和大小可以帮助我们理解变量之间的关系。

3. 逻辑回归的优缺点

优点：

简单且易于实现：逻辑回归的数学模型非常简单，易于理解和实现。
计算效率高：尤其在数据集不是非常庞大的情况下，训练速度非常快。
可解释性强：逻辑回归模型的系数可以帮助我们解释每个特征如何影响结果的概率。
概率输出：逻辑回归不仅仅给出分类结果，还能输出属于某个类别的概率，这对于决策制定非常有用。

缺点：

线性假设：逻辑回归假设特征与类别之间存在线性关系，因此对于特征之间的非线性关系处理能力较差。
容易过拟合：如果数据中的特征过多，或者数据量较少，模型容易产生过拟合。
不适用于多类别问题：虽然可以通过扩展（如 One-vs-Rest 或 Softmax回归）来处理多分类问题，但标准的逻辑回归仅适用于二分类问题。

方法一

1. 加载所需的库

library(openxlsx)
library(gtsummary)
library(dplyr)
library(tidyr)
library(car)
library(rms)
library(ResourceSelection)
library(pROC)
library(caret)
library(ggplot2)
library(lattice)

openxlsx: 用于读取和写入 Excel 文件。
gtsummary: 用于创建数据摘要表格，生成漂亮的总结统计。
dplyr: 数据处理的核心库，用于数据的筛选、分组、汇总等操作。
tidyr: 用于数据整理，特别是将数据从长格式转换为宽格式。
car: 用于计算多重共线性诊断（如 VIF）等。
rms: 用于拟合回归模型，进行模型评估等。
ResourceSelection: 用于执行Hosmer-Lemeshow拟合优度检验。
pROC: 用于绘制和分析ROC曲线。
caret: 用于分类模型的评价，特别是混淆矩阵和性能度量。
ggplot2 和 lattice: 用于数据可视化。

2. 读取数据并进行预处理

mydata <- read.xlsx("gpt.xlsx")  # 导入Excel数据
mydata <- subset(mydata, mydata$final != 3)
str(mydata)

read.xlsx("gpt.xlsx") 读取一个名为 gpt.xlsx 的 Excel 文件，并将数据存储在 mydata 数据框中。
subset(mydata, mydata$final != 3)：筛选掉 final 列中值为3的行。可以理解为剔除掉某些不符合条件的数据。
str(mydata)：查看数据的结构，输出数据框的每列的类型和预览。

3. 将一些变量转换为因子类型

mydata$gpt <- factor(mydata$gpt)
mydata$system <- factor(mydata$system)
mydata$query <- factor(mydata$query)
mydata$final_new <- factor(mydata$final, labels = c("不参加", "参加"))
mydata$sort_new <- factor(mydata$sort, labels = c("10%", "20%","30%","40%","50%","60%","70%","80%","90%"))
str(mydata)

将 gpt, system, query, final, 和 sort 列转换为因子（factor）类型，其中 final 列重命名为 final_new，并将其值重编码为 “不参加”和“参加”。
str(mydata)：再次检查数据的结构，确认因子的转换。

4. 重新设定参考水平

mydata$gpt <- relevel(mydata$gpt, ref = "gpt3.5")
mydata$system <- relevel(mydata$system, ref = "system2")
mydata$query <- relevel(mydata$query, ref = "query1")
mydata$sort_new <- relevel(mydata$sort_new, ref = "10%")

relevel()：将某些因子的参考水平重设为指定的类别。例如，将 gpt 因子的参考水平设为 “gpt3.5”。

5. 创建列联表

table(mydata$final_new)
table(mydata$final_new, mydata$gpt)
table(mydata$final_new, mydata$system)
table(mydata$final_new, mydata$query)
table(mydata$final_new, mydata$sort_new)

table()：生成列联表，查看 final_new 和其他变量之间的关系。列联表展示了每一对变量组合的频数。

6. 按多个变量分组并计算频数

result <- mydata %>%
  group_by(gpt, system, query, final_new, sort_new) %>%
  summarise(Frequency = n(), .groups = "drop") %>%  # 明确移除分组
  ungroup()
print(result)

使用 dplyr 包的管道操作符 %>%，将数据按多个变量（gpt, system, query, final_new, sort_new）分组，并计算每组的频数（n()）。
.groups = "drop"：在汇总后移除分组，使结果更加简洁。
ungroup()：确保数据框不再保持分组状态。

7. 将数据转换为宽格式

result_wide <- result %>%
  pivot_wider(names_from = final_new, values_from = Frequency, values_fill = list(Frequency = 0))

pivot_wider()：将 final_new 变量转换为列，并填充每组频数。values_fill = list(Frequency = 0) 表示如果某个组合的频数为空，则用 0 填充。

8. 计算参加的比例

result_wide <- result_wide %>%
  mutate(Proportion = 参加 / (不参加 + 参加))
print(result_wide)

mutate()：创建一个新的列 Proportion，表示每组参加的比例。

9. 单因素逻辑回归模型

fit1 <- glm(final_new ~ gpt, data = mydata, family = "binomial")
fit2 <- glm(final_new ~ system, data = mydata, family = "binomial")
fit3 <- glm(final_new ~ query, data = mydata, family = "binomial")
fit4 <- glm(final_new ~ sort_new, data = mydata, family = "binomial")

使用 glm() 拟合逻辑回归模型，每个模型只使用一个自变量（gpt, system, query, sort_new）。
family = "binomial"：表示这是一个二项逻辑回归模型。

10. 查看回归模型的结果

summary(fit1)
summary(fit2)
summary(fit3)
summary(fit4)

summary()：显示每个逻辑回归模型的详细结果，包括回归系数、标准误差、Z 值、P 值等。

11. 多变量逻辑回归模型

fita <- glm(final_new ~ gpt + system + query + sort_new, data = mydata, family = "binomial")
fitb <- glm(final_new ~ gpt + system + query + sort_new, data = mydata, family = "quasibinomial")

fita：拟合一个标准的逻辑回归模型。
fitb：拟合一个准二项式（quasibinomial）回归模型，用于处理过度离散的问题。

12. 检查过度离散（Overdispersion）

pchisq(summary(fitb)$dispersion * fita$df.residual,
       fita$df.residual, lower = F)

通过卡方检验检查过度离散问题，如果 p < 0.05，则说明数据可能存在过度离散。

13. Cook距离与多重共线性诊断

cook <- cooks.distance(fitb)
cook[cooks.distance(fitb) > 0.5]  # 显示Cook距离大于0.5的个案
max(cook)  # 显示最大Cook距离

vif(fitb)  # 计算方差膨胀因子（VIF）

Cook距离：用来检测回归模型中的异常值或影响力较大的观测值。如果 Cook 距离大于 0.5，通常认为该点可能是异常值。
VIF：用于检测多重共线性，如果 VIF 值大于 10，说明可能存在多重共线性问题。

14. 模型拟合与结果解释

summary(fitb)
confint(fitb)
exp(coef(fitb))  # 计算OR
exp(confint(fitb))  # 计算OR的95%CI

summary(fitb)：显示回归模型的摘要结果。
confint(fitb)：计算回归系数的置信区间。
exp(coef(fitb))：将回归系数转化为优势比（Odds Ratio, OR），并计算置信区间。

15. 模型拟合优度检验

hoslem.test(fitb$y, fitted(fitb), g = 10)
deviance(fitb) / df.residual(fitb)

Hosmer-Lemeshow检验：用于检验模型的拟合优度，p > 0.05 表示模型拟合良好。
过度离散检验：检查模型的 deviance 和残差自由度的比值，接近 1 表示没有过度离散。

16. ROC 曲线

mydata$pre <- predict(fitb, type = 'response')
modelroc <- roc(mydata$final_new, mydata$pre)
plot(modelroc, print.auc = TRUE, auc.polygon = TRUE, grid = c(0.1, 0.2), grid.col = c("green", "red"))

计算模型的预测概率并绘制 ROC 曲线，auc 是曲线下面积，反映模型的分类能力。

17. 混淆矩阵与模型评估

mydata$pred.cutoff <- ifelse(mydata$pre >= 0.836, "参加", "不参加")
confusionMatrix(data = factor(mydata$pred.cutoff), reference = factor(mydata$final_new), positive = "参加")

根据预测概率确定分类的截断值（0.836），并计算混淆矩阵，评估模型的性能（准确率、灵敏度、特异度等）。

总结

# 1. 加载所需的库
library(openxlsx)
library(gtsummary)
library(dplyr)
library(tidyr)
library(car)
library(rms)
library(ResourceSelection)
library(pROC)
library(caret)
library(ggplot2)
library(lattice)

# 2. 读取数据并进行预处理
mydata <- read.xlsx("gpt.xlsx")  # 导入Excel数据
mydata <- subset(mydata, mydata$final != 3)  # 筛选掉 final 列为 3 的数据
str(mydata)

# 3. 将一些变量转换为因子类型
mydata$gpt <- factor(mydata$gpt)
mydata$system <- factor(mydata$system)
mydata$query <- factor(mydata$query)
mydata$final_new <- factor(mydata$final, labels = c("不参加", "参加"))
mydata$sort_new <- factor(mydata$sort, labels = c("10%", "20%","30%","40%","50%","60%","70%","80%","90%"))
str(mydata)

# 4. 重新设定参考水平
mydata$gpt <- relevel(mydata$gpt, ref = "gpt3.5")
mydata$system <- relevel(mydata$system, ref = "system2")
mydata$query <- relevel(mydata$query, ref = "query1")
mydata$sort_new <- relevel(mydata$sort_new, ref = "10%")

# 5. 创建列联表
table(mydata$final_new)
table(mydata$final_new, mydata$gpt)
table(mydata$final_new, mydata$system)
table(mydata$final_new, mydata$query)
table(mydata$final_new, mydata$sort_new)

# 6. 按多个变量分组并计算频数
result <- mydata %>%
  group_by(gpt, system, query, final_new, sort_new) %>%
  summarise(Frequency = n(), .groups = "drop") %>%  # 明确移除分组
  ungroup()
print(result)

# 7. 将数据转换为宽格式
result_wide <- result %>%
  pivot_wider(names_from = final_new, values_from = Frequency, values_fill = list(Frequency = 0))

# 8. 计算参加的比例
result_wide <- result_wide %>%
  mutate(Proportion = 参加 / (不参加 + 参加))
print(result_wide)

# 9. 单因素逻辑回归模型
fit1 <- glm(final_new ~ gpt, data = mydata, family = "binomial")
fit2 <- glm(final_new ~ system, data = mydata, family = "binomial")
fit3 <- glm(final_new ~ query, data = mydata, family = "binomial")
fit4 <- glm(final_new ~ sort_new, data = mydata, family = "binomial")

# 10. 查看回归模型的结果
summary(fit1)
summary(fit2)
summary(fit3)
summary(fit4)

# 11. 多变量逻辑回归模型
fita <- glm(final_new ~ gpt + system + query + sort_new, data = mydata, family = "binomial")
fitb <- glm(final_new ~ gpt + system + query + sort_new, data = mydata, family = "quasibinomial")

# 12. 检查过度离散（Overdispersion）
pchisq(summary(fitb)$dispersion * fita$df.residual,
       fita$df.residual, lower = F)

# 13. 检查 Cook 距离与多重共线性诊断
cook <- cooks.distance(fitb)  # 计算 Cook 距离
cook[cooks.distance(fitb) > 0.5]  # 显示 Cook 距离大于 0.5 的个案
max(cook)  # 显示最大 Cook 距离

# 共线性诊断（方差膨胀因子）
vif(fitb)

# 14. 模型拟合与结果解释
summary(fitb)
confint(fitb)  # 计算回归系数的置信区间
exp(coef(fitb))  # 计算优势比（OR）
exp(confint(fitb))  # 计算优势比（OR）的 95% 置信区间

# 15. 模型拟合优度检验
hoslem.test(fitb$y, fitted(fitb), g = 10)  # Hosmer-Lemeshow 拟合优度检验
deviance(fitb) / df.residual(fitb)  # 过度离散检验

# 16. 计算模型预测值并绘制 ROC 曲线
mydata$pre <- predict(fitb, type = 'response')  # 计算预测概率
modelroc <- roc(mydata$final_new, mydata$pre)  # 计算 ROC 数据
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 = "skyblue", print.thres = TRUE)  # 绘制 ROC 曲线
plot(ci(modelroc, of = "thresholds", thresholds = "best"))  # 添加最佳截断值位置
ci(modelroc)  # 计算 AUC 的 95% CI

# 17. 模型评价的混淆矩阵
mydata$pred.cutoff <- ifelse(mydata$pre >= 0.836, "参加", "不参加")  # 由最佳截断值判断分类
confusionMatrix(data = factor(mydata$pred.cutoff), 
                reference = factor(mydata$final_new), 
                positive = "参加")  # 混淆矩阵与模型评估

# 输出指标：accuracy, sensitivity, specificity, pos, neg 等

方法二

1. 选择分析变量

logreg_uv <- mydata %>%
  select(gpt, system, query, sort_new, final_new) %>%

select(gpt, system, query, sort_new, final_new)：使用 dplyr 包的 select() 函数从 mydata 数据框中选择变量 gpt、system、query、sort_new 和 final_new。这些变量是后续回归分析中要使用的预测变量和因变量。
结果是一个包含这些变量的子数据框，并且会传递到下一个管道操作 %>%。

2. 进行单变量回归分析

tbl_uvregression(
    method = glm, # 设置回归模型
    y = final_new, # 设置因变量
    method.args = list(family = quasibinomial(link='logit')), # 设置glm参数
    exponentiate = TRUE, # 将系数进行自然指数转换
    estimate_fun = function(x){style_ratio(x, digits = 3)}, # 设置显示小数点位数
    pvalue_fun = function(x){style_pvalue(x, digits = 3)}, # 设置p值的小数点位数
    hide_n = TRUE # 不显示计数
  ) %>%

tbl_uvregression()：gtsummary 包中的一个函数，用于执行单变量回归分析并生成结果汇总表。在这里，我们进行的是 单变量逻辑回归（每个预测变量分别与因变量 final_new 进行回归）。

具体参数解释：
- method = glm：指定回归模型为广义线性模型（Generalized Linear Model，glm）。这里使用的是逻辑回归模型，适用于二项型因变量（final_new）的分析。
- y = final_new：指定因变量是 final_new，即回归的目标变量（参加/不参加）。
- method.args = list(family = quasibinomial(link='logit'))：设置 glm 函数的参数，指定使用 广义线性模型，且适用于二项分布（logistic回归）。quasibinomial(link='logit') 表示使用 quasibinomial 分布，并且连接函数为 logit（即 logistic 回归）。
- exponentiate = TRUE：回归系数默认以对数尺度显示，exponentiate = TRUE 会对回归系数进行自然指数转换，输出 优势比 (Odds Ratios, OR)，它更加容易解释（表示为每单位增加对应的变化倍数）。
- estimate_fun = function(x){style_ratio(x, digits = 3)}：定义回归系数的显示格式，通过 style_ratio() 函数格式化显示为比例（OR），并指定保留 3 位小数。
- pvalue_fun = function(x){style_pvalue(x, digits = 3)}：定义 p 值的显示格式，通过 style_pvalue() 函数格式化 p 值，保留 3 位小数。
- hide_n = TRUE：不显示每个变量的计数信息（n）。

3. 添加全局 p 值和加粗显著 p 值

  add_global_p() %>%
  bold_p()

add_global_p()：为表格添加一个全局 p 值，通常用于检验回归模型的总体显著性（如模型是否在总体上显著）。
bold_p()：将 p 值小于指定显著性水平（通常为 0.05）的回归结果加粗，以便突出显示显著的回归系数。

4. 显示结果

logreg_uv

最后，代码输出 logreg_uv 对象，显示单变量回归分析的结果表格，结果包括每个预测变量（gpt, system, query, sort_new）与因变量 final_new 的回归系数、优势比 (OR)、p 值等。回归系数如果为 0.1，则意味着相应的变量每增加一单位，final_new 为 "参加" 的可能性会增加或减少 10%（如果 exponentiate = TRUE 的话）。

总结

library(openxlsx)
library(gtsummary)
library(dplyr)

# 导入数据
mydata <- read.xlsx("gpt.xlsx")  # 导入Excel文件
mydata <- subset(mydata, mydata$final != 3)  # 过滤数据，去除final等于3的行
str(mydata)

# 处理数据类型
mydata$gpt <- factor(mydata$gpt)
mydata$system <- factor(mydata$system)
mydata$query <- factor(mydata$query)
mydata$final_new <- factor(mydata$final, labels = c("不参加", "参加"))
mydata$sort_new <- factor(mydata$sort, labels = c("10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%"))
str(mydata)

# 设置参考水平
mydata$gpt <- relevel(mydata$gpt, ref = "gpt3.5")
mydata$system <- relevel(mydata$system, ref = "system2")
mydata$query <- relevel(mydata$query, ref = "query1")
mydata$sort_new <- relevel(mydata$sort_new, ref = "10%")

# 单变量回归分析（逻辑回归）
logreg_uv <- mydata %>%
  select(gpt, system, query, sort_new, final_new) %>%
  tbl_uvregression(
    method = glm, # 设置回归模型
    y = final_new, # 设置因变量
    method.args = list(family = quasibinomial(link='logit')), # 设置glm参数
    exponentiate = TRUE, # 将系数进行自然指数转换
    estimate_fun = function(x){style_ratio(x, digits = 3)}, # 设置显示小数点位数
    pvalue_fun = function(x){style_pvalue(x, digits = 3)}, # 设置p值的小数点位数
    hide_n = TRUE # 不显示计数
  ) %>%
  add_global_p() %>%  # 添加全局p值
  bold_p()  # 显著的p值加粗

# 显示结果
logreg_uv