如何解决R与Stata中Weibull生存分析结果不一致问题？

原Stata代码 #

stset time, id(leadid) failure(c_coup)

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

我的R代码 #

# load packages
library(survival)
library(survminer)

# load data
data <- read.csv("leader_tvc_2.csv", encoding = "latin1")

# create survival object
surv_object <- 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

# run survreg
fit <- survreg(surv_object, data = data, dist = "weibull")

# examine output 
summary(fit)

结果差异 #

Stata中legislature的估计值为2.057，而R中对应值为2.48223，其他系数也存在类似差异。R输出结果如下：

Call:
survreg(formula = surv_object, data = data, dist = "weibull")
                           Value Std. Error     z       p
(Intercept)              3.69422    0.62813  5.88 4.1e-09
legislature              2.48223    0.21580 11.50 < 2e-16
leg_growth_2             4.06103    1.81609  2.24  0.0253
gdp_1k                  -0.02217    0.03678 -0.60  0.5467
chgdpen_fearonlaitin     0.49579    1.11792  0.44  0.6574
Oil_fearonlaitin         0.24194    0.23799  1.02  0.3094
postcoldwarlag           0.55797    0.31626  1.76  0.0777
civiliandictatorshiplag -1.87331    0.32142 -5.83 5.6e-09
militarydictatorshiplag -1.34963    0.29635 -4.55 5.3e-06
communist                0.88426    0.28211  3.13  0.0017
lpopl1_fearonlaitin     -0.00252    0.06636 -0.04  0.9697
ethfrac_fearonlaitin     0.47922    0.28625  1.67  0.0941
relfrac_fearonlaitin     0.67571    0.37836  1.79  0.0741
age                      0.01261    0.00704  1.79  0.0733
Log(scale)              -0.12430    0.06339 -1.96  0.0499

Scale= 0.883 

Weibull distribution
Loglik(model)= -829   Loglik(intercept only)= -978
    Chisq= 297.98 on 13 degrees of freedom, p= 6.4e-56 
Number of Newton-Raphson Iterations: 18 
n=3774 (5911 observations deleted due to missingness)

解决建议 #

1. 验证数据与变量一致性 #

• 检查所有自变量的编码、数值范围：比如在R中用summary(data$legislature)查看变量分布，对比Stata的sum legislature结果，确保0/1分类、连续变量取值完全一致。
• 确认time变量的定义：是否为统一的时间单位（如年/月），是否为从进入研究到事件发生或删失的时长。
• 核对事件变量c_coup：Stata中failure(c_coup)指定的事件触发值（如1=政变发生）需与R中Surv(time, c_coup)的规则一致（R默认非0值为事件）。

2. 处理个体聚类效应 #

Stata的stset id(leadid)标识了同一领导人（leadid）的多条观测，需在R中显式处理个体内相关性：

方法1：计算聚类标准误

library(sandwich)
library(lmtest)

# 拟合基础模型
fit <- survreg(surv_object, data = data, dist = "weibull")

# 基于leadid计算聚类稳健标准误
vcov_cluster <- vcovCL(fit, cluster = data$leadid)

# 输出调整后的结果
coeftest(fit, vcov = vcov_cluster)

方法2：拟合共享frailty模型（若Stata隐含使用）

如果原Stata模型实际包含个体脆弱性（可能遗漏了frailty(leadid)选项），则在R中拟合frailty模型：

fit_frailty <- survreg(surv_object + frailty(leadid), data = data, dist = "weibull")
summary(fit_frailty)

3. 参数化转换匹配Stata结果 #

Stata与R的Weibull AFT模型参数化存在缩放差异，需将R的系数转换为Stata格式：

• Stata系数 = R系数 / R输出的Scale值（此处为0.883）
• 示例：R中legislature系数2.48223转换后为2.48223 / 0.883 ≈ 2.81，若仍不匹配，需排查模型设定差异。

4. 核对样本与缺失值处理 #

• 确认Stata是否使用了if/in筛选样本，R中需对应添加数据子集筛选（如data <- subset(data, ...)）。
• 对比Stata与R的有效样本量：Stata的count结果需与R的nrow(na.omit(data))一致，确保缺失值处理规则相同。

5. 对比对数似然值 #

若R的Loglik(model)与Stata模型的对数似然值接近，说明仅为参数化或标准误差异；若差异较大，则存在数据或模型设定的本质区别。

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