Calculating and plotting graph on relative risk of death over time

[Napat Lee]

Hello all,

I am trying to analyze patient survival (outcome variable: DEATH) of patients who received kidney transplant compared to patients who remained on dialysis (treatment variable: TRANSPLANT).

KM curve is below. As you can see that transplanted patients appear to have worse survival initially. The benefit of transplant started to show around more than a year from transplant.



Next, what I would like to do is to look at relative risk of death at time intervals (0-3 months, 3 months-1year, 1-2 years, 2-3 years, 3-5 years, and more than 5 years).

I would like to get a graph that looks similar to this (see below, I cropped this from a paper I find online).



I don't really know who to perform relative risk of death. I spent a lot of time on here trying to figure it out. What I have done so far was that I did stssplit like this (see below).



Am I doing it the right way? So if this is the right way, does it mean that the benefit of transplant only started to show after 3 years post-transplant? And how do I plot a graph that look like the graph above?

Thank you so much!

NL

[Clyde Schechter]

Well, I'm not sure exactly what you want, as the language gets a bit confusing. People often use the terms hazard ratio and relative risk as if they are the same thing, but they are not. And I can't tell how you are using the terms.

Strictly speaking, the relative risk of death at, say, 730 days would be calculated by looking at all the people who enter the cohort at time 0, determining the proportion of those who died up to time 730 days in each arm (the risk), and calculating the ratio of those two risks. Similarly for the other time points.

The Cox proportional hazards analysis does not give you relative risks. It gives you hazard ratios, which are different. The hazard ratio at 730 days looks at all the people who survived up to day 730, then calculates the proportion of those in each arm who died exactly at time 730, or, for practical purposes, within a very short time after time 730 (the hazard), and calculating the ratio of those two hazards.

In general the numbers are different, often very different. Notice that the numerator and the denominator of the hazard are both different from the numerator and denominator of the risk. So the ways of calculating them are different, as would be the results.

So you need to clarify which you want: relative risk or hazard ratio.

[Napat Lee]

Apologies for the confusion. I guess I did not know what I wanted, and thank you for pointing that out. Papers that I have found seem to use relative risk and hazard ratios interchangeably, which make me even more confused.

I think what I am looking for is hazard ration. Here is another sample for a different paper that I found. They compared patient survival and then looked at adjusted HR over time.



In the method section, it says, "Stratified Cox regression model was used to compare hazards by the following time periods post- transplant: 0 to <3 months, 3 to <12 months, 1 to <5 years, and 5 years."

So basically, what I am asking for is how to perform stratified cox regression and how to plot the graph. I am fairly new to biostatistics, so thank you so much for your patience.

[Clyde Schechter]

In the method section, it says, "Stratified Cox regression model was used to compare hazards by the following time periods post- transplant: 0 to <3 months, 3 to <12 months, 1 to <5 years, and 5 years."

Still more confusing language! In a true stratified Cox regression model it is mathematically impossible to compare the group hazards at different time periods. That's because in a stratified Cox regression model, the differences between the group is all attributed to different baseline hazard functions and the hazards are assumed to be the same in the groups. I think what they meant was that they made the group variable into a time-varying covariate, which is precisely what you did in what you show in #1. It is also, I would say, the right way to approach this.

So what you need now is some code to get you from that Cox regression output to the time-specific group hazard ratios, and to make a graph of them.

frame create results int time float hr

foreach t of numlist 0 90 365 730 1095 1825 {
    lincom 1.TRANSPLANT#`t'.timeband - 0.TRANSPLANT#`t'.timeband, eform
    frame post results (`t') (`r(estimate)')
}

frame results: graph twoway line hr line, sort

Note: As no example data was provided, this code is untested and may contain errors.

[Napat Lee]

OMG. I worked! You are a life saver! I am so grateful for your help.

[Napat Lee]

Sir,

I am not sure how to interpret the results.

frame create results int time float hr


foreach t of numlist 0 90 180 365 730 1095 1825 2555 {
    lincom 1.TRANSPLANT#`t'.timeband - 0.TRANSPLANT#`t'.timeband, eform
    frame post results (`t') (`r(estimate)')
}


 ( 1)  - 0b.TRANSPLANT#0b.timeband + 1.TRANSPLANT#0b.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   1.858652   .4801445     2.40   0.016     1.120231    3.083816
------------------------------------------------------------------------------

 ( 1)  - 0b.TRANSPLANT#90.timeband + 1.TRANSPLANT#90.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   1.254299   .4469706     0.64   0.525     .6238423    2.521898
------------------------------------------------------------------------------

 ( 1)  - 0b.TRANSPLANT#180.timeband + 1.TRANSPLANT#180.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   .7504125    .176492    -1.22   0.222     .4732637    1.189863
------------------------------------------------------------------------------

 ( 1)  - 0b.TRANSPLANT#365.timeband + 1.TRANSPLANT#365.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   .4833205    .090499    -3.88   0.000     .3348517    .6976183
------------------------------------------------------------------------------

 ( 1)  - 0b.TRANSPLANT#730.timeband + 1.TRANSPLANT#730.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   .7178118   .1495201    -1.59   0.111     .4772054    1.079732
------------------------------------------------------------------------------

 ( 1)  - 0b.TRANSPLANT#1095.timeband + 1.TRANSPLANT#1095.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   .7547792   .1495668    -1.42   0.156      .511854    1.112996
------------------------------------------------------------------------------

 ( 1)  - 0b.TRANSPLANT#1825.timeband + 1.TRANSPLANT#1825.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   .5921839   .1834338    -1.69   0.091     .3226934    1.086734
------------------------------------------------------------------------------

 ( 1)  - 0b.TRANSPLANT#2555.timeband + 1.TRANSPLANT#2555.timeband = 0

------------------------------------------------------------------------------
          _t |     exp(b)   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   .4218264   .2526088    -1.44   0.149     .1304355    1.364181
------------------------------------------------------------------------------

frame results: graph twoway line hr time, sort

Does this mean hazard ratio = 1.858652 at time 0? Or does it mean hazard ratio is 1.858652 from the period of 0 to 3 months?

Then, hazard ratio =n1.254299 at 3-month mark, or it means hazard ratio of 1.254299 from the period of 3 to 6 months?

[Clyde Schechter]

It depends on how the variable timeband was defined. If timeband itself corresponds to a single point in time, then the same is true of the hazard ratios you are working with here. If timeband is a variable that defines an interval , then the hazard ratios also refer to that interval.

[Napat Lee]

stset PATIENT_FU_TIME_NEW, failure(DEATH) id(WL_ID_CODE)

stsplit timeband, at(0 90 180 365 730 1095 1825 2555) 


stcox TRANSPLANT#i.timeband, strata(timeband)

I used stsplit to create timeband. So this means what I've got are time intervals and not single points, correct?

[Clyde Schechter]

Yes, the variable -stsplit- creates defines intervals.

[Nazli Uludere Aragon]

Dear Statalisters, and Clyde and Napat,

I realize the above thread is not very recent, but just wanted to follow up with a brief question.
If I understand correctly Clyde confirmed in #4 that Napat's approach in #1 to run:

stcox TRANSPLANT#i.timeband, strata(timeband)

where timeband is the variable created after stsplit is valid (for what Napat wanted to accomplish). I am not following why it was necessary to specify timeband as strata?

Separately, is it acceptable to run:

stcox TRANSPLANT if timeband==0

and

stcox TRANSPLANT if timeband!=0

if we were interested in comparing survival outcomes for the first 90 days, and then, conditional on surviving those first 90 days, for the remaining group? I hope my question is meaningful. Thank you.

[Clyde Schechter]

In #1, O.P. said

Next, what I would like to do is to look at relative risk of death at time intervals (0-3 months, 3 months-1year, 1-2 years, 2-3 years, 3-5 years, and more than 5 years).

and we subsequently clarified that he really meant hazard ratio, not relative risk.

So it is the very statement of the problem that calls for using the interaction term. It allows you to estimate a hazard ratio for the TRANSPLANT effect separately in each of those time intervals.

In non-survival (time-to-event) analysis, this kind of stratification can also be accomplished by simply carrying out the regression separately in each stratum. That yields different standard errors from the interaction approach, but the estimated hazard ratios will be the same either way. But in survival analysis you cannot do this. The reason is that if you restrict the analysis to any of the intervals other than the initial one (0-3 months in Napat's case), then the fitted survival functions will start from 1, and they will not reflect either the true survival up to the start of the interval, nor the fact that the two TRANSPLANT groups' survival functions are already separated from each other.

[Nazli Uludere Aragon]

Dear Clyde,

Thanks for the response. The time to event vs. survival distinction is still not very clear to me.

In my analysis failure is not death, it is some undesirable outcome, and treatment is a policy intervention. A large proportion of the treated sample "fails" early on (year 1) after the end of multi-year policy "treatment". Then the remaining, smaller portion does not experience the failure event (or survives) for longer after receiving the same treatment. The treated K-M curve looks like a hockey stick; steep early on and rather flat in latter periods.

With respect to the control (untreated), it looks like everyone is doing worse with treatment. (Let us assume the control is well-matched. Their analysis time is longer because they don't spend time in treatment).

With the split at year 1, the second group (whose members by construct have not failed in year 1), the policy intervention seems to prolong survival. Am I thinking correctly?

I hypothesize that despite being selected for the same treatment, these two groups are different. I also hypothesize that the same explanatory variables --e.g., climate (30-yr average annual temperature, "temp")-- influence outcomes differently for the group that fails early on versus the group that survives (time varying "effects").

So I planned to use stsplit

stsplit split, at(1)

and running a regression similar to O.P.'s:

stcox treatment treatment#i.split temp c.temp#i.split, strata(split)

I also considered carrying out the regression separately, as I explained in #10, for the two groups. Am I understanding your reply correctly that if I run separate regressions:

stcox treatment temp, if split==1

and

stcox treatment temp, if split!=1

is also ok?

Thanks again for your time, help, and insights!

[Clyde Schechter]

Survival analysis is the same thing as time to event analysis. It was first invented and used in mechanical engineering as a way of studying the durability of parts in mechanical systems. It was usually called failure analysis. When it got adopted by medicine, it was first used in mortality studies and so got the name of survival analysis. But it was later broadened to look at incidence rates as well. And the more general term that covers all of these applications is time-to-event analysis.

In fact, one could be more general than that. The statistical techniques are applicable to any outcome variable whose range is contained in the positive real numbers. And occasionally one will see economic or financial data analyzed in this way. But this is unusual, and calling it "Positive real number outcome analysis" probably would never catch on.

As for your question about doing the two separate -stcox- regressions, I think you misunderstood my final paragraph in #11. It is not OK to do that. You have to use the interaction method to get things right in this situation: this is one of the ways in which time-to-event analysis differs from other kinds of regression analysis.

[Nazli Uludere Aragon]

Thank you! I get it now: You are correct, I absolutely misunderstood #11 last paragraph.

If you have any thoughts on graphing and interpreting K-M curves after stsplit as I've shown, please let me know.