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:
perform cox regression and test non-linearity:
Create estimates for graph:
Graph the line:
Graph:
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)
Code:
stcox age_spline* treatment var1 var2 var3 var4 var5 testparm age_spline2 age_spline3 age_spline4
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)
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:
