Cox比例風險模型的Stata VS

我正嘗試使用以下數據http://iojournal.org/wp-content/uploads/2015/05/FortnaReplicationData.dta Cox比例風險模型的Stata VS

在Stata該命令中的R中複製的Stata從一個Cox比例風險模型估計如下：

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，I'm使用下列內容：

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結果之間爲什麼會有所不同？

來源

2015-11-27 user2246905

幾個意見（最有可能無益）。對於R結果，漸近和穩健的se很接近，我傾向於找到讓人放心的地方，並且z統計量可以看作是從coef/rob.se中計算出來的。我似乎無法計算stata結果中的z-stat（log（HR）/ rob.se不是） - 你知道爲什麼/如何？建議st.errors已被轉換maybies？ – user20650

我認爲在某種程度上，se可能會發生變化，但我真的不清楚它們是如何或是否真的轉化了。 – user2246905

林猜測狂放，但你有沒有嘗試指定'nohr'到你的靜態代碼.. – user20650

正如user20650注意到的，在Stata選項中包含「nohr」時，您將得到與R中完全相同的標準錯誤。使用羣集時仍然存在標準錯誤的小差異。 user20650再次注意到，由於Stata默認的標準錯誤與g /（g-1）相乘，因此給出了區別，其中g是羣集數量，而R不調整這些標準錯誤。因此，一個解決方案只是noadjust包括在Stata或者R中調整的標準誤差做：

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

如果我們仍然希望R中有從塔塔相同的標準誤差，因爲當不指定nohr，我們需要知道當nhr停止時，我們得到$ exp（\ beta）$與由這些比例擬合模型產生的標準誤差。特別是通過將delta方法應用於原始標準誤差估計而獲得。「德爾塔方法通過計算相應的一階泰勒展開的方差來獲得變換後的變量的標準誤差，對於變換$ exp（\ beta）$等於將標準誤差減去$ exp（\ hat { \ beta}）$。這種計算技巧會產生與在估計之前轉換參數相同的rsults，然後重新估計「（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

來源

2015-11-28 06:24:36 user2246905

非常感謝，對我非常有用。我有STATA的標準錯誤（0.7），HR（1.88），但是，由於我沒有這些數據，我如何使用R來獲得如R中的標準錯誤。羣集數爲182. – user2669497

i已經使用「（SE/HR）*（g-1/g）」來直接計算從STATA到SE的SE。以HCTrebels爲例，（0.1299916/0.4089758）*（48/49）= 0.31136，這在R中非常接近0.3146。 – user2669497

Cox比例風險模型的Stata VS

回答

相關問題