version 9.0
set more on
set scheme s2color
log using SUPPORTlogisticRCS.log, replace
* SUPPORTlogisticRC.log
*
*  Regress mortal status at discharge against MAP 
*  in the SUPPORT data set
*
set more on
use 3.25.2.SUPPORT.dta , replace
*
*  Calculate the proportion of patients who die in hospital
*   stratified by MAP.
*
generate map_gr = round(map,5)
sort map_gr
label variable map_gr "Mean Arterial Pressure (mm Hg)"
by map_gr: egen proportion = mean(fate)
generate rate = 100*proportion
label variable rate "Observed In-Hospital Mortality Rate (%)"
generate deaths = map_gr if fate
*
*  Draw an exploratory graph showing the number of patients, 
*  the number of deaths and the mortality rate for each MAP.
*
twoway histogram map_gr, discrete frequency color(gs13) gap(20) ///
    || histogram deaths, discrete frequency color(red) gap(20)	///
    || scatter rate map_gr, yaxis(2) symbol(Oh) color(blue)	///
    ,  xlabel(20 (20) 180) ylabel(0(10)100, angle(0))		///
       xmtick(25 (5) 175) ytitle(Number of Patients)		///
       ylabel(0 (10) 100, angle(0) labcolor(blue) axis(2)) 	///
       ytitle(,color(blue) axis(2))				///
       legend(order(1 "Total" 2 "Deaths" 3 "Mortality Rate" )	///
       rows(1)) name(g1) 
more
*
*  Regress in-hospital mortality against MAP using simple 
*  logistic regression.
*
logistic fate map
estimates store simple
predict p,p
label variable p "Probabilty of In-Hospital Death"
line p map, ylabel(0(.1)1, angle(0)) xlabel(20(20)180) name(g2)
more
drop p
*
*  Repeat the preceding model using restricted cubic splines 
*  with 5 knots at their default locations.
*
mkspline _Smap = map, cubic displayknots
logistic fate _S*
*
*  Test null hypotheses that the logit of the probability of
*  in-hospital death is a linear function of MAP.
*
lrtest simple .
display 2*(545.25721 -498.65571)
*
*  Plot the estimated probability of death against MAP together 
*  with the 95% confidence interval for this curve.  Overlay 
*  the MAP-specific observed mortality rates.
*
predict p,p
predict logodds, xb
predict stderr, stdp
generate p2 = exp(logodds)/(1+exp(logodds))
*
*  The values of p and p2 are identical
*
scatter p p2
generate lodds_lb = logodds - 1.96*stderr
generate lodds_ub = logodds + 1.96*stderr
generate ub_p= exp(lodds_ub)/(1+exp(lodds_ub))
generate lb_p= exp(lodds_lb)/(1+exp(lodds_lb))
twoway rarea lb_p ub_p map, color(yellow)				///
    || line p map, lwidth(medthick) color(red)				///
    || scatter proportion map_gr,  symbol(Oh) color(blue)		///
    ,  ylabel(0(.1)1, angle(0)) xlabel(20 (20) 180)  			///
       xmtick(25(5)175) ytitle(Probabilty of In-Hospital Death) 	///
       legend(order(3 "Observed" "Mortality" 2 "Expected" "Mortality" 	///
           1 "95% Confidence" "Interval") rows(1)) name(g3)
more
*
*  Determine the spline covariates at MAP = 90
*
list _S* if map == 90
*
*  Let or1 = _Smap1 minus the value of _Smap1 at 90.
*  Define or2, or3 and or3 in a similar fashion.
*
generate or1 = _Smap1 - 90
generate or2 = _Smap2 - 11.82436 
generate or3 = _Smap3 -  2.055919
generate or4 = _Smap4 -   .2569899
*
*  Calculate the log odds ratio for in-hospital death  
*  relative to patients with MAP = 90.
*
predictnl log_or = or1*_b[_Smap1] + or2*_b[_Smap2] ///
    + or3*_b[_Smap3] +or4*_b[_Smap4], se(se_or)
generate lb_log_or = log_or - 1.96*se_or
generate ub_log_or = log_or + 1.96*se_or
generate or = exp(log_or)
generate lb_or = exp(lb_log_or)
generate ub_or = exp(ub_log_or)

twoway rarea lb_or ub_or map, color(yellow)		 ///
    || line or map, lwidth(medthick) color(red)		 ///		
    ,  ylabel(1 (3) 10 40(30)100 400(300)1000, angle(0)) ///
       ymtick(2(1)10 20(10)100 200(100)900) yscale(log)  ///
       xlabel(20 (20) 180) xmtick(25 (5) 175) 		 ///
       ytitle(In-Hospital Mortal Odds Ratio)  		 ///
       legend(ring(0) position(2) order(2 "Odds Ratio" 	 ///
           1 "95% Confidence Interval") cols(1)) name(g4)
more
log close