Generating predicted hazard ratio's after Cox regression models

[Ahmed Ibrahim]

Greetings good people,

I am trying to generate a figure with predicted hazard ratio's with 95%Ci's after running a Cox model.

Here is the sample data

* Example generated by -dataex-. For more info, type help dataex
clear
input float(predictor followup outcome)
.08   64 0
.08 2130 0
.07  590 0
.07  736 1
.09   43 1
.09 1894 0
 .1 1845 0
 .1    4 0
.12  452 1
.11 1810 0
.12  957 0
.11 1709 0
.11 1496 0
.13  814 0
.12 1837 0
.13 1933 0
.12 1260 0
.09 1771 0
.12 2025 0
.14    5 0
.11 1016 0
.17 1914 0
.18 1712 0
.18 1645 0
.12  533 0
.15 1826 0
.17 1829 0
.15 1403 0
 .2 2227 0
.11 1864 0
.21 1721 0
.21   80 0
.21  554 0
.17 1964 0
.22  245 0
.21 1785 0
 .2 1846 0
.22 2014 0
.15   45 0
.19  197 0
.24 1935 0
.18 1770 0
.23   53 1
.18 1865 0
.25    1 0
 .2    7 0
.25 2240 0
.22  903 1
.25  188 0
.28  147 1
 .3 1748 0
.28 1791 0
 .3 1578 0
.29 1288 1
.22    1 0
.24 2079 0
.32 1727 0
.34 1695 0
.28    1 0
.35  315 1
 .3  972 1
.26    9 0
.22  726 0
.24  177 0
 .3 1112 1
.33   16 0
.15    1 0
.37    1 0
.25 1320 1
.12 1797 0
.38  400 1
.38  624 1
.37 1525 0
.38 1502 0
.23 1561 1
.29    3 0
.31 2216 0
.36 2075 1
.42  140 0
.27 1401 0
.32    1 0
.08 1651 0
 .2 1797 0
.18 1911 0
.28  221 0
.47 1377 1
.44    3 1
.39  518 0
.27 1837 0
.45  768 0
.39 1734 0
.11    8 0
.27   63 0
.31 2065 0
 .5 2072 0
.51 2140 0
.14  564 1
.51 2138 0
.25  473 0
.38  293 1
end
label values outcome yn
label def yn 0 "No", modify
label def yn 1 "Yes", modify

When i run the code below - what i am getting seems to be predicted risk rather than predicted hazard ratio. Please advise.

            *customs check
                stset followup, failure(outcome) show scale(30.4)                
                            
            *run cox regression
                 stcox predictor
                 
            *predict hazard 
                 predict hr_predictor, hr

            *generate 95% ci
                 predict seyhat_predictor, stdp
                 gen lowerb_predictor = hr_predictor - invt(1205, 0.975)*seyhat_predictor
                 gen upperb_predictor = hr_predictor + invt(1205, 0.975)*seyhat_predictor
                 
            *sort
                 sort predictor

            *generate figure
                 twoway  (rarea lowerb_predictor upperb_predictor predictor, ///
                                 fcolor(red) fintensity(inten10) ///
                                 lcolor(black) lwidth(medium) lpattern(dash))  || ///
                          (line hr_predictor predictor, lcolor(black) lwidth(thick)) || ///
                          (scatteri 1 0 1 1, recast(line) ///
                                    lcolor(red) lwidth(thin) ///
                                    lpattern(dash) lpattern(dash)), ///
                            legend(off)

[Clyde Schechter]

Why do you think you are getting predicted risk? The labels on the vertical axis of the graph range from 0 to 8, and all of the plotted values are above 1. The predicted risk would necessarily be between 0 and 1. Moreover, a few back of the envelope calculations suggest that the plotted graph does agree with the hazard ratios predicted by your model and data. I think this is all working correctly. I don't see the problem.

[Ahmed Ibrahim]

Hey Clyde, belated happy new year

Thank you very much for the clarification. You are right, figure is correct.

The next issue I have is trying to plot hazard ratio's adjusted to a confounder.

Here is a new batch of the sample data

* Example generated by -dataex-. For more info, type help dataex
clear
input float predictor byte confounder float(followup outcome)
.39 81 1734 0
.84 69 2214 0
.91 74 2079 0
.89 59 2061 0
.71 54 1991 0
.76 70 2135 0
.08 66 2130 0
.88 73  865 1
.71 59 1728 0
.29 87 2021 0
.71 69 1652 0
.81 65 1717 0
.82 64  530 0
.87 85  903 1
.41 54 2131 0
.94 62 1752 0
.75 61 1808 0
 .4 84 1632 0
.94 56 2222 0
.83 59 1917 0
.84 71 1458 0
.33 67 1525 0
.29 86 1288 1
.81 61 2215 0
.73 55 1712 0
.42 52  140 0
.86 69 1862 0
.72 63   21 0
.95 81 1885 0
.25 68 2240 0
.87 68 1954 0
.54 79    2 0
.95 66  513 0
.24 77 2079 0
 .9 69 2205 0
.22 71  245 0
.22 88  903 1
.59 62 1408 0
.75 76 2189 0
.65 87  450 0
.55 74  884 0
.98 81 2079 0
 .8 75   72 0
 .9 72 1881 0
.75 83 2123 0
.35 70  315 1
.94 62 2172 0
.95 80 1766 0
.72 69 2043 0
.71 71 1956 0
.85 60 2033 0
 .8 70 2045 0
.84 81 2173 0
 .7 74  549 1
.07 69  736 1
.51 84 2138 0
.66 70 1759 0
.85 69 1742 0
 .6 69  335 1
.94 65 1787 0
.37 69 1816 0
.26 68  131 1
.09 71 1604 0
.72 73 1943 0
.78 76 2184 0
.69 82 1999 0
.86 48 1295 0
.89 63 1611 0
.54 63 1861 0
.82 69 1877 0
.27 74  996 1
.59 76    1 0
.89 64 1776 0
.75 59  699 0
.19 59 1513 0
.28 50  147 1
.93 69 1705 0
.88 82  347 0
.93 77 1705 0
.94 71 1752 0
 .8 88  999 1
.47 85   35 1
.84 78 1748 0
.86 73 1864 0
.13 78  814 0
.78 61    1 0
.18 85 1708 0
.94 63    1 0
.76 78 1939 0
.53 48 1984 0
.68 91 2174 0
.48 82 1720 0
.38 76 1833 0
.57 64 1922 1
.89 65 2141 0
.94 81 1969 0
.49 64    7 0
.77 59 1952 0
.89 76 1747 0
.83 64 1919 0
end
label values outcome yn
label def yn 0 "No", modify
label def yn 1 "Yes", modify

All is well when I run the univariable model

*univariable 
    *customs check
        stset followup, failure(outcome) show scale(30.4)                
        
    *run cox regression
        stcox predictor
        
    *predict hazard 
        predict xb_predictor, xb
        predict hr_predictor, hr
    
    *generate 95% ci
        predict se_predictor, stdp
        gen lowerb_predictor = exp(xb_predictor- invt(1205, 0.975)*se_predictor) 
        gen upperb_predictor = exp(xb_predictor + invt(1205, 0.975)*se_predictor)
        
    *sort
        sort predictor
    
    *generate figure
        twoway  (rarea lowerb_predictor upperb_predictor predictor, ///
                        fcolor(red) fintensity(inten10) ///
                        lcolor(black) lwidth(medium) lpattern(dash))  || ///
                (line hr_predictor predictor, lcolor(black) lwidth(thick)) || ///
                (scatteri 1 0 1 1, recast(line) ///
                            lcolor(red) lwidth(thin) ///
                            lpattern(dash) lpattern(dash)), ///
                    legend(off) ///
                    name(hr_univar, replace)
                    
        graph export ".\Incremental prognostic role of ML-FFRCT\Graphs/hr_univar.jpg", replace width(500)
                    
        capture drop xb* 
        capture drop hr* 
        capture drop lowerb* 
        capture drop upperb* 
        capture drop se*

this is the figure i get

however, when i try running a multivariable model adjusted to the confounder, this is what i get

*multivariable 
    *customs check
        stset followup, failure(outcome) show scale(30.4)                
        
    *run cox regression
        stcox predictor confounder
        
    *predict hazard 
        predict xb_predictor, xb
        predict hr_predictor, hr
    
    *generate 95% ci
        predict se_predictor, stdp
        gen lowerb_predictor = exp(xb_predictor- invt(1205, 0.975)*se_predictor) 
        gen upperb_predictor = exp(xb_predictor + invt(1205, 0.975)*se_predictor)
        
    *sort
        sort predictor
    
    *generate figure
        twoway  (rarea lowerb_predictor upperb_predictor predictor, ///
                        fcolor(red) fintensity(inten10) ///
                        lcolor(black) lwidth(medium) lpattern(dash))  || ///
                (line hr_predictor predictor, lcolor(black) lwidth(thick)) || ///
                (scatteri 1 0 1 1, recast(line) ///
                            lcolor(red) lwidth(thin) ///
                            lpattern(dash) lpattern(dash)), ///
                    legend(off) ///
                    name(hr_multivar, replace)
                    
        graph export ".\Incremental prognostic role of ML-FFRCT\Graphs/hr_multivar.jpg", replace width(500)
                    
        capture drop xb* 
        capture drop hr* 
        capture drop lowerb* 
        capture drop upperb* 
        capture drop se*

I've tried sorting with the predictor, confounder and the hazard ratio's but nothing seems to change

Please advise

[Clyde Schechter]

Well, your hazard ratio now depends on two variables, predictor and confounder, and these will occur in various combinations with each other. Some observations will have predictor equal to, say, 0.4, and confounder at some low value. But other observations will have predictor equal to 0.4, or nearly so, and confounder at some high value. The associated hr's will be quite different. So you no longer have hr as a function of predictor alone, and the graph is showing you that.

So, if you are looking for an average hr as a function of predictor, averaged over the joint distribution of predictor and confounder, you can do that with the -margins- command, and then get the graph with -marginsplot-. (The default graph after -marginsplot- may not be to your liking, but you can apply almost any -graph twoway- option to -marginsplot- to customize the look. If that doesn't get you what you want (I'm not sure whether it will give you an area plot like that) -marginsplot- also has a -saving()- option that you will save all its results in a Stata data set, and you can then -use- that data set with -graph twoway- commands and option. Note that for uncertainty around the hr's that -margins- gives you, -margins- computes its own standard error. You will not be able to use the results of -predict, stdp- with -margins-.

If you are not familiar with the -margins- command, I suggest you start out by reading the excellent Richard Williams' https://www3.nd.edu/~rwilliam/stats/Margins01.pdf. I think it provides the clearest explanation of how -margins- works and what it does.

[Ahmed Ibrahim]

Thank you very much Clyde. This was very helpful.

I used the margins command and was able to generate the graph.

However, the 95% CI is very important to the analysis I'm running as I want to show at what point the lower 95%CI crosses a hazard ratio of 1.

Is there any other way I would be able to get the 95%CI's from the margins command?

[Andrew Lover]

Check out flex parametric models (stpm2 from SSC), which allow predictions at covariates with 95% CIs.

https://www.stata.com/meeting/nordic...18_Lambert.pdf