R语言中为每个rs开头列批量拟合两种逻辑回归模型的方法

需求说明 #

针对数据集中所有以rs开头的列，分别拟合以下两种逻辑回归模型：

• 模型1：AD ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + [rs列]（依次替换每个rs列）
• 模型2：AD ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + [rs列] * alcohol_intake（依次替换每个rs列）

修正后的示例数据集 #

原示例代码存在行长度不匹配问题，以下是可直接运行的修正版：

set.seed(543)
mydata <- data.frame(ID = 1:500, 
                     sex = sample(1:2, size = 500, replace = TRUE), 
                     age = runif(500, min= 35, max = 70), 
                     bmi = runif(500, min= 15, max = 35),
                     smoker = rep(c('smoker', 'not smoker'), 250),
                     alcohol_intake = rep(c('regular', 'not regular'), 250),
                     AD = sample(0:1, size = 500, replace = TRUE),
                     cov1 = runif(500, min= 0.0000, max = 1.0000),
                     cov2 = runif(500, min= -3.013e-02, max = 1.818e-02),
                     cov3 = runif(500, min= -3.562e-02, max = 1.540e-02),
                     cov4 = runif(500, min= -2.356e-02, max = 1.685e-02),
                     cov5 = runif(500, min= -1.392e-02, max = 2.894e-02),
                     cov6 = runif(500, min= -1.896e-02, max = 2.136e-02),
                     rs1 = sample(0:2, size = 500, replace = TRUE), 
                     rs2 = sample(0:2, size = 500, replace = TRUE),
                     rs3 = sample(0:2, size = 500, replace = TRUE),
                     rs4 = sample(0:2, size = 500, replace = TRUE),
                     rs5 = sample(0:2, size = 500, replace = TRUE),
                     rs6 = sample(0:2, size = 500, replace = TRUE),
                     rs7 = sample(0:2, size = 500, replace = TRUE),
                     rs8 = sample(0:2, size = 500, replace = TRUE),
                     rs9 = sample(0:2, size = 500, replace = TRUE),
                     rs10 = sample(0:2, size = 500, replace = TRUE),
                     rs11 = sample(0:2, size = 500, replace = TRUE),
                     rs12 = sample(0:2, size = 500, replace = TRUE),
                     rs13 = sample(0:2, size = 500, replace = TRUE),
                     rs14 = sample(0:2, size = 500, replace = TRUE),
                     rs15 = sample(0:2, size = 500, replace = TRUE),
                     rs16 = sample(0:2, size = 500, replace = TRUE),
                     rs17 = sample(0:2, size = 500, replace = TRUE),
                     rs18 = sample(0:2, size = 500, replace = TRUE),
                     rs19 = sample(0:2, size = 500, replace = TRUE),
                     rs20 = sample(0:2, size = 500, replace = TRUE),
                     rs21 = sample(0:2, size = 500, replace = TRUE),
                     rs22 = sample(0:2, size = 500, replace = TRUE),
                     rs23 = sample(0:2, size = 500, replace = TRUE),
                     rs24 = sample(0:2, size = 500, replace = TRUE),
                     rs25 = sample(0:2, size = 500, replace = TRUE),
                     rs26 = sample(0:2, size = 500, replace = TRUE),
                     rs27 = sample(0:2, size = 500, replace = TRUE),
                     rs28 = sample(0:2, size = 500, replace = TRUE),
                     rs29 = sample(0:2, size = 500, replace = TRUE),
                     rs30 = sample(0:2, size = 500, replace = TRUE)
)

解决方案 #

1. 提取所有rs开头的列名 #

先筛选出目标列：

# 基础R方法
rs_cols <- grep("^rs", colnames(mydata), value = TRUE)

# 可选：用stringr包（需先安装install.packages("stringr")）
# library(stringr)
# rs_cols <- str_subset(colnames(mydata), "^rs")

2. 用sapply批量拟合模型 #

模型1：无交互项

# 存储所有模型结果
model_list1 <- sapply(rs_cols, function(col) {
  formula <- as.formula(paste0("AD ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + ", col))
  glm(formula, data = mydata, family = binomial(link = "logit"))
}, simplify = FALSE)

# 查看单个模型结果，比如rs1的模型
summary(model_list1[["rs1"]])

模型2：带交互项

model_list2 <- sapply(rs_cols, function(col) {
  formula <- as.formula(paste0("AD ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + ", col, " * alcohol_intake"))
  glm(formula, data = mydata, family = binomial(link = "logit"))
}, simplify = FALSE)

# 查看rs1的交互模型结果
summary(model_list2[["rs1"]])

3. 用for循环批量拟合模型 #

如果更习惯循环写法：

# 初始化空列表存储模型
model_list1_loop <- list()
model_list2_loop <- list()

for (col in rs_cols) {
  # 模型1
  formula1 <- as.formula(paste0("AD ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + ", col))
  model_list1_loop[[col]] <- glm(formula1, data = mydata, family = binomial)
  
  # 模型2
  formula2 <- as.formula(paste0("AD ~ cov1 + cov2 + cov3 + cov4 + cov5 + cov6 + ", col, " * alcohol_intake"))
  model_list2_loop[[col]] <- glm(formula2, data = mydata, family = binomial)
}

4. 提取关键结果（OR值和p值） #

快速提取每个模型中rs列的核心统计量，以模型1为例：

# 定义提取函数
extract_results <- function(model, rs_col) {
  coefs <- summary(model)$coefficients
  rs_row <- grep(rs_col, rownames(coefs))
  or <- exp(coefs[rs_row, 1])
  p_val <- coefs[rs_row, 4]
  data.frame(rs_col = rs_col, OR = round(or, 3), p_value = round(p_val, 4))
}

# 批量提取模型1的结果
results1 <- do.call(rbind, lapply(names(model_list1), function(col) {
  extract_results(model_list1[[col]], col)
}))

print(results1)

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