R和Stata中的Cox比例风险模型

Question

R和Stata中的Cox比例风险模型

5

我正在尝试使用R复制Stata中的Cox比例风险模型估计，使用以下数据： http://iojournal.org/wp-content/uploads/2015/05/FortnaReplicationData.dta 在Stata中使用的命令如下：

stset enddate2009, id(VPFid) fail(warends) origin(time startdate)
stcox HCTrebels o_rebstrength demdum independenceC transformC lnpop lngdppc africa diffreligion warage if keepobs==1, cluster(js_country)

Cox regression -- Breslow method for ties

No. of subjects      =          104                Number of obs   =       566
No. of failures      =           86
Time at risk         =       194190
                                               Wald chi2(10)   =     56.29
Log pseudolikelihood =   -261.94776                Prob > chi2     =    0.0000

                           (Std. Err. adjusted for 49 clusters in js_countryid)
-------------------------------------------------------------------------------
              |               Robust
           _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
    HCTrebels |   .4089758   .1299916    -2.81   0.005     .2193542    .7625165
o_rebstrength |   1.157554   .2267867     0.75   0.455     .7884508    1.699447
       demdum |   .5893352   .2353317    -1.32   0.185     .2694405    1.289027
independenceC |   .5348951   .1882826    -1.78   0.075      .268316    1.066328
   transformC |   .5277051   .1509665    -2.23   0.025     .3012164    .9244938
        lnpop |   .9374204   .0902072    -0.67   0.502     .7762899    1.131996
      lngdppc |   .9158258   .1727694    -0.47   0.641     .6327538    1.325534
       africa |   .5707749   .1671118    -1.92   0.055     .3215508    1.013165
 diffreligion |   1.537959   .4472004     1.48   0.139      .869834    2.719275
       warage |   .9632408   .0290124    -1.24   0.214     .9080233    1.021816
-------------------------------------------------------------------------------

我使用R语言完成以下操作:

data <- read.dta("FortnaReplicationData.dta")
data4 <- subset(data, keepobs==1)
data4$end_date <- data4$`_t`
data4$start_date <- data4$`_t0`
levels(data4$o_rebstrength) <- c(0:4)
data4$o_rebstrength <- as.numeric(levels(data4$o_rebstrength[data4$o_rebstrength])
data4 <- data4[,c("start_date", "end_date","HCTrebels",  "o_rebstrength", "demdum", "independenceC", "transformC", "lnpop", "lngdppc", "africa", "diffreligion", "warage", "js_countryid", "warends")]
data4 <- na.omit(data4)
surv <- coxph(Surv(start_date, end_date, warends) ~ HCTrebels+ o_rebstrength +demdum + independenceC+ transformC+ lnpop+ lngdppc+ africa +diffreligion+ warage+cluster(js_countryid), data = data4, robust = TRUE, method="breslow")

                 coef exp(coef) se(coef) robust se     z      p
HCTrebels     -0.8941    0.4090   0.3694    0.3146 -2.84 0.0045
o_rebstrength  0.1463    1.1576   0.2214    0.1939  0.75 0.4505
demdum        -0.5288    0.5893   0.4123    0.3952 -1.34 0.1809
independenceC -0.6257    0.5349   0.3328    0.3484 -1.80 0.0725
transformC    -0.6392    0.5277   0.3384    0.2831 -2.26 0.0240
lnpop         -0.0646    0.9374   0.1185    0.0952 -0.68 0.4974
lngdppc       -0.0879    0.9158   0.2060    0.1867 -0.47 0.6377
africa        -0.5608    0.5708   0.3024    0.2898 -1.94 0.0530
diffreligion   0.4305    1.5380   0.3345    0.2878  1.50 0.1347
warage        -0.0375    0.9632   0.0405    0.0298 -1.26 0.2090

Likelihood ratio test=30.1  on 10 df, p=0.000827
n= 566, number of events= 86

我得到相同的风险比系数，但标准误差看起来不一样。Z 值和 p 值接近但不完全相同。为什么 R 和 Stata 的结果会有差异？

- user2246905

一些注释（可能不太有用）。对于R的结果，渐近和鲁棒标准误非常接近，这让我感到放心，而z统计量可以从coef / rob.se计算得出。我似乎无法从stata的结果中计算出z统计量（log（HR）/ rob.se不是吗）-您知道为什么/如何吗？可能表明标准误已被转换？ - user20650

我认为在某种程度上，这些可能会发生变化，但我真的不清楚它们是如何变化的，或者它们是否真的发生了变化。 - user2246905

我猜测一下，你是否尝试在Stata代码中指定“nohr”？ - user20650

1

哈哈...搞定了!!! 找到一台旧笔记本电脑，用Stata添加“noadjust”。手册中有几句话。 - user20650

1

使用手册第3页中的调整公式... sqrt(diag(vcov(surv))* (49/48)) - 自动化聚类数量可能是值得的。 - user20650

显示剩余5条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user2246905 · Accepted Answer

正如用户20650所注意到的那样，当在Stata选项中包含“nohr”时，您会得到与R中完全相同的标准误差。但是，在使用群集时，标准误差仍然存在一些小差异。用户20650再次注意到，这种差异是因为Stata默认标准误差被乘以g /（g-1），其中g是群集数，而R不调整这些标准误差。因此，解决方案就是在Stata中包含noadjust或通过执行以下操作来调整R中的标准误差：

sqrt(diag(vcov(surv))* (49/48))

如果我们希望在R中获得与未指定nohr时相同的Stata标准误差，则需要知道，当省略nhr时，我们会获得$exp(\beta)$，其标准误差是由在该比例下拟合模型产生的。特别地，通过将delta方法应用于原始标准误差估计而获得。 "Delta方法通过计算相应一阶泰勒展开的方差来获得变换变量的标准误差，对于变换$exp(\beta)$而言，这相当于将原始标准误差乘以$exp(\hat{\beta})$。这种计算技巧可以获得与在估计之前转换参数然后重新估计相同的结果"（Cleves et al 2010）。在R中，我们可以使用以下代码：

library(msm)
se <-diag(vcov(surv)* (49/48))
sapply(se, function(x) deltamethod(~ exp(x1), coef(surv)[which(se==x)], x))

     HCTrebels o_rebstrength    demdum independenceC transformC     lnpop   lngdppc    africa diffreligion     warage
     0.1299916     0.2267867 0.2353317     0.1882826  0.1509665 0.0902072 0.1727694 0.1671118    0.4472004 0.02901243