在R中复现Stata生存分析结果出现偏差的原因排查求助

问题背景 #

我正尝试复现某期刊发表的生存分析结果，Stata运行原代码后结果与期刊完全一致，但R中复现结果偏差较大，求排查原因。

原Stata代码 #

* COUPS
gen c_coup=c
replace c_coup=0 if exit!="coup"
stset time, id(leadid) failure(c_coup)

* REVOLT ENTRY LEFT OUT BECAUSE IT IS A PERFECT PREDICTOR
streg  legislature  leg_growth_2 gdp_1k chgdpen_fearonlaitin Oil_fearonlaitin postcoldwarlag civiliandictatorshiplag militarydictatorshiplag communist lpopl1_fearonlaitin ethfrac_fearonlaitin relfrac_fearonlaitin  age, distribution(weibull) time
outreg2 using survival, replace ctitle(coups, partial) tex nonotes bdec(3) e(all)
stcurve, hazard

我的R代码 #

# load required libraries
library(survival)
library(haven)

# load data 
leader_tvc_2 <- read_dta("leader_tvc_2.dta", encoding = "latin1")

# create survival object 
surv_obj_coup <-
  Surv(time = leader_tvc_2$time, event = leader_tvc_2$c_coup)

# fit a survival regression model in R
surv_model <- survreg(
  surv_obj_coup ~ legislature + leg_growth_2 + gdp_1k + chgdpen_fearonlaitin + Oil_fearonlaitin + postcoldwarlag + civiliandictatorshiplag + militarydictatorshiplag + communist + lpopl1_fearonlaitin + ethfrac_fearonlaitin + relfrac_fearonlaitin + age,
  data = leader_tvc_2,
  dist = "weibull"
)

# summarize results
summary(surv_model)

核心差异点及修正方案 #

1. 模型参数化的尺度差异（time选项） #

Stata代码中streg命令的time选项是关键：它将模型的因变量从对数时间（log(time)）改为原始时间（time），直接拟合线性形式的加速失效时间（AFT）模型：

T = Xβ + σε

而R的survreg默认拟合的是对数时间尺度的AFT模型：

log(T) = Xβ + σε

两种参数化方式的系数量级和解释完全不同，这是结果偏差的核心原因之一。

2. 个体聚类的处理差异（id(leadid)） #

Stata的stset, id(leadid)指定leadid为个体标识符，用于识别同一个体的多条记录（比如时变协变量场景），会自动调整基线风险和标准误。而你的R代码未处理聚类，将每条记录视为独立个体，导致估计偏差。

3. 修正后的R代码 #

推荐使用flexsurv包实现更灵活的参数化，匹配Stata的设置：

library(survival)
library(haven)
library(dplyr)
library(flexsurv)

# 加载数据
leader_tvc_2 <- read_dta("leader_tvc_2.dta", encoding = "latin1")

# 严格匹配Stata的c_coup生成逻辑
leader_tvc_2 <- leader_tvc_2 %>%
  mutate(c_coup = case_when(exit == "coup" ~ c, TRUE ~ 0))

# 拟合匹配Stata的Weibull模型：原始时间尺度 + 个体聚类调整
surv_model <- flexsurvreg(
  formula = Surv(time, c_coup) ~ legislature + leg_growth_2 + gdp_1k + chgdpen_fearonlaitin + 
    Oil_fearonlaitin + postcoldwarlag + civiliandictatorshiplag + militarydictatorshiplag + 
    communist + lpopl1_fearonlaitin + ethfrac_fearonlaitin + relfrac_fearonlaitin + age +
    cluster(leadid),  # 匹配Stata的id(leadid)
  data = leader_tvc_2,
  dist = "weibull",
  scale = "identity"  # 匹配Stata的time选项，用原始时间尺度
)

# 查看结果
summary(surv_model)

4. 额外验证项 #

• 确认c_coup变量在R和Stata中完全一致：检查是否有缺失值、编码错误
• 验证数据格式：如果是时变协变量的长格式数据，需改用Surv(start, stop, event, type = "counting")构建生存对象
• 检查Stata模型的输出类型：若实际为比例风险（PH）形式的Weibull模型，R中需用coxph替代survreg/flexsurvreg

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