version 9.0
set more on
log using 7.9.4.Framingham.modelfit2.log, replace
*  Framingham.ModelFit.log
*
*  Check the proportional hazards model fit for the effect of gender on
*  CHD with time measured as either time since recruitment.
*
use C:\WDDtext\2.20.Framingham.dta, clear
more
generate time= followup/365.25
label variable time "Follow-up in Years"
stset time,  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))	/// 
*    clcolor(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(.003 .004 (.002) .01 .02 .03, angle(0)) 	    ///
      ytick(.004(.001).009 ) yscale(log)		    ///
      xtitle(Years since recruitment) ytitle(Smoothed hazard function) 
    
log close