version 9.0
set more on
log using SepsisExactCI.log, replace
*  in progress
*  4.18.Sepsis.do
*
*  Simple logistic regression of mortality against APACHE score in
*  in the Ibuprofen in Sepsis Study.  Each record of 4.18.Sepsis.dta
*  gives the number of patients and number of deaths among patients
*  with a specified apache score. These variables are named patients
*  deaths and apache, respectively.
*
use 4.18.Sepsis.dta, clear
*
*  Calculate 95% confidence intervals for observed mortality rates
*
generate p = deaths/patients
generate m = patients
generate q = 1-p
generate c = 1.96
generate ci95lb = (2*m*p+c^2-1-c*sqrt(c^2-(2+1/m)+4*p*(m*q+1)))/(2*(m+c^2)) ///
    if deaths ~=0 & deaths ~= m
generate ci95ub = (2*m*p+c^2+1+c*sqrt(c^2+(2-1/m)+4*p*(m*q-1)))/(2*(m+c^2)) ///
    if deaths ~=0 & deaths ~= m
label variable ci95ub "95% CI ub"
label variable ci95lb "lb obs."
*
*  Regress deaths against apache
*
glm deaths apache, family(binomial patients) link(logit)
vce
predict logodds, xb
generate e_prob = exp(logodds)/(1+exp(logodds))
label variable e_prob "Expected Mortality Rate"
*
*  Calculate 95% confidence region for e_prob
*
predict stderr, stdp
generate lodds_lb = logodds - 1.96*stderr
generate lodds_ub = logodds + 1.96*stderr
generate prob_lb = exp(lodds_lb)/(1+exp(lodds_lb))
generate prob_ub = exp(lodds_ub)/(1+exp(lodds_ub))
label variable p "Observed Mortality Rate"
label variable prob_ub "95% CI ub"
label variable prob_lb "lb exp."
list p e_prob prob_lb prob_ub ci95lb ci95ub apache if apache == 20
twoway rarea prob_ub prob_lb apache, color(yellow) 			///
    || scatter p apache, color(blue)					///
    || rcap ci95ub ci95lb apache, color(blue)			  	///
    || line e_prob apache, yaxis(2)  clwidth(medthick) color(red)	///
    ,  ylabel(0(.2)1, labcolor(blue) angle(0)) 				///
       ytick(0(.1)1, tlcolor(blue))					///
       ylabel(0(.2)1, axis(2) labcolor(red) angle(0)) 			///
       ytick(0(.1)1, axis(2) tlcolor(red))				///
       xlabel(0(5)40) xtick(0(1)40)					///
       ytitle(,axis(2) color(red)) 					///
       ytitle(Observed Mortality Rate, color(blue))			///
       legend(order(1 "95% CI from model" 2 3 "95% CI from this obs." 4))
more
histogram apache [freq=patients], discrete frequency  xlabel(0(5)40) 	///
    ylabel(0(5)30, angle(0)) ytitle(Number of Patients)
log close