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