-------------------------------------------------------------------------------------------------------------------
       log:  C:\MyDocs\MPH\LectureNotes\ClassDoLogData\logistic regression\SepsisSplineCheck.log
  log type:  text
 opened on:  14 Feb 2008, 14:26:20

. *  SepsisSplineCheck.do
. *
. *  Logistic regression of mortality against APACHE score in
. *  in the Ibuprofen in Sepsis Study using restricted cubic splines.  
. *
. use C:\WDDtext\4.18.Sepsis.dta, clear

. 
. generate p = deaths/patients

. generate obs_logodds = logit(p) if p ~=0 & p ~= 1
(11 missing values generated)

. *
. *  Regress deaths against apache using simple logistic regression
. *
. 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
------------------------------------------------------------------------------

. predict logodds, xb

. *
. *  Calculate 95% confidence region for logodds
. *
. predict stderr, stdp

. generate lodds_lb = logodds - 1.96*stderr

. generate lodds_ub = logodds + 1.96*stderr

. *
. *  Graph observed logodds with regression line & confidence region
. *
. twoway rarea lodds_lb lodds_ub apache,  bcolor(yellow)                  ///
>     || line logodds apache                                              ///
>     || scatter obs_logodds apache, ytitle("Log Odds of Death")          ///
>     ,  legend(order(2 "Simple log odds" 3 "Observed log odds"           ///
>             1 "95% CI") rows(1))

. more

. *
. *  Calculate rc_spline covariates with 3 knots
. *
. mkspline _Sapache=apache [weight=patients], nknots(3) cubic displayknots
(frequency weights assumed)

             |     knot1      knot2      knot3 
-------------+---------------------------------
      apache |         7         14         25 

. *
. *  Regress deaths against spline covariates
. *
. glm deaths _S*, family(binomial patients) link(logit)

Iteration 0:   log likelihood = -61.241346  
Iteration 1:   log likelihood = -60.915581  
Iteration 2:   log likelihood = -60.914709  
Iteration 3:   log likelihood = -60.914709  

Generalized linear models                          No. of obs      =        38
Optimization     : ML                              Residual df     =        35
                                                   Scale parameter =         1
Deviance         =  43.96057815                    (1/df) Deviance =  1.256017
Pearson          =  45.82765759                    (1/df) Pearson  =  1.309362

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

                                                   AIC             =  3.363932
Log likelihood   = -60.91470918                    BIC             = -83.35494

------------------------------------------------------------------------------
             |                 OIM
      deaths |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   _Sapache1 |   .1237794   .0447174     2.77   0.006      .036135    .2114238
   _Sapache2 |  -.0116944   .0596984    -0.20   0.845     -.128701    .1053123
       _cons |  -2.375381   .5171971    -4.59   0.000    -3.389069   -1.361694
------------------------------------------------------------------------------

. predict logodds3, xb

. line logodds apache,  lwidth(thick) color(blue)                 ///
>     || line logodds3 apache,  lwidth(thick) color(red)                  ///
>     || scatter obs_logodds apache, ytitle("Log Odds of Death")          ///
>     ,  legend(order(1 "Simple log odds" 2 "Spline log odds"             ///
>             3 "Observed log odds") rows(1))

. more

. predict e_deaths, mu

. generate prob = e_deaths/patients

. predict se, stdp

. generate lb_logodds3= logodds3-1.96* se

. generate ub_logodds3= logodds3+1.96* se

. generate ub_prob= exp( ub_logodds3)/(1+exp( ub_logodds3))

. generate lb_prob= exp( lb_logodds3)/(1+exp( lb_logodds3))

. twoway rarea lb_prob ub_prob apache, bcolor(yellow)                     ///
>     || line prob apache, color(blue) lwidth(medthick)                   ///
>     || scatter p apache, ytitle("Probability of Death by 30 Days")      ///
>     legend(order(2 "Expected mortality rate"                            ///
>            3 "Observed mortality rate" 1 "95% CI") rows(1)) subtitle(   ///
>        "Restricted Cubic Spline Model with 3 Knots",ring(0) position(11))

. log close
       log:  C:\MyDocs\MPH\LectureNotes\ClassDoLogData\logistic regression\SepsisSplineCheck.log
  log type:  text
 closed on:  14 Feb 2008, 14:26:39
----------------------------------------------------------------------------------------------------------------