log using support.log, replace
*
*  Perform linear and logistic regressions of log length of stay and 
*  hospdead against meanbp, respectively.  Use restricted cubic splines 
*  for these regessions.
*
use support.dta, replace
set more on
*
*  Analyze relationship between entry mean BP and length of stay
*
scatter los meanbp							///
    , msymbol(Oh) xlabel(25(25)175) xtick(30 (5) 170) name(scatter, replace)
more
*
*  Define logarithmic tick marks and labels for y-axis of length of stay
*
local ylogtick = " "
foreach i of numlist 3/10 20(10)100 200 { 
    local log`i' = log(`i')
    local ylogtick "`ylogtick' `log`i''"
    }
local yloglabel " "
foreach i of numlist 4(2)10 20(20)100 200 {
    local yloglabel "`yloglabel' `log`i'' "`i'""
    }
gen log_los = log(los)
more
rc_spline meanbp
regress log_los _S*
test _Smeanbp2  _Smeanbp3 _Smeanbp4
predict y_hat, xb
scatter  log_los meanbp ,msymbol(Oh) 					///	
    || line y_hat meanbp						///
    ,  xlabel(25 (25) 175) xtick(30 (5) 170) clcolor(red) 		///
       clwidth(thick) xline(47 66 78 106 129, lcolor(blue))		///
       ylabel(`yloglabel', angle(0)) ytick(`ylogtick') 			///
       ytitle("Length of Stay (days)") 					///
       legend(order(1 "Observed" 2 "Expected")) name(knot5, replace)
more
drop _S* y_hat
rc_spline meanbp, nknots(7)
regress log_los _S*
test _Smeanbp2  _Smeanbp3 _Smeanbp4  _Smeanbp5  _Smeanbp6
predict y_hat, xb
scatter  log_los meanbp ,msymbol(Oh) 					///
    || line y_hat meanbp						///
    ,  xlabel(25 (25) 175) xtick(30 (5) 170) clcolor(red)  		///
       clwidth(thick) xline(41 60 69 78 101 113 138, lcolor(blue))	///
       ylabel(`yloglabel', angle(0)) ytick(`ylogtick') 			///
       ytitle("Length of Stay (days)") 					///
       legend(order(1 "Observed" 2 "Expected")) name(knot7, replace)
more
drop _S* y_hat
rc_spline meanbp, nknots(7) knots(40(17)142)
regress log_los _S*
predict y_hat, xb
scatter  log_los meanbp ,msymbol(Oh) 					///
    || line y_hat meanbp						///
    ,  xlabel(25 (25) 175) xtick(30 (5) 170) clcolor(red) 		///
       clwidth(thick) xline(40(17)142, lcolor(blue))			///
       ylabel(`yloglabel', angle(0)) ytick(`ylogtick') 			///
       ytitle("Length of Stay (days)") 					///
       legend(order(1 "Observed" 2 "Expected")) name(evenknots, replace)
more
drop _S* y_hat
rc_spline meanbp, nknots(7)
regress log_los _S*
predict se, stdp
predict y_hat, xb
generate lb = y_hat - invttail(_N-7, 0.025)*se
generate ub = y_hat + invttail(_N-7, 0.025)*se
twoway rarea lb ub meanbp , bcolor(gs6) lwidth(none)			///
    || scatter  log_los meanbp ,msymbol(Oh) mcolor(blue) 		///
    || line y_hat meanbp, xlabel(25 (25) 175) xtick(30 (5) 170) 	///
       clcolor(red) clwidth(thick) ytitle("Length of Stay (days)") 	///
       ylabel(`yloglabel', angle(0)) ytick(`ylogtick') name(ci,replace)	///
       legend(rows(1) order(2 "Observed" 3 "Expected" 1 "95% CI" ))        
more
predict rstudent, rstudent
lowess rstudent meanbp							///
    , yline(-2 0 2) msymbol(Oh) rlopts(clcolor(green) clwidth(thick))	///
    xlabel(25 (25) 175) xtick(30 (5) 170) 
more
gen big = abs(rstudent)>2
tabulate big
*
*  Analyze relationship between hospital death and length of stay
*  First, calculate rate to equal the proportion of patients with 
*  a given meanbp who die in hospital.
*
sort meanbp
by meanbp: egen rate = mean(hospdead)
label variable rate "Hospital Mortality Rate"
generate deads = meanbp if hospdead
*
*  Draw an exploratory graph showing the # of patients, the # of deaths
*  and the mortality rate for each mean BP
*
twoway histogram meanbp, discrete frequency blwidth(none) gap(20)	///
    || histogram deads, discrete frequency bcolor(red) blwidth(none) 	///
           gap(20)							///
    || scatter rate meanbp, yaxis(2) mcolor(blue) msymbol(Oh) 		///
    ,  ylabel(0 (5) 30, angle(0) axis(1)) xlabel(20 (20) 180)		///
       xtick(25 (5) 175) ytitle("Number of Patients")			///
       ytitle(, axis(2) color(blue) )  					///
      ylabel(0 (.1) 1, angle(0) axis(2) labcolor(blue)) 		///
      ytick(0(.05)1, axis(2) tlcolor(blue)) 				///
      legend(order(1 "Total" 2 "Deaths" 3 "Mortality Rate") rows(1)) 
more
*
*  Regress hospdead against meanbp using simple logistic regression
*
logistic hospdead meanbp
predict p,p
line p meanbp, ylabel(0 (.1) 1) ytitle(Probabilty of Hospital Death)
more
drop p
*
*  Repeat the preceding model using restricted cubic splines with 5 knots
*
drop _S* 
rc_spline meanbp
logistic hospdead _S*, coef
*
*  Plot the predicted probability of death against meanbp together with
*  the 95% confidence interval for this curve and the observed mortality
*  for each mean bp.
*
predict p,p
predict logodds, xb
predict stderr, stdp
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 meanbp, bcolor(gs14)				///
    || line p meanbp, clcolor(red) clwidth(medthick)			///
    || scatter rate meanbp, msymbol(Oh)	mcolor(blue)			///
    ,  ylabel(0 (.1) 1, angle(0)) xlabel(20 (20) 180) 			///
       xtick(25 (5) 175) ytitle(Probabilty of Hospital Death)  		///
       legend(order(3 "Observed Mortality"  				///
            2 "Expected Mortality"  1 "95% CI") rows(1))
test _Smeanbp2 _Smeanbp3 _Smeanbp4
*
*  Calculate odds ratios for meanbp =  60 vs meanbp = 90
*                            meanbp = 120 vs meanbp = 90
*
list _S* logodds p 							///
    if (meanbp==60 | meanbp==90 | meanbp==120) & meanbp ~= meanbp[_n-1]
    
lincom (5.531072 +       60*_Smeanbp1 +   .32674*_Smeanbp2  		///
                 +        0*_Smeanbp3 +       0 *_Smeanbp4)		///
                 							///
      -(5.531072 +       90*_Smeanbp1 + 11.82436*_Smeanbp2 		///
                 + 2.055919*_Smeanbp3 + .2569899*_Smeanbp4)
                 
lincom (5.531072 +      120*_Smeanbp1 + 56.40007*_Smeanbp2 		///
		 + 22.30039*_Smeanbp3 + 10.11355*_Smeanbp4) 		///
		 							///
      -(5.531072 +       90*_Smeanbp1 + 11.82436*_Smeanbp2 		///
                 + 2.055919*_Smeanbp3 + .2569899*_Smeanbp4)
log close