Roland Andersson wants to make a graph similar to figures 3-4 in:
Tekkis, P. P et al. Mortality control charts for comparing
performance of surgical units: validation study using hospital
mortality data. BMJ 2003;326:786 ( 12 April )
He asks:
How can I draw confidence limits for a reference proportion
(proportion of deaths) at increasing number of operations, from 30 to
400? How can I combine such a graph with plots of observed mortality?
He assumes the overall mortality to be 0.15%
--------------------------------------------------------------
1.
Is the average mortality really 0.15% (1 in 667). In that case graphs
are hardly useful since the vast majority of sites will experience no
deaths or a single death. In the following I assume that this was a
typo, the average mortality being 15%.
2.
As I read the paper from Tekkis et al. they calculate the control
intervals as exact binomial confidence intervals as if they had an
observation corresponding to the average mortality, e.g. 15 deaths
out of 100 observations. Tekkis et al. cite Stark et al.
(Stark J et al. Mortality rates after surgery for congenital heart
defects in children and surgeon's performance. Lancet 2000;355:1004-7);
Stark et al describe the procedure. They write:
"This method was used for illustrative purposes only and did not
represent a formal statistical hypothesis test"
I am a bit worried that mistakes may occur, because the intervals
constructed in this way are not confidence intervals proper, but I
accept the figures in Stark et al's paper as illustrative. Tekkis
et al. show both 90%, 95%, and 99% "confidence intervals"; to me
that is overkill, considering the illustrative purpose.
3.
Despite my reservations, the following do-file illustrates how
a graph can be constructed.
HTH
Svend
-------------------------------------------------------------------
clear
set seed 12345
// Creating artificial dataset: 30 hospitals
set obs 30
gen hospital=_n
gen ops = round(300*uniform()) + 30 // Number of operations
expand ops // One observation per operation
gen died = uniform()<0.15 // Mortality 15%
// Back to one observation per hospital; calculate mortality
collapse (sum) died=died (count) operations=died , by(hospital)
gen mortality = died/operations
save mortality.dta , replace
// Create coordinates for control intervals. With a mortality of 15%
// we can use 20,40,60,... operations (N of deaths must be integer)
clear
set obs 20
egen operations=fill(20 40)
gen died= operations*0.15
set level 95
statsby l95=r(lb) u95=r(ub) , by(operations) saving(r95.dta , replace) ///
: cii operations died , binomial
set level 90
statsby l90=r(lb) u90=r(ub) , by(operations) clear ///
: cii operations died , binomial
set level 95
gen mean=0.15
// Combine datasets
append using r95.dta
append using mortality.dta
twoway ///
(line mean l95 u95 l90 u90 operations , ///
lpattern(l dash dash shortdash shortdash)) ///
(scatter mortality operations) ///
, ///
ylabel(0 "0" .1 "10" .2 "20" .3 "30" .4 "40") ///
plotregion(margin(b=0 l=0)) ///
legend(order(1 "Mean" 2 "95% CI" 4 "90% CI")) ///
xtitle("Annual number of operations") ///
ytitle("Mortality (%)")
-------------------------------------------------------------------
________________________________________________________
Svend Juul
Institut for Folkesundhed, Afdeling for Epidemiologi
(Institute of Public Health, Department of Epidemiology)
Vennelyst Boulevard 6
DK-8000 Aarhus C, Denmark
Phone, work: +45 8942 6090
Phone, home: +45 8693 7796
Fax: +45 8613 1580
E-mail: [email protected]
_________________________________________________________
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
