--------------------------------------------------------------------------------------------------------------------------
       log:  C:\MyDocs\MPH\Text\SecondEdition\WebDoFiles\4.18.Sepsis.log
  log type:  text
 opened on:  14 Jan 2008, 12:18:34

. *  4.18.Sepsis.do
. *
. *  Simple logistic regression of mortality against APACHE score 
. *  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 C:\WDDtext\4.18.Sepsis.dta, clear

. *
. *  Calculate observed mortality rate and 
. *  exact 95% confidence intervals for observed mortality rates.
. *
. generate p = deaths/patients

. generate ci95ub = invbinomial(patients, deaths, 0.025)
(7 missing values generated)

. generate ci95lb = 1 - invbinomial(patients, patients-deaths, 0.025)
(4 missing values generated)

. *
. *  Regress deaths against apache
. *
. glm deaths apache, family(binomial patients) link(logit)

Iteration 0:   log likelihood = -61.188504  
Iteration 1:   log likelihood = -60.934124  
Iteration 2:   log likelihood = -60.933906  
Iteration 3:   log likelihood = -60.933906  

Generalized linear models                          No. of obs      =        38
Optimization     : ML                              Residual df     =        36
                                                   Scale parameter =         1
Deviance         =  43.99897135                    (1/df) Deviance =  1.222194
Pearson          =  46.72842945                    (1/df) Pearson  =  1.298012

Variance function: V(u) = u*(1-u/patients)         [Binomial]
Link function    : g(u) = ln(u/(patients-u))       [Logit]

                                                   AIC             =  3.312311
Log likelihood   = -60.93390578                    BIC             = -86.95413

------------------------------------------------------------------------------
             |                 OIM
      deaths |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      apache |   .1156272   .0159997     7.23   0.000     .0842684    .1469861
       _cons |  -2.290327   .2765286    -8.28   0.000    -2.832314   -1.748341
------------------------------------------------------------------------------

. estat vce

Covariance matrix of coefficients of glm model

             | deaths                 
        e(V) |     apache       _cons 
-------------+------------------------
deaths       |                        
      apache |  .00025599             
       _cons | -.00410256   .07646805 

. predict logodds, xb

. generate e_prob = invlogit(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 = invlogit(lodds_lb)

. generate prob_ub = invlogit(lodds_ub)

. list p e_prob prob_lb prob_ub ci95lb ci95ub apache if apache == 20

     +------------------------------------------------------------------------+
     |        p    e_prob    prob_lb   prob_ub     ci95lb     ci95ub   apache |
     |------------------------------------------------------------------------|
 20. | .4615385   .505554   .4462291   .564723   .1922324   .7486545       20 |
     +------------------------------------------------------------------------+

. twoway rarea prob_ub prob_lb apache, color(yellow)                      ///
>    || scatter p apache, symbol(circle) color(blue)                      ///        
>    || rcap ci95ub ci95lb apache ,color(blue)                            ///
>    || line e_prob apache, yaxis(2) lwidth(thick) color(red)             ///
>    , xlabel(0(5)40) xmtick(0(1)41)                                      ///
>      ylabel(0(.1)1,angle(0) labcolor(blue) tlcolor(blue))               ///
>      ytitle("Observed Mortality Rate", color(blue))                     ///
>      ylabel(0(.1)1, axis(2) angle(0) labcolor(red) tlcolor(red))        ///
>      ytitle(,axis(2) color(red) ) legend(off)                   

. more

. histogram apache [freq=patients]                                        ///
>     , discrete frequency ytitle(Number of Patients)                     ///
>     xlabel(0(5)40) xmtick(1(1)41) ylabel(0(5)30) ymtick(1(1)33) 
(start=0, width=1)

. more

. log close
       log:  C:\MyDocs\MPH\Text\SecondEdition\WebDoFiles\4.18.Sepsis.log
  log type:  text
 closed on:  14 Jan 2008, 12:18:43
-----------------------------------------------------------------------------------------------------------------------