window push version 8.0
window push set more on
window push log using SepsisSplineCheck.log, replace
window push *  SepsisSplineCheck.do
window push *
window push *  Logistic regression of mortality against APACHE score in
window push *  in the Ibuprofen in Sepsis Study using restricted cubic splines.  
window push *
window push use C:\WDDtext\4.18.Sepsis.dta, clear
window push 
window push generate p = deaths/patients
window push generate obs_logodds = logit(p) if p ~=0 & p ~= 1
window push *
window push *  Regress deaths against apache using simple logistic regression
window push *
window push glm deaths apache, family(binomial patients) link(logit)
window push predict logodds, xb
window push *
window push *  Calculate 95% confidence region for logodds
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 *
window push *  Graph observed logodds with regression line & confidence region
window push *
window push twoway rarea lodds_lb lodds_ub apache,  bcolor(yellow)			///
window push     || line logodds apache 						///
window push     || scatter obs_logodds apache, ytitle("Log Odds of Death")		///
window push     ,  legend(order(2 "Simple log odds" 3 "Observed log odds" 		///
window push             1 "95% CI") rows(1))
window push more
window push *
window push *  Calculate rc_spline covariates with 3 knots
window push *
window push mkspline _Sapache=apache [weight=patients], nknots(3) cubic displayknots
window push *
window push *  Regress deaths against spline covariates
window push *
window push glm deaths _S*, family(binomial patients) link(logit)
window push predict logodds3, xb
window push line logodds apache,  lwidth(thick) color(blue)			///
window push     || line logodds3 apache,  lwidth(thick) color(red)			///
window push     || scatter obs_logodds apache, ytitle("Log Odds of Death")		///
window push     ,  legend(order(1 "Simple log odds" 2 "Spline log odds" 		///
window push             3 "Observed log odds") rows(1))
window push more
window push predict e_deaths, mu
window push generate prob = e_deaths/patients
window push predict se, stdp
window push generate lb_logodds3= logodds3-1.96* se
window push generate ub_logodds3= logodds3+1.96* se
window push generate ub_prob= exp( ub_logodds3)/(1+exp( ub_logodds3))
window push generate lb_prob= exp( lb_logodds3)/(1+exp( lb_logodds3))
window push twoway rarea lb_prob ub_prob apache, bcolor(yellow) 			///
window push     || line prob apache, color(blue) lwidth(medthick)			///
window push     || scatter p apache, ytitle("Probability of Death by 30 Days") 	///
window push     legend(order(2 "Expected mortality rate" 				///
window push            3 "Observed mortality rate" 1 "95% CI") rows(1)) subtitle(	///
window push        "Restricted Cubic Spline Model with 3 Knots",ring(0) position(11))
window push log close