Dear Statalist,

I want to plot the hazard ratio for survival over age using a restricted cubic splines model. However, rather than the line representing the HR of each age on survival, I want the line to represent the hazard ratio of treatment A vs. treatment B at each age. In other words, the line should represent at which treatment should be preferred at every age (as < 1 favors A, ≥ 1 favors treatment B). I know this should be possible with Stata, since I found an article that describes in the methods that they use the postrcspline package (article: https://doi.org/10.1016/j.jvs.2018.07.062). I have been trying for weeks now, but I'm only able to plot the HRs of age on survival.

My guess was I should use an interaction term with age and treatment, but the packages I tried do not seem to support interactions. If anyone would have a suggestion on how to edit the code below (or how to do it in a different way), that would be greatly appreciated.

This is the code I used to come to the HR for age on survival:

Create splines:
Code:
mkspline age_spline = age, cubic displayknots
mat knots = r(knots)
perform cox regression and test non-linearity:
Code:
stcox age_spline* treatment var1 var2 var3 var4 var5 
testparm  age_spline2 age_spline3 age_spline4
Create estimates for graph:
Code:
levelsof age
xbrcspline age_spline, values(`r(levels)') ref(65) matknots(knots) eform
predictnl xb= _b[age_spline1]*(age_spline1-65)+ _b[age_spline2]*(age_spline2-0)+ _b[age_spline3]*(age_spline3-0)+ _b[age_spline4]*(age_spline4-0), ci(lo hi)

gen hr = exp(xb)
gen lb = exp(lo)
gen ub = exp(hi)
Graph the line:
Code:
sort age
graph tw line hr age, color(navy) ylabel(0.2 0.3 0.4 0.5 0.7 1 1.5 2 3 5) /// 
    || rarea lb ub age, color(navy%30) ///
    yline(1) yscale(log) ytitle("Hazard Ratio") xtitle("Age") ///
    xlabel(20 (5) 90) title({bf:Effect of age on survival})


Graph:
Array