How to Replicate Survival Analysis Conducted in Stata to R

Question

I am curious on how to replicate the results of a Cox model estimated in Stata in R. I get completely different results, and I think I understand why (but I am unsure of how to fix it), but before I discuss that, I will present a glimpse of the data I am working with along with the simplified syntax to illustrate my problem.

An idea of the data:

"Failure" denotes an instance of a country experiencing conflict recurrence (this is just dummy data here so its not accurate) and "peace_sp_id" is a unique numeric ID value for a spell of peace. So now, I will run a model in both Stata and R (with the real data), provide the syntax and model outputs:

Stata: stset year, id(peace_sp_id) failure(failure) stcox GDP

R: coxm1 <- coxph(Surv(year, failure) ~ GDP, data = chapter1data) coxm1

A big difference that I notice is that I'm not sure how to recreate the id(peace_sp_id) aspect in R. This is important since each row in this data set is not an independent observation from the other. Is there a way to re-create the "multiple-record ID variable" option in R that Stata provides?

Answer 1

(year might not be the correct time variable to stset on, but can't say for sure because your posted data snippet's dummy data. Let's bracket that issue, for now.)

The id() option for stset doesn't adjust the standard errors or coefficients in the way you're potentially thinking. id() only impacts how _t0 and _t are defined:

webuse diet2, clear

stset dox, id(id) failure(fail)
list id _st-_t0 in 1/5

stset dox, failure(fail)
list id _st-_t0 in 1/5

_t0 and _t are often referred to as "start" and "stop" variables, from counting-process notation.

In R, you need to generate the start and stop variables yourself. When you specify your Cox model, you'll then pass three arguments to Surv() instead of two:

coxm1 <- coxph(Surv(YOUR_NEW_START_VARIABLE_NAME_HERE, year, failure) ~ GDP, 
               data = chapter1data, ties="breslow") 

Note, too, the addition of ties="breslow" to coxph(), to match stcox's default tie correction.

With your dataset snippet (and continuing with the initial stset syntax you provided), I suspect chapter1data$year0 <- year - 1 will get you what you need. Verify that by checking _t0's values in Stata after running your current stset statement, though.

(Or, alternatively: save your dataset in Stata after you stset, import that dataset into R, then:

coxph(Surv(`_t0`, `_t`, `_d`) ~ GDP, data = chapter1data_stset, ties="breslow")

)