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