Stratified Cox regression R vs. Stata - can R do more?

[Fabian Fortner]

Hi
I am a regular Stata user, but I'll also looking into solving problems in R. Here, I have problem I can't solve in Stata, but in R. Does anyone know how to do this in Stata?

Those are equal commands:

   stcox sex_rf i.bsi, strata(timeband)

fit <- coxph(Surv(t0,t, event==1) ~ sex_rf
                        +bsi_active +strata(timeband), data=data )

However, according to this article (subchapter "Step function to explore time-varying coefficient"), I would like to fit something accoring to (note the: "+bsi_active:strata(timeband)")

fit <- coxph(Surv(t0,t, event==1) ~ sex_rf
                        +bsi_active:strata(timeband), data=data )

Is there a way to do this in Stata as well?

[Carlo Lazzaro]

Fabian:
have you already taken a look at the community-contributed module -scurve_tvc-?

[Paul Dickman]

I am not extremely familiar with R syntax, but it appears that's a Cox model stratified on timeband and with a bsi by timeband interaction. That is, I think the corresponding model can be estimated in Stata using

1.

stcox sex_rf i.bsi i.bsi#i.timeband, strata(timeband)

This is slightly off-topic, but I reflected over the difference between that model and this one:

2.

stcox sex_rf i.bsi i.bsi#i.timeband

Note that timeband (if OP has followed the same procedure as the cited R code) was created using -stsplit-. A stratified Cox model, when stratified on a covariate that is fixed at baseline, conceptually has multiple baseline hazards. That is, if we stratify on sex then we conceptually have a separate baseline for each sex. Model 1, however, only has one baseline hazard for each value of time and the effect of strata(timeband) would seem to be "make the baseline hazard discontinuous at each of the cutpoints of timeband". Since the baseline is not actually a continuous function that is estimated, it's not obvious to me what the effect of "forcing discontinuities" would be (in theory and in practice). I realise model 1 makes more sense when estimating the time varying effects of bsi (which are by definition discontinuous) but I have never given much thought to the differences between these two models.

[Fabian Fortner]

Hi
Thank you very much for your replys!

@carlo:
I read the documentation but I'm not sure it would help in this case.

@Paul:
I tried your suggestions:

 stcox sex_rf i.bsi_active#i.timeband, strata(timeband)

-------------------------------------------------------------------------------------
                 _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
             sex_rf |   1.031496   .0898074     0.36   0.722     .8696758    1.223425
                    |
bsi_active#timeband |
          0#28-72d  |   .1163187   .0462414    -5.41   0.000     .0533657    .2535343
         0#72-365d  |   2.043443   .2970399     4.92   0.000     1.536844    2.717035
           0#365+d  |   .0262389   .0048015   -19.89   0.000     .0183309    .0375584
           1#0-28d  |   6.787313   4.204989     3.09   0.002     2.015335    22.85854
          1#28-72d  |   .5785286          .        .       .            .           .
         1#72-365d  |   38.82113          .        .       .            .           .
           1#365+d  |          1          .        .       .            .           .
-------------------------------------------------------------------------------------
                                                        Stratified by timeband

However,

 
 fit <- coxph(Surv(t0,t, event==1) ~ sex_rf    +bsi_active:strata(timeband), data=data )

Results in this, whereas tgroups are equal to the 4 timebands in the stata output:

--> Some coefficients are equal, others not at all and Stata seems not to estimate the standard errors for some coefficients.

[Paul Dickman]

I believe we are on the right track, but we have to get Stata to use the same paramerisation as R. I think the following might work.

stcox sex_rf i.bsi_active#i.timeband i.timeband, strata(timeband) efron

Sorry I got it wrong the first time; it's tough without data so you might need to experiment a bit.

Stata seems not to estimate the standard errors for some coefficients

This is to be expected, because you're modelling time in two ways (as the underlying timescale and explicitly as a step function with the variable timescale). That is, there is collinearity. I think you can expect the effects of i.timeband to look strange, but hopefully the estimates for i.bsi_active#i.timeband (interpreted as the HRs for BSI within each timeband) will look OK.

I added the efron option so that the Stata results will (hopefully) mirror the R results. R, by default, uses the Efron method for ties and Stata uses the Breslow approximation.

I also recommend

efron approximation for

[Fabian Fortner]

I'm so sorry, I thought I already replied to your post Paul!
I'm very grateful, as this worked like a charm! thanks so much!