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

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)
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))

# 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 ***

• This can be done in R with the following formula: survObject ~ covar + cluster(id) for the coxph function. Oct 19, 2015 at 1:09
• That clusters the standard errors by the id, but it doesn't affect the coefficients. See gist.github.com/6430d4d65e69109c03de Oct 19, 2015 at 1:17
• Ah. I'm not terribly familiar with exactly what Stata is doing. I think there is a good chance that it is fitting a frailty model. This can be fit in R, see ?frailty in the survival package. Perhaps unfortunately, you have a few choices of the distribution of the frailty terms, so it might be worth checking exactly what Stata is doing. Oct 19, 2015 at 1:30
• It doesn't look like frailty or any sort of mixed effects—it looks like Stata is trying to do some sort of panel correction or something (stata.com/manuals13/ststset.pdf) Oct 19, 2015 at 1:43
• Where exactly is that shown? Stata certainly fits frailty models (see stata.com/manuals13/ststcox.pdf on page 30). But this particular manual does not seem to say what calls actually fit the frailty model. Oct 19, 2015 at 1:52