在R的Cox风险比函数中添加控制变量job_site的技术求助

核心报错原因：加入job_site控制变量后，每个Cox模型包含两个变量，但原代码仍用单变量模型的固定索引提取结果，导致维度混乱，触发行数不匹配错误。另外，公式中直接用breakthrough$引用列的写法易引发环境问题。

修改后的完整代码 #

covariates <- c("outpatient_settings",
                "icu",
                "non_icu",
                "emergency_department",
                "sex")

# 构造带控制变量的公式，用列名替代直接数据框引用
univ_formulas <- sapply(covariates, function(x) {
  as.formula(paste('Surv(person_time, n_result) ~ job_site +', x))
})

# 批量拟合Cox模型
univ_models <- lapply(univ_formulas, function(x) {
  coxph(x, data = breakthrough)
})

# 精准提取每个目标协变量的结果
univ_results <- Map(function(model, target_var) {
  summ <- summary(model)
  # 提取目标变量的beta、HR及置信区间
  beta <- signif(summ$coef[target_var, "coef"], digits = 2)
  hr <- signif(summ$coef[target_var, "exp(coef)"], digits = 2)
  hr_ci_lower <- signif(summ$confint[target_var, "lower .95"], digits = 2)
  hr_ci_upper <- signif(summ$confint[target_var, "upper .95"], digits = 2)
  hr_with_ci <- paste0(hr, " (", hr_ci_lower, "-", hr_ci_upper, ")")
  
  # 提取模型整体wald检验统计量和p值（与原逻辑保持一致）
  wald_test <- signif(summ$wald["test"], digits = 2)
  p_value <- signif(summ$wald["pvalue"], digits = 2)
  
  res <- c(beta, hr_with_ci, wald_test, p_value)
  names(res) <- c("beta", "HR (95% CI for HR)", "wald.test", "p.value")
  return(res)
}, univ_models, covariates)

# 整理为最终数据框
res <- t(as.data.frame(univ_results, check.names = TRUE))
as.data.frame(res)

关键修改说明 #

• 公式优化：移除breakthrough$前缀，依赖coxph的data参数读取数据，避免环境变量冲突；固定加入job_site作为控制变量。
• 结果提取逻辑修正：用Map将模型与目标协变量一一对应，通过变量名（而非固定索引）提取结果，确保拿到的是目标协变量的统计值，而非控制变量的结果。
• 维度一致性保障：统一每个模型的结果提取规则，保证所有结果向量长度一致，解决转置时的行数不匹配问题。

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