How can I define a multiple-record ID variable for a Cox model in R?

Question

I'm reproducing and translating data from an old Stata file into R for future expansion and analysis. The original data uses a Cox proportional hazard model with stset + stcox. Running Cox models in R is trivial with survival::coxph. However, the models in Stata define a multiple-record ID variable for the time series, and R seems to not be able to do the same.

Here's some example data. Running a model in Stata without defining an id variable is trivial:

import delimited "https://gist.githubusercontent.com/andrewheiss/c60cb1c9e55aaf42585f/raw/9c5568d7af43c9cf523cf72db2591e80b466b69a/df.csv", clear
stset years_since_event fail
stcox x1 x2, nohr

This produces these results:

------------------------------------------------------------------------------
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |   .1963139   .0375578     5.23   0.000      .122702    .2699259
          x2 |  -.5599549   .0642704    -8.71   0.000    -.6859225   -.4339872
------------------------------------------------------------------------------

I can produce the same results in R with survival::coxph:

library(survival)
library(readr)
df <- read_csv("https://gist.githubusercontent.com/andrewheiss/c60cb1c9e55aaf42585f/raw/9c5568d7af43c9cf523cf72db2591e80b466b69a/df.csv", na=".")
model1 <- coxph(Surv(years_since_event, fail) ~ x1 + x2, data=df, ties="breslow")  # Breslow because Stata
summary(model1)

Results:

       coef exp(coef) se(coef)      z Pr(>|z|)    
x1  0.19631   1.21691  0.03756  5.227 1.72e-07 ***
x2 -0.55995   0.57123  0.06427 -8.712  < 2e-16 ***

In Stata, I can define a multiple-record ID variable for the time series by appending id(name) to stset, yielding results that are fairly substantially different from the first model, like so:

stset years_since_event fail, id(name)
stcox x1 x2, nohr

Results:

------------------------------------------------------------------------------
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |   .1584929   .0438302     3.62   0.000     .0725874    .2443984
          x2 |  -.3068401     .06471    -4.74   0.000    -.4336694   -.1800108
------------------------------------------------------------------------------

However, I have not found a way to define a multiple-record ID variable in R—there are no obvious arguments to adjust in either coxph() or Surv().

So, is there a way reproduce the stset years_since_event fail, id(name) Stata command in R using survival::coxph?

UPDATE: Following Terry Therneau's answer, the way to fix this is to explicitly add a variable indicating the start time. Stata does this behind the scenes when you include id(name), while R does not.

Here's what that looks like:

library(dplyr)
df <- df %>%
  group_by(name) %>%
  mutate(start_time = 0:(n() - 1))
head(df)

# Source: local data frame [6 x 6]
# Groups: name [3]
# 
#          name years_since_event  fail       x1    x2 start_time
#         (chr)             (int) (int)    (dbl) (int)      (int)
# 1 Afghanistan                 1     0 17.07171     6          0
# 2 Afghanistan                 2     1 17.10008     5          1
# 3     Albania                 1     1 14.93779     1          0
# 4     Algeria                 1     0 17.23435     2          0
# 5     Algeria                 2     0 17.24892     2          1
# 6     Algeria                 3     0 17.26365     2          2

model2 <- coxph(Surv(start_time, years_since_event, fail) ~ x1 + x2, 
                data=df, ties="breslow") 
summary(model2)

#        coef exp(coef) se(coef)      z Pr(>|z|)    
# x1  0.15849   1.17174  0.04383  3.616 0.000299 ***
# x2 -0.30684   0.73577  0.06471 -4.742 2.12e-06 ***

Answer 1

It looks like Stata has implemented a variant of the counting process style of survival model, which uses (time1, time2, status) in R. See chapter 3.7 of Therneau and Grambsch (Springer 2000) for an explanation of how this works. The start,stop form is much more flexible, and appears to be what the enter/exit variables do in Stata. If you use "id" then it builds the data set for you. In R you need to create it yourself.