version 10.0 set more on log using 4.18.Sepsis.log, replace * 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) generate ci95lb = 1 - invbinomial(patients, patients-deaths, 0.025) * * Regress deaths against apache * glm deaths apache, family(binomial patients) link(logit) estat vce 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 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) log close