--------------------------------------------------------------------------------- log: C:\Documents and Settings\dupontwd\My Documents\Office\NASUG2005\sup > port.log log type: text opened on: 10 Jul 2005, 21:12:18 . * . * 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 { 2. local log`i' = log(`i') 3. local ylogtick "`ylogtick' `log`i''" 4. } . local yloglabel " " . foreach i of numlist 4(2)10 20(20)100 200 { 2. local yloglabel "`yloglabel' `log`i'' "`i'"" 3. } . gen log_los = log(los) . more . rc_spline meanbp number of knots = 5 value of knot 1 = 47 value of knot 2 = 66 value of knot 3 = 78 value of knot 4 = 106 value of knot 5 = 129 . regress log_los _S* Source | SS df MS Number of obs = 996 -------------+------------------------------ F( 4, 991) = 24.70 Model | 60.9019393 4 15.2254848 Prob > F = 0.0000 Residual | 610.872879 991 .616420665 R-squared = 0.0907 -------------+------------------------------ Adj R-squared = 0.0870 Total | 671.774818 995 .675150571 Root MSE = .78512 ------------------------------------------------------------------------------ log_los | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _Smeanbp1 | .0296009 .0059566 4.97 0.000 .017912 .0412899 _Smeanbp2 | -.3317922 .0496932 -6.68 0.000 -.4293081 -.2342762 _Smeanbp3 | 1.263893 .1942993 6.50 0.000 .8826076 1.645178 _Smeanbp4 | -1.124065 .1890722 -5.95 0.000 -1.495092 -.7530367 _cons | 1.03603 .3250107 3.19 0.001 .3982422 1.673819 ------------------------------------------------------------------------------ . test _Smeanbp2 _Smeanbp3 _Smeanbp4 ( 1) _Smeanbp2 = 0 ( 2) _Smeanbp3 = 0 ( 3) _Smeanbp4 = 0 F( 3, 991) = 30.09 Prob > F = 0.0000 . 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) number of knots = 7 value of knot 1 = 41 value of knot 2 = 60 value of knot 3 = 69 value of knot 4 = 78 value of knot 5 = 101.3251 value of knot 6 = 113 value of knot 7 = 138.075 . regress log_los _S* Source | SS df MS Number of obs = 996 -------------+------------------------------ F( 6, 989) = 16.92 Model | 62.5237582 6 10.4206264 Prob > F = 0.0000 Residual | 609.25106 989 .616027361 R-squared = 0.0931 -------------+------------------------------ Adj R-squared = 0.0876 Total | 671.774818 995 .675150571 Root MSE = .78487 ------------------------------------------------------------------------------ log_los | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _Smeanbp1 | .0389453 .0092924 4.19 0.000 .0207101 .0571804 _Smeanbp2 | -.3778786 .12678 -2.98 0.003 -.6266673 -.12909 _Smeanbp3 | .9316267 .8933099 1.04 0.297 -.8213739 2.684627 _Smeanbp4 | .1269005 1.58931 0.08 0.936 -2.991907 3.245708 _Smeanbp5 | -.7282771 1.034745 -0.70 0.482 -2.758824 1.30227 _Smeanbp6 | -.3479716 .4841835 -0.72 0.473 -1.298117 .6021733 _cons | .6461153 .4496715 1.44 0.151 -.2363046 1.528535 ------------------------------------------------------------------------------ . test _Smeanbp2 _Smeanbp3 _Smeanbp4 _Smeanbp5 _Smeanbp6 ( 1) _Smeanbp2 = 0 ( 2) _Smeanbp3 = 0 ( 3) _Smeanbp4 = 0 ( 4) _Smeanbp5 = 0 ( 5) _Smeanbp6 = 0 F( 5, 989) = 18.59 Prob > F = 0.0000 . 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) number of knots = 7 value of knot 1 = 40 value of knot 2 = 57 value of knot 3 = 74 value of knot 4 = 91 value of knot 5 = 108 value of knot 6 = 125 value of knot 7 = 142 . regress log_los _S* Source | SS df MS Number of obs = 996 -------------+------------------------------ F( 6, 989) = 17.16 Model | 63.3425146 6 10.5570858 Prob > F = 0.0000 Residual | 608.432304 989 .615199498 R-squared = 0.0943 -------------+------------------------------ Adj R-squared = 0.0888 Total | 671.774818 995 .675150571 Root MSE = .78435 ------------------------------------------------------------------------------ log_los | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _Smeanbp1 | .0416447 .0090891 4.58 0.000 .0238086 .0594808 _Smeanbp2 | -.4433587 .1124576 -3.94 0.000 -.6640417 -.2226757 _Smeanbp3 | .9071589 .352071 2.58 0.010 .2162668 1.598051 _Smeanbp4 | -.1009243 .5310043 -0.19 0.849 -1.142949 .9411001 _Smeanbp5 | -.8156462 .6103555 -1.34 0.182 -2.013387 .3820943 _Smeanbp6 | .454242 .5740784 0.79 0.429 -.6723098 1.580794 _cons | .5387057 .4427316 1.22 0.224 -.3300955 1.407507 ------------------------------------------------------------------------------ . 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) number of knots = 7 value of knot 1 = 41 value of knot 2 = 60 value of knot 3 = 69 value of knot 4 = 78 value of knot 5 = 101.3251 value of knot 6 = 113 value of knot 7 = 138.075 . regress log_los _S* Source | SS df MS Number of obs = 996 -------------+------------------------------ F( 6, 989) = 16.92 Model | 62.5237582 6 10.4206264 Prob > F = 0.0000 Residual | 609.25106 989 .616027361 R-squared = 0.0931 -------------+------------------------------ Adj R-squared = 0.0876 Total | 671.774818 995 .675150571 Root MSE = .78487 ------------------------------------------------------------------------------ log_los | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- _Smeanbp1 | .0389453 .0092924 4.19 0.000 .0207101 .0571804 _Smeanbp2 | -.3778786 .12678 -2.98 0.003 -.6266673 -.12909 _Smeanbp3 | .9316267 .8933099 1.04 0.297 -.8213739 2.684627 _Smeanbp4 | .1269005 1.58931 0.08 0.936 -2.991907 3.245708 _Smeanbp5 | -.7282771 1.034745 -0.70 0.482 -2.758824 1.30227 _Smeanbp6 | -.3479716 .4841835 -0.72 0.473 -1.298117 .6021733 _cons | .6461153 .4496715 1.44 0.151 -.2363046 1.528535 ------------------------------------------------------------------------------ . 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 big | Freq. Percent Cum. ------------+----------------------------------- 0 | 949 95.28 95.28 1 | 47 4.72 100.00 ------------+----------------------------------- Total | 996 100.00 . * . * 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 (747 missing values generated) . * . * 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 Logistic regression Number of obs = 996 LR chi2(1) = 29.66 Prob > chi2 = 0.0000 Log likelihood = -545.25721 Pseudo R2 = 0.0265 ------------------------------------------------------------------------------ hospdead | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- meanbp | .9845924 .0028997 -5.27 0.000 .9789254 .9902922 ------------------------------------------------------------------------------ . 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 number of knots = 5 value of knot 1 = 47 value of knot 2 = 66 value of knot 3 = 78 value of knot 4 = 106 value of knot 5 = 129 . logistic hospdead _S*, coef Logistic regression Number of obs = 996 LR chi2(4) = 122.86 Prob > chi2 = 0.0000 Log likelihood = -498.65571 Pseudo R2 = 0.1097 ------------------------------------------------------------------------------ hospdead | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _Smeanbp1 | -.1055538 .0203216 -5.19 0.000 -.1453834 -.0657241 _Smeanbp2 | .1598036 .1716553 0.93 0.352 -.1766345 .4962418 _Smeanbp3 | .0752005 .6737195 0.11 0.911 -1.245265 1.395666 _Smeanbp4 | -.4721096 .6546662 -0.72 0.471 -1.755232 .8110125 _cons | 5.531072 1.10928 4.99 0.000 3.356923 7.705221 ------------------------------------------------------------------------------ . * . * 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 ( 1) _Smeanbp2 = 0 ( 2) _Smeanbp3 = 0 ( 3) _Smeanbp4 = 0 chi2( 3) = 80.69 Prob > chi2 = 0.0000 . * . * 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] +------------------------------------------------------------------+ | _Smean~1 _Smean~2 _Smean~3 _Smean~4 logodds p | |------------------------------------------------------------------| 178. | 60 .32674 0 0 -.7499399 .3208344 | 575. | 90 11.82436 2.055919 .2569899 -2.045913 .114466 | 893. | 120 56.40007 22.30039 10.11355 -1.220147 .2279105 | +------------------------------------------------------------------+ . . lincom (5.531072 + 60*_Smeanbp1 + .32674*_Smeanbp2 /// > + 0*_Smeanbp3 + 0 *_Smeanbp4) /// > /// > -(5.531072 + 90*_Smeanbp1 + 11.82436*_Smeanbp2 /// > + 2.055919*_Smeanbp3 + .2569899*_Smeanbp4) ( 1) - 30 _Smeanbp1 - 11.49762 _Smeanbp2 - 2.055919 _Smeanbp3 - .2569899 _Smea > nbp4 = 0 ------------------------------------------------------------------------------ hospdead | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | 3.65455 1.044734 4.53 0.000 2.086887 6.399835 ------------------------------------------------------------------------------ . . 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) ( 1) 30 _Smeanbp1 + 44.57571 _Smeanbp2 + 20.24447 _Smeanbp3 + 9.85656 _Smeanb > p4 = 0 ------------------------------------------------------------------------------ hospdead | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- (1) | 2.283625 .5871892 3.21 0.001 1.379606 3.780023 ------------------------------------------------------------------------------ . log close log: C:\Documents and Settings\dupontwd\My Documents\Office\NASUG2005\s > upport.log log type: text closed on: 10 Jul 2005, 21:12:44 -------------------------------------------------------------------------------