log using support.log, replace * * Perform linear and logistic regressions of log length of stay and * hospdead against meanbp, respectively. Use restricted cubic splines * for these regessions. * use support.dta, replace set more on * * Analyze relationship between entry mean BP and length of stay * scatter los meanbp /// , msymbol(Oh) xlabel(25(25)175) xtick(30 (5) 170) name(scatter, replace) more * * Define logarithmic tick marks and labels for y-axis of length of stay * local ylogtick = " " foreach i of numlist 3/10 20(10)100 200 { local log`i' = log(`i') local ylogtick "`ylogtick' `log`i''" } local yloglabel " " foreach i of numlist 4(2)10 20(20)100 200 { local yloglabel "`yloglabel' `log`i'' "`i'"" } gen log_los = log(los) more rc_spline meanbp regress log_los _S* test _Smeanbp2 _Smeanbp3 _Smeanbp4 predict y_hat, xb scatter log_los meanbp ,msymbol(Oh) /// || line y_hat meanbp /// , xlabel(25 (25) 175) xtick(30 (5) 170) clcolor(red) /// clwidth(thick) xline(47 66 78 106 129, lcolor(blue)) /// ylabel(`yloglabel', angle(0)) ytick(`ylogtick') /// ytitle("Length of Stay (days)") /// legend(order(1 "Observed" 2 "Expected")) name(knot5, replace) more drop _S* y_hat rc_spline meanbp, nknots(7) regress log_los _S* test _Smeanbp2 _Smeanbp3 _Smeanbp4 _Smeanbp5 _Smeanbp6 predict y_hat, xb scatter log_los meanbp ,msymbol(Oh) /// || line y_hat meanbp /// , xlabel(25 (25) 175) xtick(30 (5) 170) clcolor(red) /// clwidth(thick) xline(41 60 69 78 101 113 138, lcolor(blue)) /// ylabel(`yloglabel', angle(0)) ytick(`ylogtick') /// ytitle("Length of Stay (days)") /// legend(order(1 "Observed" 2 "Expected")) name(knot7, replace) more drop _S* y_hat rc_spline meanbp, nknots(7) knots(40(17)142) regress log_los _S* predict y_hat, xb scatter log_los meanbp ,msymbol(Oh) /// || line y_hat meanbp /// , xlabel(25 (25) 175) xtick(30 (5) 170) clcolor(red) /// clwidth(thick) xline(40(17)142, lcolor(blue)) /// ylabel(`yloglabel', angle(0)) ytick(`ylogtick') /// ytitle("Length of Stay (days)") /// legend(order(1 "Observed" 2 "Expected")) name(evenknots, replace) more drop _S* y_hat rc_spline meanbp, nknots(7) regress log_los _S* predict se, stdp predict y_hat, xb generate lb = y_hat - invttail(_N-7, 0.025)*se generate ub = y_hat + invttail(_N-7, 0.025)*se twoway rarea lb ub meanbp , bcolor(gs6) lwidth(none) /// || scatter log_los meanbp ,msymbol(Oh) mcolor(blue) /// || line y_hat meanbp, xlabel(25 (25) 175) xtick(30 (5) 170) /// clcolor(red) clwidth(thick) ytitle("Length of Stay (days)") /// ylabel(`yloglabel', angle(0)) ytick(`ylogtick') name(ci,replace) /// legend(rows(1) order(2 "Observed" 3 "Expected" 1 "95% CI" )) more predict rstudent, rstudent lowess rstudent meanbp /// , yline(-2 0 2) msymbol(Oh) rlopts(clcolor(green) clwidth(thick)) /// xlabel(25 (25) 175) xtick(30 (5) 170) more gen big = abs(rstudent)>2 tabulate big * * Analyze relationship between hospital death and length of stay * First, calculate rate to equal the proportion of patients with * a given meanbp who die in hospital. * sort meanbp by meanbp: egen rate = mean(hospdead) label variable rate "Hospital Mortality Rate" generate deads = meanbp if hospdead * * Draw an exploratory graph showing the # of patients, the # of deaths * and the mortality rate for each mean BP * twoway histogram meanbp, discrete frequency blwidth(none) gap(20) /// || histogram deads, discrete frequency bcolor(red) blwidth(none) /// gap(20) /// || scatter rate meanbp, yaxis(2) mcolor(blue) msymbol(Oh) /// , ylabel(0 (5) 30, angle(0) axis(1)) xlabel(20 (20) 180) /// xtick(25 (5) 175) ytitle("Number of Patients") /// ytitle(, axis(2) color(blue) ) /// ylabel(0 (.1) 1, angle(0) axis(2) labcolor(blue)) /// ytick(0(.05)1, axis(2) tlcolor(blue)) /// legend(order(1 "Total" 2 "Deaths" 3 "Mortality Rate") rows(1)) more * * Regress hospdead against meanbp using simple logistic regression * logistic hospdead meanbp predict p,p line p meanbp, ylabel(0 (.1) 1) ytitle(Probabilty of Hospital Death) more drop p * * Repeat the preceding model using restricted cubic splines with 5 knots * drop _S* rc_spline meanbp logistic hospdead _S*, coef * * Plot the predicted probability of death against meanbp together with * the 95% confidence interval for this curve and the observed mortality * for each mean bp. * predict p,p predict logodds, xb predict stderr, stdp generate lodds_lb = logodds - 1.96*stderr generate lodds_ub = logodds + 1.96*stderr generate ub_p= exp(lodds_ub)/(1+exp(lodds_ub)) generate lb_p= exp(lodds_lb)/(1+exp(lodds_lb)) twoway rarea lb_p ub_p meanbp, bcolor(gs14) /// || line p meanbp, clcolor(red) clwidth(medthick) /// || scatter rate meanbp, msymbol(Oh) mcolor(blue) /// , ylabel(0 (.1) 1, angle(0)) xlabel(20 (20) 180) /// xtick(25 (5) 175) ytitle(Probabilty of Hospital Death) /// legend(order(3 "Observed Mortality" /// 2 "Expected Mortality" 1 "95% CI") rows(1)) test _Smeanbp2 _Smeanbp3 _Smeanbp4 * * Calculate odds ratios for meanbp = 60 vs meanbp = 90 * meanbp = 120 vs meanbp = 90 * list _S* logodds p /// if (meanbp==60 | meanbp==90 | meanbp==120) & meanbp ~= meanbp[_n-1] lincom (5.531072 + 60*_Smeanbp1 + .32674*_Smeanbp2 /// + 0*_Smeanbp3 + 0 *_Smeanbp4) /// /// -(5.531072 + 90*_Smeanbp1 + 11.82436*_Smeanbp2 /// + 2.055919*_Smeanbp3 + .2569899*_Smeanbp4) lincom (5.531072 + 120*_Smeanbp1 + 56.40007*_Smeanbp2 /// + 22.30039*_Smeanbp3 + 10.11355*_Smeanbp4) /// /// -(5.531072 + 90*_Smeanbp1 + 11.82436*_Smeanbp2 /// + 2.055919*_Smeanbp3 + .2569899*_Smeanbp4) log close