在R循环中为逻辑回归结果添加变量各水平计数

需求说明 #

通过循环批量执行逻辑回归，提取汇总统计量（OR值、置信区间、P值），同时自动添加每个自变量各水平的样本计数，最终输出到Excel文件，避免手动统计计数的重复工作。

现有代码基础 #

用户已实现批量逻辑回归及结果整理，但缺少计数列，现有代码如下：

library(tidyverse)
install.packages("AER")
library("AER")
data(Affairs, package="AER")
Affairs$ynaffair[Affairs$affairs >  0] <- 1
Affairs$ynaffair[Affairs$affairs == 0] <- 0

Affairs <- Affairs %>% 
  mutate_at(c("affairs", "religiousness", "occupation", "rating", "ynaffair"), as.factor) 

frmlas <- list(ynaffair~gender,
               ynaffair~children,
               ynaffair~religiousness,
               ynaffair~occupation,
               ynaffair~rating)
output <- list() 
output_df_list <- list() 

for(i in 1:length(frmlas)){
  output[[i]] <- glm(frmlas[[i]], data = Affairs, family=binomial)
  names(output)[i] <- capture.output(frmlas[[i]]) 
  output_df_list[[i]] <- data.frame("Model"=capture.output(frmlas[[i]]),
                                    "Covariate"=names(output[[i]]$coefficients),
                                    "Beta_Estimate"=output[[i]]$coefficients,
                                    "P_val"=summary(output[[i]])$coefficients[,4],
                                    "CI_95_LL"=confint(output[[i]])[,1],
                                    "CI_95_UL"=confint.default(output[[i]])[,2])
    rownames(output_df_list[[i]]) <- NULL             
}  

output_df_full_t1 <- do.call("rbind", output_df_list)

output_df_full_t1 <- output_df_full_t1 %>% 
  mutate(OR = exp(Beta_Estimate)) %>% 
  mutate(CI_95_LL = exp(CI_95_LL),
         CI_95_UL = exp(CI_95_UL))  %>% 
  select(Model, Covariate, OR, CI_95_LL, CI_95_UL, P_val) %>% 
  filter(!(Covariate %in% '(Intercept)'))  %>% 
  mutate_if(is.numeric, round, digits = 3)   %>% 
  unite(CI, c(CI_95_LL, CI_95_UL), sep = ", ", remove = TRUE) 

head(output_df_full_t1)

修改方案：添加变量水平计数 #

核心思路是在循环中，针对每个模型的自变量，统计其各水平的样本量，再合并到结果数据框中。

修改后的完整代码 #

library(tidyverse)
install.packages(c("AER", "writexl")) # 新增writexl用于输出Excel
library("AER")
library("writexl")

# 数据预处理
data(Affairs, package="AER")
Affairs$ynaffair <- as.factor(ifelse(Affairs$affairs > 0, 1, 0))
Affairs <- Affairs %>% 
  mutate_at(c("religiousness", "occupation", "rating"), as.factor) 

# 定义模型公式列表
frmlas <- list(ynaffair~gender,
               ynaffair~children,
               ynaffair~religiousness,
               ynaffair~occupation,
               ynaffair~rating)

output_df_list <- list() 

for(i in 1:length(frmlas)){
  # 拟合逻辑回归模型
  model <- glm(frmlas[[i]], data = Affairs, family=binomial)
  model_formula <- capture.output(frmlas[[i]])
  
  # 提取自变量名称
  covar_name <- as.character(frmlas[[i]][[3]])
  
  # 统计自变量各水平的样本计数
  count_df <- Affairs %>% 
    count(.data[[covar_name]], name = "Count") %>% 
    rename(Covariate = !!sym(covar_name))
  
  # 提取模型结果并整理
  model_result <- data.frame(
    Model = model_formula,
    Covariate = names(model$coefficients),
    Beta_Estimate = model$coefficients,
    P_val = summary(model)$coefficients[,4],
    CI_95_LL = confint(model)[,1],
    CI_95_UL = confint.default(model)[,2]
  ) %>% 
    filter(Covariate != "(Intercept)") %>% # 去掉截距项
    left_join(count_df, by = "Covariate") # 合并计数列
  
  # 加入结果列表
  output_df_list[[i]] <- model_result
}  

# 合并所有结果并整理格式
output_df_full <- do.call("rbind", output_df_list) %>% 
  mutate(OR = exp(Beta_Estimate)) %>% 
  mutate(CI_95_LL = exp(CI_95_LL),
         CI_95_UL = exp(CI_95_UL)) %>% 
  select(Model, Covariate, Count, OR, CI_95_LL, CI_95_UL, P_val) %>% 
  mutate_if(is.numeric, round, digits = 3) %>% 
  unite(CI, c(CI_95_LL, CI_95_UL), sep = ", ", remove = TRUE)

# 查看结果
head(output_df_full)

# 输出到Excel文件
write_xlsx(output_df_full, "逻辑回归结果带计数.xlsx")

关键修改点说明 #

1. 统计计数：在循环中通过count()函数统计当前自变量各水平的样本量，生成count_df数据框。
2. 合并计数：使用left_join()将计数列合并到模型结果中，确保每个自变量水平对应正确的样本数。
3. 输出Excel：使用writexl包的write_xlsx()函数直接将结果导出到Excel，无需手动复制。
4. 简化代码：去掉了不必要的output列表，直接在循环中整理结果，代码更简洁。

内容的提问来源于stack exchange，提问作者Newtostats_24