在R中复现Stata的Cox生存模型结果出现差异，如何解决？

Stata原代码 #

stset t, id(leadid) failure(c_coup)

stcox legislature  lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age

R复现代码 #

# Load survival package
library(survival)

# Set the survival object
surv_obj <- Surv(data$t, data$c_coup)

# Run model
m1 <- coxph(surv_obj ~ legislature + lgdp_1 + growth_1 + exportersoffuelsmainlyoil_EL2008 + ethfrac_FIXED + communist + mil + cw + age, data = data, method = "breslow")

# Examine hazard ratios
exp(coef(m1))

差异现象 #

Stata中legislature的风险比（HR）估计值为0.298，但R中得到的对应值为0.1688371，结果无法匹配。

数据集预览 #

structure(list(t = structure(c(1, 2, 3, 4, 5, 6), label = "Current time in office", format.stata = "%9.0g"), 
    c_coup = structure(c(0, 0, 0, 0, 0, 0), format.stata = "%9.0g"), 
    leadid = structure(c("A2.2-208", "A2.2-208", "A2.2-208", 
    "A2.2-208", "A2.2-208", "A2.2-208"), label = "Leader ID", format.stata = "%13s"), 
    legislature = structure(c(1, 1, 1, 1, 1, 1), format.stata = "%9.0g"), 
    lgdp_1 = structure(c(7.68524360656738, 7.69938945770264, 
    7.54960918426514, 7.57916784286499, 7.6033992767334, 7.67089462280273
    ), format.stata = "%9.0g"), growth_1 = structure(c(6.35386085510254, 
    1.42463231086731, -13.910285949707, 3, 2.45273375511169, 
    6.98254346847534), label = "annual growth, t-1, Maddison", format.stata = "%9.0g"), 
    exportersoffuelsmainlyoil_EL2008 = structure(c(0, 0, 0, 0, 
    0, 0), format.stata = "%8.0g"), ethfrac_FIXED = structure(c(NA_real_, 
    NA_real_, NA_real_, NA_real_, NA_real_, NA_real_), label = "eth. frac", format.stata = "%8.0g"), 
    communist = structure(c(0, 0, 0, 0, 0, 0), label = "Communist Leader", format.stata = "%8.0g"), 
    mil = structure(c(1, 1, 1, 1, 1, 1), format.stata = "%9.0g"), 
    cw = structure(c(1, 1, 1, 1, 1, 1), format.stata = "%9.0g"), 
    age = structure(c(52, 53, 54, 55, 56, 57), label = "Current age", format.stata = "%9.0g")), row.names = c(NA, 
-6L), class = c("tbl_df", "tbl", "data.frame"))

解决建议 #

• 匹配计数过程数据结构：Stata的stset t, id(leadid)是将每个leadid的多条观测视为计数过程（重复时间点）数据，对应R中需要构造区间型Surv对象。修改代码：

  data <- data %>% mutate(tstart = t - 1)
  surv_obj <- Surv(data$tstart, data$t, data$c_coup)
  

  这样才和Stata的时间点定义逻辑一致。

• 加入聚类稳健标准误：Stata的stset id(leadid)默认会对leadid组内的相关性做调整，R中需要在coxph里明确指定cluster(leadid)参数，确保估计逻辑匹配：

  m1 <- coxph(surv_obj ~ legislature + lgdp_1 + growth_1 + exportersoffuelsmainlyoil_EL2008 + ethfrac_FIXED + communist + mil + cw + age + cluster(leadid), 
              data = data, method = "breslow")
  

• 统一数据排序：Stata会自动按id和时间排序，R中需手动对数据按leadid和t排序，避免顺序影响模型计算：

  data <- data[order(data$leadid, data$t), ]
  

• 验证缺失值处理：检查ethfrac_FIXED的NA在两个软件中是否都做了列表删除，可在R中用data_clean <- na.omit(data)后重新建模，确认和Stata的分析数据集一致。

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