--------------------------------------------------------------------------------------------------------
      name:  <unnamed>
       log:  C:\MyDocs\MPH\LectureNotes\ClassDoLogData\logistic regression\4.18.Sepsis.log
  log type:  text
 opened on:  23 Apr 2010, 17:00:42

. *  4.18.Sepsis.do
. *
. *  Simple logistic regression of mortality against APACE score in the 
. *  Ibuprofen in Sepsis study (Bernard et al. 1997).  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 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 = 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"

. 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, 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)
(start=0, width=1)

. 
. log close
      name:  <unnamed>
       log:  C:\MyDocs\MPH\LectureNotes\ClassDoLogData\logistic regression\4.18.Sepsis.log
  log type:  text
 closed on:  23 Apr 2010, 17:00:54
--------------------------------------------------------------------------------------------------------