version 9.0
set more on
log using 7.9.4.Framingham.modelfit.log, replace
*  Framingham.ModelFit.log
*
*  Check the proportional hazards model fit for the effect of gender on
*  CHD with time measured as age.
*
use C:\WDDtext\2.20.Framingham.dta, clear
more
generate time= followup/365.25
label variable time "Follow-up in Years"
generate exitage = age+time
stset exitage, enter(time age) failure(chdfate) 
sts graph, by(sex) failure ytitle(Cumulative CHD Morbidity) xtitle(Age)	///
    ylabel(0(.1).8, angle(0)) legend(ring(0) position(11) col(1))	/// 
    color(blue pink) clwidth(medthick medthick) xlabel(30(10)90) noorigin
more
stphplot, by(sex) nolntime 						///
    legend(ring(0) position(2) col(1) ) 				///
    msymbol(none none) color(red blue )
more
stcoxkm, by(sex) legend(ring(0) position(7) col(1)) 			///
        msymbol(none none none none )					///
        lpattern(solid solid dash dash)					///
        color(red blue red blue)
more
stcox sex,  basehc(h)
stcurve, at1(sex=1) at2(sex=2) hazard yscale(log) outfile(plot, replace)
more
drop h
stcox if sex==1, estimate  basehc(h)
stcurve, hazard yscale(log) outfile(plotmen, replace)
more
drop h
stcox if sex==2, estimate  basehc(h)
stcurve, hazard yscale(log) outfile(plotwomen, replace)
use plotwomen, clear
list in 1/3
rename haz1 hazwomen
save plotwomen, replace
use plotmen, clear
list in 1/3
rename haz1 hazmen
save plotmen, replace
use plot, clear
list in 1/3
merge _t using plotmen plotwomen, sort
sort _t
list in 1/3
line        haz2 _t, color(blue) lpattern(solid)	    ///
    || line haz3 _t, color(pink) lpattern(solid)	    ///
    || line hazmen _t, color(blue) lpattern(dash)	    ///
    || line hazwomen _t, color(pink) lpattern(dash)	    ///
    , legend(ring(0) position(5) col(1) 		    ///
          order(1 "Men:	      proportional hazards" 	    ///
                2 "Women: proportional hazards" 	    ///
                3 "Men:       unrestricted"  		    ///
                4 "Women: unrestriced")) 		    ///
      ylabel(.001 .002 .004 .006 .01 .02 .04 .06, angle(0)) ///
      ytick(.003(.001).009 .03 .05) yscale(log)		    ///
      xtitle(Age) ytitle(Smoothed hazard function) 
more
*
*  Analyze patients less than 47 years old
*
use C:\WDDtext\2.20.Framingham.dta, clear
generate time= followup/365.25
label variable time "Follow-up in Years"
generate exitage = age+time
stset exitage, enter(time age) failure(chdfate) 
generate chdfate47 = chdfate & exitage <= 47
generate exitage47 = min(47,exitage)
generate chdfate57 = chdfate & exitage <= 57
generate exitage57 = min(57, exitage)
generate age47 = max(age, 47)
generate age57 = max(age, 57)
stset exitage47, enter(time age) failure(chdfate47)
stcox sex if age <47,  basehc(h)
stcurve, at1(sex=1) at2(sex=2) hazard yscale(log) outfile(plot47, replace)
more
stset exitage57, enter(time age47) failure(chdfate57)
drop h
stcox sex if age < 57 & exitage >47,  basehc(h)
stcurve, at1(sex=1) at2(sex=2) hazard yscale(log) outfile(plot57, replace)
more
stset exitage, enter(time age57) failure(chdfate)
drop h
stcox sex if exitage >=57,  basehc(h)
stcurve, at1(sex=1) at2(sex=2) hazard yscale(log) outfile(plotold, replace)
more
use plotold, replace
merge _t using plotmen plotwomen plot47 plot57, sort
sort _t
line        haz2 _t, color(blue) lpattern(solid)	    ///
    || line haz3 _t, color(pink) lpattern(solid)	    ///
    || line hazmen _t, color(blue) lpattern(dash)	    ///
    || line hazwomen _t, color(pink) lpattern(dash)	    ///
    , legend(ring(0) position(5) col(1) 		    ///
          order(1 "Men:	      time-dependent model" 	    ///
                2 "Women: time-dependent model" 	    ///
                3 "Men:       unrestricted"  		    ///
                4 "Women: unrestriced")) 		    ///
      ylabel(.001 .002 .004 .006 .01 .02 .04 .06, angle(0)) ///
      ytick(.003(.001).009 .03 .05) yscale(log)		    ///
      xtitle(Age) ytitle(Smoothed hazard function) 

log close