Risk Ratios from Survival Data

[Sean Hardiman]

Hi all,

Does anyone have any experience calculating risk ratios from survival data generated sts list? I am using the following code to generate a table that gives me probability of death at 36 months.

sts list, failure risktable(1 12 24 36) by(study_grp)

What I'd like to do next is calculate risk ratios comparing each of the three of the study groups to the one reference group. This code gives me the outcome probability at the end of my study period in each of the relevant groups. What I wonder is whether there is a way to do this automatically or whether I need to manually take the results and do the calculation? I'm running Stata 17.

Thanks for any suggestions you can provide!

Sean

[Paul Dickman]

Is it correct that by "risk ratio" you mean the ratio of the outcome probabilities?

Do you want standard errors or confident intervals?

How many times do you need to do this?

If you don't want measures of uncertainty and you are doing it once then manually is probably easiest.

If you don't want measures of uncertainty, but want something "automatic" then you can use the -saving- option and then do the calculations.

If you want measures of uncertainty then you could fit a model and take the ratios of the predicted survival probabilities. I would fit a flexible parametric model (stpm2 from SSC).

[Sean Hardiman]

Hi Paul,

Thank you for your reply! Answers to your questions follow.

Is it correct that by "risk ratio" you mean the ratio of the outcome probabilities?

Yes, the ratio of outcome probabilities is what I want to estimate.

Do you want standard errors or confident intervals?

Confidence intervals.

How many times do you need to do this?

I have a primary analysis and three sub-group analyses, so four times in total.

If you don't want measures of uncertainty and you are doing it once then manually is probably easiest.

If you don't want measures of uncertainty, but want something "automatic" then you can use the -saving- option and then do the calculations.

If you want measures of uncertainty then you could fit a model and take the ratios of the predicted survival probabilities. I would fit a flexible parametric model (stpm2 from SSC).

I am unfamiliar with stpm2, so will give it a look! Am I correct that you're its author? If so, any tips you have for this particular application are welcomed!

Thanks again!

[Paul Dickman]

No, I am not the author of stpm2. Here's a simple example that illustrates the basic approach.

clear
use rott2
stset rf, f(rfi==1) scale(12) exit(time 120)


sts gen S_km = s, by(hormon)

// Fit model
stpm2 hormon, tvc(hormon) dftvc(2) scale(hazard) df(3) eform
predict s,s

// Compare model-based estimates to Kaplan-Meier estimates
twoway    (line S_km _t if hormon == 0, sort lcolor(black) lpattern(dash) connect(stepstair)) ///
        (line S_km _t if hormon == 1, sort lcolor(red)  lpattern(dash) connect(stepstair)) ///
        (line s _t if hormon==0,sort lcolor(black) lwidth(thick)) ///
        (line s _t if hormon==1, sort lcolor(red) lwidth(thick)) ///
        , xtitle("Years from surgery") ///
        ytitle("S(t)") ///
        legend(order(3 "No hormonal therapy" 4 "hormonal therapy") ring(0) pos(1) cols(1)) ///
        caption("Dashed lines show KM estimates")
        
// relative risk at time 5
generate t5=5

predictnl relrisk = predict(survival at(hormon 1) timevar(t5)) / ///
                predict(survival at(hormon 0) timevar(t5)) ///
                in 1, ci(rr_lci rr_uci)
                
list relrisk rr_lci rr_uci in 1

[Sean Hardiman]

Hi Paul,

Thanks for the sample code, this worked really well! Here's how I translated it:

sts generate failure_km = f, by(study_group)

// Fit model
stpm2 study_group, tvc(study_group) dftvc(2) scale(hazard) df(3) eform
predict failure, s

// Compare model-based estimates to Kaplan-Meier estimates
twoway  (line failure_km _t if study_group == 1, sort lcolor(black) lpattern(dash) connect(stepstair)) ///
        (line failure_km _t study_group == 3, sort lcolor(red)  lpattern(dash) connect(stepstair)) ///
        (line failure_km _t if study_group == 4, sort lcolor(blue)  lpattern(dash) connect(stepstair)) ///
        (line failure_km _t ifstudy_group == 5, sort lcolor(green)  lpattern(dash) connect(stepstair)) ///
        (line failure _t if study_group == 1, sort lcolor(black) lwidth(thick)) ///
        (line failure _t if study_group == 3, sort lcolor(red) lwidth(thick)) ///
        (line failure _t if study_group == 4, sort lcolor(blue) lwidth(thick)) ///
        (line failure _t if study_group == 5, sort lcolor(green) lwidth(thick)) ///
        , xtitle("Months from Start") ///
        ytitle("Proportion Surviving") ///
        legend(order(1 "Study Group 1" 3 "Study Group 2" 4 "Study Group 3" 5 "Study Group 4") ring(0) pos(1) cols(1)) ///
        caption("Dashed lines show KM estimates")

The good news is that the lines in the two-way plot that use failure_km are plotting correctly as failure rates. However, the other lines that use 'failure' instead show up as survival. I think that's because the 'predict failure, s' line tells Stata to calculate them as survival. Looking at the documentation for predict at https://www.stata.com/manuals/ststregpostestimation.pdf, there is no 'failure' option for predict that I can see. Do you know if there's a way for me to generate the complement of the survival curves? Otherwise, stpm2 is working well!

Edit: I've also noticed that the two-way plot isn't labelling the curves correctly. I should think I would have 8 entries in the legend, but I only have four, and one is mislabelled. I'm thinking I've not labelled the groups correctly in the legend line...will continue to work on this!

Thanks!

[Bjarte Aagnes]

for predict after stpm2 see

help stpm2 postestimation

[Paul Dickman]

Yes, as Bjarte said.

The syntax is

predict newvar, statistic

failure is one of the statistics you can predict.

There are 8 lines, but you have specified

 
  legend(order(1 "Study Group 1" 3 "Study Group 2" 4 "Study Group 3" 5 "Study Group 4")

Only lines 1, 3, 4, and 5 will be shown in the plot and they will be labelled as you have instructed.

[Sean Hardiman]

Hi everyone,

Thanks for the help on this. I'm back to it after some other work took precedence.

I've been able to implement these suggestions with success, so that's the good news. However, when I run the following code:

predictnl relrisk = predict(survival at(study_grp 3) timevar(t36)) / ///
                predict(survival at(study_grp 1) timevar(t36)) ///
                in 1, ci(rr_lci rr_uci)
                
list relrisk rr_lci rr_uci in 1

I get the following error:

(34,803 missing values generated)
warning: prediction constant over observations; perhaps you meant to run nlcom.
note: confidence intervals calculated using Z critical values.

It looks like only one row in my data table has the relative risk values. I'm not sure if I should run nlcom or not? Any advice on whether or not this is a problem? Still learning, so appreciate your help.

[Paul Dickman]

That is as it should be. By default, predictnl performs the calculation for every observation in the data set. The risk ratio will be constant, so we use the "in 1" qualifier to only perform the calculation for the first observation (to save time).

[Sean Hardiman]

Thank you, Paul! Makes sense to me now, very much appreciate your help!

[Sean Hardiman]

Is there a way to adjust the risk ratio results by inverse probability weights using this method? I want to estimate unadjusted (done) and adjusted risk ratios and have estimated weights to apply. Am looking for guidance on this in the documentation and not finding it. Any references I can read through? Thanks!

[Paul Dickman]

Is there a way to adjust the risk ratio results by inverse probability weights using this method? I want to estimate unadjusted (done) and adjusted risk ratios and have estimated weights to apply. Am looking for guidance on this in the documentation and not finding it. Any references I can read through? Thanks!

Where do you want to apply the weights?

-stpm2- supports fweights, iweights, and pweights (specified using stset). -standsurv- (available from SSC) allows individual weights to be applied when estimating marginal probabilities (i.e., the risks), but off the top of my hand it's not obvious how to get the ratios (but it may be obvious from the documentation).

[Sean Hardiman]

Hi Paul,

When you ask 'where' I want to apply the weight, I'm not certain. I've used stpm2 to fit the model here, so I'm thinking I'll use the weights in stset. You've given me a few things to check out, so I'll do the reading and report back. Thanks!

[Sean Hardiman]

Hi Paul,

Just following up - I've been able to make the weighted analyses work properly, so that's great.

If you don't want measures of uncertainty, but want something "automatic" then you can use the -saving- option and then do the calculations.

Could you guide me on how to do this? My committee would like to see estimates from both the KM curve and the STPM model, so I'm hoping to accommodate that request.

Thanks!

Sean

[Paul Dickman]

Hi Paul,

Just following up - I've been able to make the weighted analyses work properly, so that's great.

"If you don't want measures of uncertainty, but want something "automatic" then you can use the -saving- option and then do the calculations."

Could you guide me on how to do this? My committee would like to see estimates from both the KM curve and the STPM model, so I'm hoping to accommodate that request.

Thanks!

Sean

I was just thinking of something like this:

webuse stan3
sts list, by(posttran) risktable(1400) failure saving(foo, replace)
use foo, clear
generate ratio=failure/failure[1]
list failure ratio