4
$\begingroup$

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 ***
$\endgroup$
6
  • $\begingroup$ This can be done in R with the following formula: survObject ~ covar + cluster(id) for the coxph function. $\endgroup$
    – Cliff AB
    Oct 19, 2015 at 1:09
  • $\begingroup$ That clusters the standard errors by the id, but it doesn't affect the coefficients. See gist.github.com/6430d4d65e69109c03de $\endgroup$
    – Andrew
    Oct 19, 2015 at 1:17
  • $\begingroup$ 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. $\endgroup$
    – Cliff AB
    Oct 19, 2015 at 1:30
  • $\begingroup$ 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) $\endgroup$
    – Andrew
    Oct 19, 2015 at 1:43
  • $\begingroup$ 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. $\endgroup$
    – Cliff AB
    Oct 19, 2015 at 1:52

1 Answer 1

8
$\begingroup$

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.

$\endgroup$
3
  • 5
    $\begingroup$ Welcome to CV! We hope you'll stay around & answer more questions from time to time. $\endgroup$ Oct 19, 2015 at 14:43
  • $\begingroup$ Perfect! This fixes it. I've included an example of this in the original question. $\endgroup$
    – Andrew
    Oct 19, 2015 at 15:20
  • $\begingroup$ Please consider registering your account. $\endgroup$ Oct 19, 2015 at 15:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.