在R中复现Stata生存分析时如何获取完全一致的标准误？

我正在复现某期刊用Stata生成的生存分析结果，目前在R中得到的系数（风险比HR）和显著性水平（除轻微舍入差异外）与原结果完全一致，但标准误不匹配。例如原结果中模型1（Natural Causes）的Legislature系数（0.456*）标准误为0.198，而我在R中得到的是0.414。推测差异源于标准误的转换方式（可能和delta方法有关），求解决建议。

原Stata结果说明 #

原结果为包含三个Cox回归模型的表格：Natural Causes、Coups、Revolts，输出各变量的风险比（HR）、标准误及显著性标记。

我的R代码 #

# 加载包
library(dplyr)
library(foreign)
library(msm)
library(stargazer)

# 导入原始数据
data <- read_stata("leaders, institutions, covariates, updated tvc.dta")

# 为每个leadid生成t0
data <- mutate(data, t0 = lag(t, default = 0), .by = leadid)

# 构建生存对象与模型：Coups
survobj_coup <- Surv(data[["t0"]], data[["_t"]], data$c_coup)
coups_original <- coxph(survobj_coup ~ legislature + lgdp_1 + growth_1 + 
                          exportersoffuelsmainlyoil_EL2008 + ethfrac_FIXED + 
                          communist + mil + cw + age, 
                        data = data, ties = "breslow")

# 构建生存对象与模型：Revolts
survobj_revolt <- Surv(data[["t0"]], data[["_t"]], data$c_revolt)
revolt_original <- coxph(survobj_revolt ~ legislature + lgdp_1 + growth_1 + 
                           exportersoffuelsmainlyoil_EL2008 + ethfrac_FIXED + 
                           mil + cw + age, 
                         data = data, ties = "breslow")

# 构建生存对象与模型：Natural Causes
survobj_natural <- Surv(data[["t0"]], data[["_t"]], data$c_natural)
natural_original <- coxph(survobj_natural ~ legislature + lgdp_1 + growth_1 + 
                            exportersoffuelsmainlyoil_EL2008 + ethfrac_FIXED + 
                            communist + mil + cw + age, 
                          data = data, ties = "breslow")

# 定义系数指数化函数
exp_coef <- function(x) {exp(x)}

# 使用stargazer生成表格
stargazer(natural_original, coups_original, revolt_original, 
          apply.coef = exp_coef, p.auto = FALSE)

问题核心与解决方案建议 #

1. 标准误的尺度转换问题 #

你当前用apply.coef = exp_coef仅将线性系数转换为HR，但stargazer默认输出的是原始线性系数的标准误，而非HR尺度的标准误。Stata输出的是HR对应的标准误，这需要用delta方法计算：

HR的标准误 = 原始线性系数的标准误 × HR值

比如原结果中Legislature的HR为0.456，标准误0.198，对应原始线性系数为log(0.456)≈-0.783，原始标准误应为0.198 / 0.456≈0.434，和你得到的0.414接近（舍入差异）。

2. 修改stargazer输出，手动传入HR尺度的标准误 #

不要依赖apply.coef，而是提前计算HR、HR的标准误和p值，手动传给stargazer：

# 整理所有模型
models <- list(natural_original, coups_original, revolt_original)

# 计算HR、HR的标准误、p值
hr_list <- lapply(models, function(x) exp(coef(x)))
se_hr_list <- lapply(models, function(x) sqrt(diag(vcov(x))) * exp(coef(x)))
p_list <- lapply(models, function(x) coef(summary(x))[, "Pr(>|z|)"])

# 生成匹配Stata格式的表格
stargazer(models, coef = hr_list, se = se_hr_list, p = p_list, 
          title = "Cox Regression Results (HR Scale)",
          column.labels = c("Natural Causes", "Coups", "Revolts"),
          type = "text")

3. 核对模型拟合细节 #

• 确认Stata与R的ties处理一致：你的R代码用了ties="breslow"，需核对原Stata代码是否用了stcox, breslow（Stata默认是efron）。不过你系数一致，大概率设置匹配，但仍需确认。
• 验证数据导入一致性：确保read_stata正确读取了Stata数据的变量类型、缺失值编码，与原分析使用的数据集完全一致。

4. 手动验证delta方法计算 #

运行summary(natural_original)提取Legislature的原始系数和标准误，计算原始SE × exp(原始系数)，看结果是否接近原Stata的0.198，以此确认转换逻辑正确。

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