window push version 11 window push set more on window push log using Framingham.Spline.log, replace window push * Framingham.Spline.log window push * window push * Proportional hazards regression analysis of the effect of gender and window push * baseline diastolic blood pressure (DBP) on coronary heart disease (CHD) window push * using restricted cubic splines to model the effect of DBP on CHD risk. window push * We will use a DBP of 60 as the denominator of our relative risk estimates. window push * window push set more on window push use C:\WDDtext\2.20.Framingham.dta, clear window push generate time= followup/365.25 window push label variable time "Follow-up in Years" window push stset time, failure(chdfate) window push sort dbp window push generate dbp60 = dbp - 60 window push mkspline _Sdbp60 = dbp60, cubic displayknots window push stcox _S*, nohr window push test _Sdbp602 _Sdbp603 _Sdbp604 window push predict relhaz5, hr window push more window push * window push * Experiment with fewer knots window push * window push drop _S* window push mkspline _Sdbp60 = dbp60, nknots(3) cubic displayknots window push stcox _S*, nohr window push predict relhaz3, hr window push * window push * How about no knots? window push * window push stcox dbp60 window push predict relhaz0, hr window push drop _S* window push summarize dbp60, detail window push * window push * Add a knot at DBP60 = 60 and remove the knot at DBP60 = 8 window push * window push mkspline _Sdbp60 = dbp60, knots(20 40 60) cubic displayknots window push stcox _S*, nohr window push predict relhaz3a, hr window push * window push * Calculate the relative hazard from model 7.12 in the text window push * window push gen relhazcat = 1 window push replace relhazcat = 1.97 if dbp > 60 window push replace relhazcat = 2.56 if dbp > 70 window push replace relhazcat = 3.06 if dbp > 80 window push replace relhazcat = 4.54 if dbp > 90 window push replace relhazcat = 6.29 if dbp > 100 window push replace relhazcat = 9.46 if dbp > 110 window push * window push * Plot relative hazards estimated so far window push * window push line relhaz0 relhaz3 relhaz3a relhaz5 dbp /// window push , color(blue green red purple) /// window push || line relhazcat dbp, connect(stepstair) color(gray) /// window push , legend(ring(0) position(11) col(1) /// window push order(1 "No knots" 2 "3 default knots" /// window push 3 "3 special knots" 4 "5 default knots" /// window push 5 "Categorical")) ytitle(Relative risk) /// window push ytick(1(1)16) ylabel(0(3)15, angle(0)) window push more window push predict loghaz, xb window push predict se, stdp window push generate logcil = loghaz - 1.96*se window push generate logciu = loghaz +1.96*se window push twoway rarea logcil logciu dbp, color(yellow) /// window push || line loghaz dbp, color(red) /// window push , legend(ring(0) position(11) col(1) /// window push order(1 "95% CI: 3 special knots" /// window push 2 "3 special knots")) ytitle(Log relative risk) window push more window push generate cil3a = exp(logcil) window push generate ciu3a = exp(logciu) window push twoway rarea cil3a ciu3a dbp, color(yellow) /// window push || line relhaz3a dbp, color(red) /// window push , legend(ring(0) position(11) col(1) /// window push order(1 "95% CI: 3 special knots" /// window push 2 "3 special knots" )) ytitle(Relative risk) window push more window push * window push * Plot results from the no knot model and the preceding window push * model together. Truncate the upper error bounds. window push * window push stcox dbp60 window push drop loghaz se logcil logciu window push predict loghaz, xb window push predict se, stdp window push generate logcil = loghaz - 1.96*se window push generate logciu = loghaz + 1.96*se window push generate cil0 = exp(logcil) window push generate ciu0 = exp(logciu) window push egen maxhaz = max(relhaz0) window push generate ciu3a_chop = min(ciu3a,maxhaz) window push generate ciu0_chop = min(ciu0,maxhaz) window push twoway rarea cil3a ciu3a_chop dbp, color(yellow) /// window push || rarea cil0 ciu0_chop dbp, color(gs14) /// window push || line relhaz3a dbp, color(red) /// window push || line relhaz0 dbp, color(blue) /// window push , legend(ring(0) position(11) col(1) /// window push order(1 "95% CI: 3 special knots" /// window push 2 "95% CI: no knots" 3 "3 special knots" /// window push 4 "No knots")) ytitle(Relative risk of CHD) /// window push ytick(1(1)16) ylabel(0(3)15, angle(0)) window push more window push * window push * In our final graphs we will want to truncate the upper window push * error bands at the top of the graph. This can cause window push * linear extrapolation errors due to sparce blood pressures window push * at the extreme upper range. To correct this we add window push * dummy records to fill in some of these blood pressures. window push * window push set obs 4739 window push replace dbp = 135 +(_n - 4699)*0.1 if _n > 4699 window push replace dbp60 = dbp - 60 window push sort dbp window push drop loghaz se logciu maxhaz ciu0 window push predict loghaz, xb window push predict se, stdp window push generate logciu = loghaz +1.96*se window push generate ciu0 = exp(logciu) window push egen maxhaz = max(relhaz0) window push replace ciu0_chop = min(ciu0,maxhaz) window push twoway rarea cil3a ciu3a_chop dbp, color(yellow) /// window push || rarea cil0 ciu0_chop dbp, color(gs14) /// window push || line relhaz3a dbp, color(red) /// window push || line relhaz0 dbp, color(blue) /// window push , legend(ring(0) position(11) col(1) /// window push order(1 "95% CI: 3 special knots" /// window push 2 "95% CI: no knots" 3 "3 special knots" /// window push 4 "No knots")) ytitle(Relative risk of CHD) /// window push ytick(1(1)16) ylabel(0(3)15, angle(0)) window push log close