Methods for Hospital

Epidemiology Quality

Improvement

Using R

Basic Statistics

Chapter 3.

1
 Introduction

We describe some statistical methods that we have found

to be useful in infection management (IM) and quality
improvement (QI) departments in hospitals. We find that
staff working in these areas understand proportion and rate
methods that employ differences better than the ratio
estimators that are more familiar in Epidemiology.

An important Evidence Based Medicine estimtor, the

Number Needed to Treat (NNT) and its confidence interval are
derived from the difference between proportions (Sackett,
Strauss, Richardson, Rosenberg, and Haynes 2000). For these
reasons, we describe the calculation of proportions and
rates and their differences together with corresponding
confidence intervals.

Odds ratios are described elsewhere, for example in the

book by Kirkwood and Sterne (2003). Sackett, Strauss,
Richardson, Rosenberg, and Haynes (2000) describe the
relationship between the odds ratio and NNT. Davies, Crombie
and Tavakoli (1998) show how to calculate risk ratios from
odds ratios. Other frequently used statistical methods such
as t-tests, analysis of variance and linear regression are
covered elsewhere, for example in Kirkwood and Sterne
(2003).

Frequently, in Hospital Epidemiology work, samples are

small. The well-known simple methods for calculating
confidence intervals and performing significance tests for
proportions and rates may perform poorly with small samples,
and so we have concentrated on methods which are more
accurate in this situation.

Recently, there has been an increasing interest in the

likelihood approach (Clayton and Hills 1993, Goodman 1999,
van der Tweel 2005). We believe that this has particular
appeal in the areas of QI and IM. Mainstream biostatistics
is predominantly frequentist although Bayesian methods are
being used with increasing frequency. Frequentist inference
is based on the idea that it is possible to think of a large
number of repetitions of an experiment.

In surveillance work in a hospital, there is no large

population from which random samples can be drawn for an
experiment. Often there is an expected value such as an
expected surgical site infection (SSI) rate of, say, 3%.
Surveillance yields data, often the entire sample of
relevant patients in a hospital at a particular time. It is
then natural to use the likelihood approach that is based on

2
the support that data provide for the expected value. If
they provide little support it can be inferred that a
probable difference exists.

We have begun to employ likelihood ideas and we believe

that their use will increase. In particular, the likelihood
counterpart of the confidence interval, the supported range,
has a much simpler and more logical definition than the
frequentist confidence interval (Goodman 1999). However, we
mostly continue to use the more familiar term confidence
interval.

We provide functions in R (R is available at

[Link] that perform the
computer-intensive small sample calculations frequently
required in QI and IM work. R (Ihaka and Gentleman 1996) is
a wonderfully versatile and powerful open source statistical
environment but it is not overly user-friendly for a
hospital scientist who may have little formal statistical
training and who is used to attractive graphical interfaces.
The R functions we provide have been written in a way that
we hope will minimise this difficulty. In some cases, R
commands in the text are shown preceded by “> “. These can
usually be run as a group if the “> “ are removed, for
example by employing Edit Replace and then copying and
pasting at the R prompt.

A graphical user interface is available for R that may

improve user-frendliness for routine statistical testing.
Rcmdr is available at cran and from the author’s Internet
site [Link]

There is excellent R documentation. For example there

are “R for Beginners” by Emmanuel Paradis and “Simple R” by
John Verzani, both available at cran, and “A Quick
Introduction to R” by Lawrence Joseph
([Link] In addition,
excellent introductory books are available such as
Dalgaard’s “Introductory Statistics with R” (2002), “Using R
for Introductory Statistics” by Verzani (2004) and “Data
Analysis and Graphics Using R” by Maindonald and Braun
(2003). Myatt’s “Open Source Solutions–R” (2005) explains
the use of R for epidemiological analyses, including advice
on the use of logistic regression for analysing outbreaks
and epidemics ([Link]). “epi”
([Link] and the R “meta” and
“epitools” packges (cran), have several functions that may
be useful to Hospital Epidemiologists. Other very useful
programs are PEPI 4 ([Link]/[Link])

3
and WINPEPI (brixtonealth), CIA (Altman, Machin, Bryant and
Gardner 2000 and EpiInfo (Dean, Arner, Sangam, Sunki,
Friedman, Lantinga, Zubieta, Sullivan, and Smith 2000.

Getting data into R can be accomplished with the

function [Link](). In all cases a text file can be loaded
and .csv comma delimited text files are convenient for this
purpose. In Microsoft Windows, it is also possible to copy a
file, for example using EDIT COPY in a spreadsheet, and
[Link]() will read the data into R. This does not work in
Linux or Mac OSX. However, in both Windows and Linux it is
possible to use x<-[Link](file=”clipboard”), where x is
the name of a suitable [Link]. This will place the
headings in the data so to use this do not copy the headings
with EDIT COPY. Also, make sure that numerical data are
specified as such, for example using FORMAT CELLS NUMBER
before copying or they may not copy correctly.

Gibberd, Pathmeswaran and Burtenshaw (2000) have

advocated the use of shrikage estimators for the analysis of
QI data. We discuss this further when we describe the use of
funnelplots. Two further issues of importance are
denominators and risk adjustment. These are dealt with in
more detail in Chapter 2.

In Section 1 of this chapter we describe proportions

and in Section 2 counts and rates. Section 3 is devoted to
ordered categorical data, Section 4 to numerical data, and
Section 5 to a brief outline of the miscellaneous topics of
decision making and rater agreement.

Section 1
Binary data.
Introduction

In this section we deal with binary data. These data

are often binomial, although risk adjusted binary data are
often not binomially distributed. First we describe methods
to calculate confidence intervals and perform a significance
test for a single proportion. Next we deal with data in two
or more independent proportions, and matched data. Following
this we examine methods for stratified data. There is then a
brief section dealing with graphical methods.

Example.

4
 A consecutive series of 106 patients undergoing upper
abdominal surgery was investigated to determine the presence
of postoperative respiratory complications (McGrath and
Morton 1986). The aim was to determine which patients might
benefit from referral for additional preoperative
respiratory preparation. The following variables were
examined - age, history of smoking, diagnosis of chronic
bronchitis or emphysema on the patient’s file, and peak
expiratory flow rate (PEFR) which can be measured in the
surgical outpatients department using simple equipment.
These patients are part of a larger group of 232 patients
that will be used to illustrate logistic regression in
Chapter 5. Initially, univariate analyses were conducted as
follows -

1. Age. There were 47 patients whose age was below 60 years

and the number of complications was 5, a rate of 11%. 59
patients were 60 or above and the number of complications in
this group was 22, a rate of 37%.

2. Smoking history. There were 65 patients who gave a

history of having smoked less than 20 pack years with 9
complications, a rate of 14%. 41 patients smoked 20 or more
pack years with 18 complications, a rate of 44%.

3. Bronchitis. In 21 cases there was a diagnosis of chronic

bronchitis, emphysema or chronic obstructive pulmonary
disease on the patient’s file and 12 of these patients had a
postoperative respiratory complication, a rate of 57%. This
diagnosis was absent in 85 cases and 15 of these patients
had complications, a rate of 18%.

4. PEFR. Finally, there were 59 patients who had a PEFR of

at least 75% of their predicted value with 9 complications,
a rate of 15%. 47 patients had PEFR of less than 75% with 18
complications giving a rate of 38%.

Using R –

> load("[Link]") # data from clipboard (Windows) or text file

> [Link]() # [Link] has 232 records; first 106 required
Loading data.
Data from clipboard (C) or file (F) c # copy and paste from clipboard
Do data column(s) have heading(s) (Y/N) y
> datain[1,] # first record
PatientNumber Age PackYrsSmkd COPD PEFR Outcome
1 1 17 0 0 93 0
> age<-datain[,1];outcome<-datain[,5] # age and outcome columns
> age1<-age[1:106];outcome1<-outcome[1:106] # selects first 106 records
> age2<-age1 # to tabulate ages
> age0<-cut(age2,c(15,30,45,60,75,90)) # ages in categories
> x<-table(age0,outcome1) # tabulates ages and outcomes

5
> x # many more complications (outcome1=1) from 60 on
outcome1
age0 0 1
(15,30] 7 0
(30,45] 22 2
(45,60] 16 4
(60,75] 30 17
(75,90] 4 4
> age1[age1<60]<-0 # selects ages below 60
> age1[age1>=60]<-1 # ages 60 and above
> k<-[Link](table(outcome1,age1)) # tabulates data
> k
outcome1 age1 Freq
1 0 0 42
2 1 0 5
3 0 1 37
4 1 1 22
>
Note that data entry from the clipboard is only available for Windows
versions of R. It is preferable to use a text file if in doubt.
Comma-delimited (.CSV) files work well in all cases.

Proportion data, single rate.

When confronted with these or similar data we wish to

make them as intelligible as possible. An estimate of each
proportion, commonly called a rate, has already been
calculated (P=X/N=5/47=0.11 or 11% for the younger age
group). Here, X is the number of positive outcomes (for the
lower age group data there were 5 complications) and N is
the sample size (there were 47 patients in the lower age
group). In addition, a measure of the precision of each
estimate is required, its confidence interval. In likelihood
terms, this tells us the range of possible values for the
estimate that are supported by the data according to some
critical value of the likelihood ratio (Clayton and Hills
1993). For example, the likelihood ratio of 1/7
approximately coincides with a 95% supported range that is
similar to a 95% confidence interval.

When samples are small, as they often are in QI and IM

studies, there is a great deal of uncertainty about whether
the calculated rate is the true rate for the patients being
studied, in this case patients undergoing upper abdominal
surgery. For example, if there were 4 patients and 2 had
complications it would be difficult to believe that these 4
patients with their 50% complication rate were
representative of all such patients. If the next 2 patients
studied did not have complications, the rate would fall from
50% to 30%. Using the likelihood approach, the supported
range (confidence interval) provides a range of values that
are supported by the data.

6
Confidence interval.

An accurate method is to use the inverse beta function;

in R the required commands for a 95% confidence interval are
qbeta(0.025,X,N-X+1) and qbeta(0.975,X+1,N-X). qbeta(0.025,5,43)=.035
or 3.5% and qbeta(0.975,6,43)=0.231 or 23.1%. We prefer mid-p
values as the above results may be slightly conservative due
to the discrete nature of the binomial distribution. To
obtain mid-p values, also select qbeta(0.025,X+1,N-X) and
qbeta(0.975,X,N-X+1). These give 0.0483 and 0.2038. The average
of 0.0355 and 0.0483 is 0.0419 (4.2%) and that of 0.2310 and
0.2038 is 0.2174 (21.7%); these are the mid-p values.

A frequent requirement is to compare an observed rate

with a reference rate. For example, one hospital’s clean
(Class 1) SSI rate for a certain procedure may be compared
with a rate obtained for a group of hospitals that is
regarded as a reference rate. Clearly a rate slightly above
the reference may be due to random variation and a
significance test is required to determine whether the
difference between the observed and reference rates is too
large for random variation to be likely to explain it. To
illustrate the method, suppose a rate of 15% was expected
for postoperative respiratory complications. We wish to
determine whether the overall observed 25% rate (X=27,
N=106) is so large that random variation is unlikely to have
produced it.

Significance test

An accurate significance test result may be obtained by

using pbeta, the beta distribution function in R. If the
observed rate exceeds the expected rate, as it does with the
postoperative respiratory function data (100×27/106=25%
observed, 15% expected), pbeta(P,X,N-X+1) is employed in R,
where P is the expected 15% rate (pbeta(0.15,27,80)=0.0034).
This is a one tailed p-value and it needs to be multiplied
by 2 to give the two-sided p-value of 0.0068. To get the
mid-p value, pbeta(P,X+1,N-X) is also employed; this gives a
one tailed p-value of 0.0016. Now add this to the one tailed
p-value of 0.0034; the result is 0.005. When the observed
proportion is less than expected, 1-pbeta(P,X+1,N-X) is used
and the result must be multiplied by 2 for a two tailed
result. To obtain the mid-p value, 1-pbeta(P,X,N-X+1) is also
employed and the result is added to the one tailed p-value
obtained with 1-pbeta(P,X+1,N-X). If X is zero, an upper mid-p
value cannot be obtained.

Likelihood ratio (Bayes factor).

7
 We have referred in the introduction to this chapter to
the usefulness of the likelihood ratio (LR). An approximate
LR can be obtained by converting the p-value to a standard
normal deviate (Z) using Z=abs(qnorm(p-value/2)) and then
employing paste("1/",[Link](round(1/exp(-z^2/2),0)),sep="").
For the mid-p value of 0.005, the LR is 1/51. This suggests
that there is 51 times the support for the observed rate of
25% as there is for the expected rate of 15%, making random
variation an unlikely explanation for the observed 25% rate
(an LR of 1/7 is equivalent approximately to a p-value of
0.05). In addition, the confidence interval for P=27/106=25%
suggests that the values supported by the data at the 95%
level range from 18% (above the expected value of 15%) to
34%.

Using R -

The above calculations are available in the R function

proportion(). It is suggested that the R functions that
accompany this book be placed in a suitable directory, for
example it could be called e:/rfunctions (We like to keep
these functions on a separate USB disk that in our system is
the e: drive. Note that R uses / rather than \ for referring
to disks and directories as it uses \ for formatting
output). Note also that in R case matters; R is not the same
as r. When R has been opened, the directory needs to be
changed, for example to e:/rfunctions using File and Change
dir in the R menu. Then at the R prompt enter
load(“[Link]”). To run the function, enter
proportion(); do not omit the (). Now enter the numerator
(X) and the denominator (N). Next, enter y or Y if an
expected proportion is available and then the expected
proportion value. The commands and output are –

> load(“[Link]”)
> proportion()
Enter numerator 27
Enter denominator 106
Is a reference proportion available? y
Enter reference proportion .15
Proportion = 0.255, lower 95% limit = 0.175, upper limit = 0.349.
Mid-P 95% limits are 0.179 and 0.343.
p-value = 0.007, LR = 1/39, mid-p value = 0.005, LR = 1/51.

Difference between two proportions.

In IM and QI studies one frequently compares two

independent proportions and the odds ratio is often used for
this purpose (Kirkwood and Sterne 2003). However, NNT is
likely to be a better estimator for most IM and QI work;

8
this is the reciprocal of the rate difference. Thus, if a
complication rate is 20% and a new process is able to reduce
it to 10%, the difference would be 10%. By using the new
method for the next 10 patients, one complication could be
averted - NNT is therefore 10. When the proportion becomes
larger following an intervention, it is referred to as the
Number Needed to Harm (NNH) (Sackett, Strauss, Richardson,
Rosenberg, and Haynes 2000).

Calculation of the confidence interval for NNT employs

the difference between proportions. It is described by
Altman in Appendix 1 of Sackett, Strauss, Richardson,
Rosenberg, and Haynes (2000).

Confidence interval.

Newcombe (1998) has proposed an accurate method that is

relatively simple to use. First, one must obtain the
confidence limits for each proportion, for example by using
the methods described above. Then the formulas are –
U=D+√[(P2-L2)2+(U1-P1)2] and L=D-√[(P1-L1)2+(U2-P2)2], where P1
and P2 are the two proportions, L1, L2 and U1, U2 are the
lower and upper confidence limits for the two proportions, D
is the difference between the 2 proportions and U and L are
the required upper and lower confidence limits for D.

Ratio of proportions.

Calculation of the ratio of two proportions and its

confidence interval are often of lesser importance in QI and
IM work. Nam (1995) describes an accurate score method. It
involves the solution of a complicated cubic equation that
we include in the R function as it does not appear to have
been implemented in the packages to which we have referred.
A difficulty occurs if there is a zero numerator; when this
occurs 0.5 can be added to each of the X’s and 1 to each of
the N’s. Nam’s formula is shown below.

Confidence interval calculation.

The method involves the solution of the following

equation –

a1P3+a2P2+a3P+a4=0
P P

where
a1=N1{N1(N1+N2)X2+N2(N1+X2)Z2},
a2=-N1{N1N2(X1+X2)+2(N1+N2)X1X2+N2(N1+X1+2X2)Z2},
a3=2N1N2X1(X1+X2)+(N1+N2)X12X2+N1N2(X1+X2)Z2 and
a4=-N2X12(X1+X2).

9
 To solve this equation define b1=a2/a1, b2=a3/a1,
b3=a4/a1, c1=b2-b12/3 and c2=b3-b1b2/3+2b13/27. The equation is
then t3+c1t+c2=0. The 3 roots of the equation are
t=-2(-c1/3)^.5cos(π/3-θ/3), t=-2(-c1/3)0.5cos(π/3+θ/3) and
t=2(-c1/3)0.5cos(θ/3), where cosθ=270.5c2/{2c1(-c1)0.5}. The 3
values of P are then P=t-b1/3. The largest of these is
discarded and then RRU,L=(1-A/B)/P, where A=(N2-X2)/(1-P) and
B=X1+N2-(N1+N2)×P.

Example.

A vascular surgical unit found that in 74 consecutive

Class 1 operations there were 14 SSIs, a rate of 19%. Since
this rate was considered unacceptably high, the processes of
wound care were carefully revised and in the following 114
operations there were 8 infections giving a rate of 7%.

The difference in the rates is 12% so the NNT is

approximately 8. The 95% confidence interval calculated
using Newcombe’s formula is 2.3% to 22.6% (NNT 4 to 43). The
ratio of the proportions or risk ratio is 2.7 and the 95%
confidence interval is 1.2 to 6.

Significance test.

Gart and Nam (1990) describe a variation of the

familiar chi-squared test that incorporates a skewness
correction and we use it in twoproportions(), the R function
for the difference and ratio of two proportions. Its formula
is Z=(P1-P2)/√[P×(1-P)×(1/N1+1/N2)], where P1=X1/N1, P2=X2/N2
and P=(X1+X2)/(N1+N2). Let S=X1+X2 and T=N1+N2, A=S/T, B=1-S/T,
C=1-2×S/T, D=S×(T-S)/T2 and E=(1/N1+1/N2). The skewness
correction is then γ=√[(A×B×C/N1-A×B×C/N2)/(D×E)]. The skew
corrected Z score value is ZS=Z-γ×(ZS2-1)/6. If N1=N2 then
ZS=Z, otherwise ZS is the solution of
ZS={-1+√(1+A×B)]}/(2×γ/6), where A=4×γ/6 and B=Z+4×γ/6. With
the above data, Z=2.42 and p-value=0.016. In addition, the
likelihood ratio of 1/19 indicates that the support for a
zero difference is only 1/19 times the support for the
observed difference of 12%.

Using R -

> load(“[Link]”)
> twoproportions()
Enter first numerator 14
Enter first denominator 74
Enter second numerator 8

10
Enter second denominator 114
Difference between proportions 0.119.
Lower 95% limit 0.023, upper limit 0.226.
Z = 2.42, P = 0.016, LR = 1/19.
Fisher Exact p-value = 0.019.
Ratio 2.7.
Ratio 95% confidence limits are 1.22 and 6.

A Fisher Exact test result is also reported and may be

more accurate than the chi-squared test when samples are
very small. This test may be conservative and a mid-p value
may be preferred.

Using R proceed as follows to get the mid-p value –

x<-c(14,8);y<-c(74-14,114-8);d<-[Link](x,y);f<-[Link](d)
> f

Fisher's Exact Test for Count Data

data: d
p-value = 0.01909
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.126773 8.972299
sample estimates:
odds ratio
3.072157

> x1<-c(15,7);y1<-c(74-15,114-7);d1<-[Link](x1,y1)
> f1<-[Link](d1)
> # add one to the numerator of the larger proportion
> # and subtract one from the numerator of the smaller proportion
> # will not work if latter is 0.
> f1

Fisher's Exact Test for Count Data

data: d1
p-value = 0.004739
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.386706 11.845307
sample estimates:
odds ratio
3.856681

> pvalue<-(f$[Link]+f1$[Link])/2
> pvalue
[1] 0.01191275

These data and their analysis suggest that the fall in

the rate of surgical site infections following the
intervention is unlikely to have been due to random
variation. The NNT of 8 suggests that for every 8 operations
performed there was one less surgical site infection than
before.

11
More than 2 independent proportions.

In some cases, there may be more than 2 proportions to

study and sometimes it may make sense to expect that a trend
might exist among the proportions. These data may be
analysed with the R function manyproportions().
Occasionally, these data may be stratified and, although
stratification methods for their analysis exist, logistic
regression will usually be preferred. Logistic regression is
described in Chapter 5.

Example.

The following data were obtained from the study of the

causes of postoperative respiratory complications among
patients undergoing abdominal surgery described above.

PEFR below or equal to 50% of predicted - 16 complications

among 21 patients.
PEFR above 50% and below or equal to 75% - 22 complications
among 73 patients.
PEFR above 75% and below or equal to 100% - 27 complications
among 123 patients.
PEFR above 100% - 2 complications among 15 patients.

Using R –

> load("[Link]")
> manyproportions()
Enter the number of proportions to compare 4

Successes Totals Proportions Expecteds Residuals

1 16 21 0.76 6.06 5.02
2 22 73 0.30 21.08 0.29
3 27 123 0.22 35.52 -2.47
4 2 15 0.13 4.33 -1.37

4-sample test for equality of proportions without continuity correction.

Chisq = 27.581, DF = 3, P-value = 0.
p-value by Monte Carlo simulation = 0.
Is a trend expected (Y/N) y
Chi-squared Test for Trend in Proportions.
Chisq = 20.192, DF = 1, P-value = 0.
Departure from trend chisq = 7.389, DF = 2, P-value = 0.025.
Warning message:
Chi-squared approximation may be incorrect in: [Link](Successes,
Totals)

There is a warning message as one of the expected

values is less than 5; however, it is over 4 and is the only
one below 5 so the chi-squared result is unlikely to be
seriously in error. A less conservative rule is that less

12
than 20% of the expected cell frequencies should be below
five and none should be below one. The Monte Carlo
simulation test result agrees with this. This test may be
run separately.

Using R -

> x<-c(16,22,27,2)
> y<-c(5,51,96,13)
# y is n-x (n is the vector 21, 73, 123 and 15)
> xy<-[Link](x,y)
> [Link](xy,[Link]=T,B=10000)

Pearson's Chi-squared test with simulated p-value (based on 10000

replicates)

data: xy
X-squared = 27.5813, df = NA, p-value = 1e-04.

The Fisher Exact test is also an alternative to the

global chi-squared test when samples are small.

Using R –

> x<-c(16,22,27,2)
> y<-c(5,51,96,13)
> xy<-[Link](x,y)
> [Link](xy)

Fisher's Exact Test for Count Data

data: xy
p-value = 1.242e-05
alternative hypothesis: [Link]

The test manyproportions() shows a strong trend in

complications as PEFR diminishes. However, there is also
some departure from trend. The patients with the lowest PEFR
values had many more complications than expected.

Funnel plots for binary data from several institutions.

Currently, the funnel plot is the preferred method for

the display of binary data from multiple institutions. These
data are increasingly risk-adjusted. It is simply a Shewhart
type control chart like those described in Chapter 4 that
has the institutions on the horizontal axis instead of
times. These are sorted from smallest to largest to give the
funnel effect as the smaller institutions have wider control
limits. When the data are risk-adjusted, it is usual to
employ indirect standardisation and the chart is then
described as being multiplicative. It has recently been
suggested that the mid-line of the chart can be the weighted

13
average of the risk-adjusted rates (Hart, Lee, Hart and
Robertson 2003) and this is used in the R function
groupfunnelb(). Comparisons with the standard rate may
sometimes be required. However, when they differ, the risk-
adjustment probabilities should probably be re-calibrated
for the group of institutions being studied. This is easily
accomplished using logistic regression described in Chapter
5.

For each hospital the observed (O) and expected (E)

outcome numbers are calculated and an SMR (O/E) is
determined. This is then multiplied by an overall average
expected value to obtain a risk-adjusted rate that is then
plotted for the hospital concerned. Exact binomial control
limits for the weighted average of the risk-adjusted rates
are calculated using the pbeta function in R. Control limits
and not confidence limits are used and this requires the
implementation of a simple numerical search. If πi is the
expected outcome probability for the ith subject, these
limits are often calculated from the variance obtained when
the πi×(1-πi) values for each institution are added.
However, this gives limits that are too narrow. The
preferred method is the subject of some controversy but the
method describd above gives useful approximations.

If Oi is the outcome for the ith subject, Ei is the

expected outcome probability for that subject and E is the
overall average expected outcome probability, the variable
for institution J is AJ=(∑JOi/∑JEi)×E} and the center line is
either E or P=∑(nJ×AJ)/∑nJ, where nJ is the sample size for
institution J. The 2 sigma equivalent control limits are
calculated for institution J by finding the values of xJ
such that qnorm(pbeta(P,xJ+1,nJ-xJ))=2 for the upper limit and
qnorm(pbeta(P,xJ,nJ-xJ+1))=-2 for the lower limit. These must be
found by trial and error or by implementing a simple
numerical search that finds the values of xJ on either side
of ±2 and then employs linear interpolation.

A limitation of this chart is that smaller institutions

may be shown in an unfavourable light due to random
variation from effects related to regression to the mean
that can be important for these institutions. To deal with
this problem, shrinkage methods are advocated (Gibberd,
Pathmeswaran and Burtenshaw 2000); here the institutions are
assumed collectively to have a distribution of their own.
Unfortunately, in the process of improving precision,
shrinkage can increase bias. The result is that, while
better control of random variation can reduce the likelihood
of a false positive result, shrinkage may suggest that a

14
performance is satisfactory when it is not; the possibility
of false negative error is thus increased. One approach
would be to perform the analysis both with and without
shrinkage. If a result without shrinkage suggests the
presence of a problem, the staff in the institution of
interest should study the relevant systems in that
institution. If no cause for the aberrant result can be
found and the shrunken value suggests acceptable
performance, confidence that the result is due to random
variation may be increased.

Shrinkage, when there is no risk-adjustment, is

straightforward to perform in the software program BUGS
(Woodworth 2004). When risk-adjustment is employed, a more
complex multi-level model is required. This can involve
fitting an hierarchical logistic regression requiring
specialised skills. However, since we are chiefly interested
in large risk-adjusted rates in smaller institutions,
employment of a simpler model analogous to the method for
the analysis of SMRs for multiple institutions may be
feasible. If this approach is valid, analysis of the
log(Odds) of the risk-adjusted rates can be employed as with
data that do not require risk-adjustment, when these rates
are too large for an SMR analysis to be feasibe. The use of
BUGS is beyond the scope of this work and hospital
scientists wishing to employ shrinkage should consult with a
statistician familiar with the BUGS software.

These ideas are summarised in the following charts.

Table 1 shows SSIs for a single operation performed in 10
hospitals. The corresponding Funnel Plot is Figure 1. This
suggests that Hospital 1 has excess SSIs (8/48=17%, where
the average for all the hospitals is just over 5%). Figure 2
shows the shrinkage-adjusted values and the 95% credible
interval for Hospital 1 now includes the average value. It
is almost certain that Hospital 1 has excessive SSIs and its
systems of wound care should be subjected to inquiry by its
staff. In the unlikely event that no system problem could be
found, it is possible that the result may be due to random
variation.

The Funnel Plot described above that employs indirect

standardisation is described as a multiplicative plot. Hart,
Lee, Hart and Robertson (2003) also describe an additive
plot for the O and E differences. However, the
multiplicative plot in various forms is in wide use so we do
not describe the additive plot. It is desirable to use the
minimum number of plots that analyse the data clearly as the
use of many charts increases the likelihood of false
positive results.

15
 Funnel plots can be used for screening institutions to
detect those requiring more rigorous study, for example with
charts like CUSUMs and O-E charts or within institution
funnel plots that are described in Chapter 4. It should be
noted that it is possible to have a run or cluster of
unsatisfactory outcomes that require investigation and that
do not show up in funnel plots, especially if outcomes are
otherwise good. The longer the period of observation, the
greater the likelihood of this occurring. It is imperative
that suitable tables and charts that display monthly or
quarterly data accompany funnel plots. These can be
accomplished with spreadsheet tabulations or in a
statistical program like R and are described in Chapter 4.

It is possible to have a result above the control limit

due to unsatisfactory outcomes that the funnel plot
identifies as an outlier early in the series and to have
more recent outcomes that are satisfactory. Funnel plots
should be based on the most recent data possible.
Unfortunately, with some outcomes, many hospitals can have
insufficient volume of relevant subjects for funnel plots to
be based only on the most recent data. It is important to
train institution staff to do their own within institution
charting, for example with CUSUM and O-E or within
institution funnel plots. Then institutions could be subject
to disciplinary action if (1) charts were not maintained and
reported, (2) signals occurred that were not investigated or
(3) systems problems were revealed by an investigation but
not dealt with. These charts tell us there may be a problem;
they do not tell us there definitely is a problem nor do
they give any indication of the cause of a possible problem.
These questions can only be answered by undertaking an
analysis of the underlying system using quality improvement
and industrial psychology methods described in Chapter 1.

Overdispersion can be an issue with Funnel plots for

very large data sets like those involving many hospitals.
Spiegelhalter (2005) has described a method to deal with
this problem that employs an adaptation of the
DerSimonian-Laird procedure for Meta-analysis. This is
beyond the scope of this work.

Using R –

> load("[Link]") # To tabulate SSI data for 10 hospitals

> [Link]()
Loading data. # File [Link], Column 1 hospital number
Data from clipboard (C) or file (F) c # Column 2 [0,1] SSI outcome
Do data column(s) have heading(s) (Y/N) y

16
> q<-table(datain[,2],datain[,1]) # Tabulates data
> hospnum<-[Link](names(q[1,])) # Extracts hospital numbers
> ssi<-[Link](q[2,]) # SSIs
> ops<-[Link](q[1,]) # Not SSIs
> ops<-ops+ssi # Total operations
> prop<-round(ssi/ops,3) # Proportion SSIs to 3 decimal places
> dataout<-[Link](hospnum,ssi,ops,prop) # Data in data frame
> dataout
hospnum ssi ops prop
1 1 8 48 0.167
2 2 4 74 0.054
3 3 2 52 0.038
4 5 3 69 0.043
5 6 1 58 0.017
6 7 3 58 0.052
7 8 6 151 0.040
8 9 1 16 0.062
9 10 3 57 0.053
10 12 4 43 0.093
> [Link](dataout,file="g:/[Link]",sep=",",[Link]=F)
# Saves tabulated data

> load("[Link]") # Funnelplot of SSI data for 10 hospitals

> groupfunnelb() # File [Link]
Do you have risk-adjusted data (Y/N) n
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
The mean outcome rate is 0.056
Do you wish to enter another value (Y/N) n
hosp outc tot exp
1 1 8 48 2.6837061
2 2 4 74 4.1373802
3 3 2 52 2.9073482
4 5 3 69 3.8578275
5 6 1 58 3.2428115
6 7 3 58 3.2428115
7 8 6 151 8.4424920
8 9 1 16 0.8945687
9 10 3 57 3.1869010
10 12 4 43 2.4041534
Chart heading.
Enter heading for chart eg Hospital/Surgeon Mortality SSI rates for 10
hospitals
X-axis heading.
Enter an X-axis heading eg Hospital/Surgeon Number Hospital number
Y-axis heading.
Enter an Y-axis heading eg Risk-adjusted rate SSI rate # Figure 1.

> load("[Link]") # Chart for shrunken SSI data in Figure 2.

> [Link]() # Data from [Link]; analysis in BUGS
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> datain
Hospital Median l2.5 u97.5
1 9 0.05657 0.02317 0.11360
2 12 0.06262 0.03389 0.11950
3 1 0.07749 0.04433 0.16330
4 3 0.05233 0.02195 0.09148

17
5 10 0.05525 0.02628 0.09629
6 6 0.04708 0.01442 0.08155
7 7 0.05508 0.02593 0.09660
8 5 0.05278 0.02378 0.09083
9 2 0.05546 0.02808 0.09495
10 8 0.04942 0.02518 0.07775
> x<-c(1,2,3,4,5,6,7,8,9,10) # Indexes hospitals for “arrows” to
> datain<-cbind(datain,x) # produce credible intervals
> datain
Hospital Median l2.5 u97.5 x
1 9 0.05657 0.02317 0.11360 1
2 12 0.06262 0.03389 0.11950 2
3 1 0.07749 0.04433 0.16330 3
4 3 0.05233 0.02195 0.09148 4
5 10 0.05525 0.02628 0.09629 5
6 6 0.04708 0.01442 0.08155 6
7 7 0.05508 0.02593 0.09660 7
8 5 0.05278 0.02378 0.09083 8
9 2 0.05546 0.02808 0.09495 9
10 8 0.04942 0.02518 0.07775 10
> d<-datain
> da<-d[,1]
> par(lab=c(length(da),5,7));par(xaxs="r")
> plot(d[,2],axes=F,ylim=c(min(d[,3]),max(d[,4])),xlab="Hospital
number",ylab="",main="Shrunken SSIs,\nHospital rates blue, 95%
credible limits red, mean rate black.",col="blue",lwd=3)
> box()
> axis(side=1,tick=T,labels=[Link](da));axis(side=c(2,3,4))
> arrows(d[,5],d[,3],d[,5],d[,4],angle=90,code=3,col="red")
> abline(h=.056)

Matched data analysis.

Matched data analysis is required infrequently with IM

and QI data. In the R function matched(), we deal only with
the case where there is 1:1 matching. A program like PEPI
(Abramson and Gahlinger 1999) or WINPEPI can be employed
when there is more than 1 control per case. Data may be
entered as the numbers of concordant and discordant pairs or
from a .csv file in 3 columns with column A containing the
case-control group (1,1;2,2;3,3;..), column B 1 for cases
and 0 for controls (0;1;0;1;0;1;..) and column C a zero or a
1 if the attribute of interest is absent or present –

A B C
1 0 0
1 1 0
2 0 0
2 1 1

We illustrate the function matched() using the first

data entry method. In a recent study the following was
found.

1. Number of concordant pairs with positive outcomes 7

18
2. Number of discordant pairs with cases positive 9
3. Number of discordant pairs with cases negative 3
4. Number of concordant pairs with negative outcomes 28

Using R –

> load(“[Link]”)
> matched()
Enter data.
Enter data from keyboard (K) or clipboard/file (F) k

Enter number of concordant pairs with positive outcomes 7

Enter number of discordant pairs with cases positive 9
Enter number of discordant pairs with cases negative 3
Enter number of concordant pairs with negative outcomes 28

Odds ratio 3.
Lower 95% limit 0.894, upper limit 11.978.
LR = 1/4, z = 1.683, p-value (mid) = 0.092.

Cases proportion 0.34, controls proportion 0.213.

Correlated proportion difference 0.128.
Lower 95% limit -0.02, upper limit 0.269.

The function employs the binomial test for a single

proportion to the discordant pairs described at the
beginning of this section (proportion()). The expected
proportion is 0.5. The correlated proportions and their
differences are analysed using the method described by
Newcombe (Altman, Machin, Bryant and Gardner 2000).

Stratified proportion data, single group rates.

Frequently it is necessary to stratify data because

there are subgroups of observations that are expected to
have different outcome rates. For example, SSIs may occur in
a clean (Class 1) surgical wound or one that is contaminated
by bacteria (Class 2), or already infected (Class 3).
Clearly, contaminated or already infected wounds are more
likely to develop an SSI than clean ones. A hospital may
appear to have a high infection rate if it performs an
unusually large number of operations on patients in the
higher risk groups; this is commonly called confounding
(Rothman 1986). It is necessary to stratify the data by
these variables to deal with this potential confounding
effect. They can then be analysed by the methods of direct
and indirect standardisation (Kirkwood and Sterne 2003).

There are other factors that increase the risk of an

SSI such as the length of the operation and the ASA class
that are described in Chapter 2. However, the data we use to
illustrate the methods to deal with this problem were

19
collected some years ago and are divided into clean and
contaminated groups only, there being very few infected
wounds in elective vascular surgery. Length of operation and
ASA class data were not collected for these operations.

Standardisation.

Standardisation may be direct or indirect. In indirect

standardisation the observed proportion of infections is
compared to an expected proportion obtained by calculating
the number of infections that would have occurred if some
reference infection rate had been applied to the observed
group of patients.

Example (continued).

Suppose that an SSI is expected to occur to 5% of

patients undergoing lower limb vascular procedures who have
clean wounds and 10% who have contaminated wounds (Ei=1..2).
The observed rates (Oi) are 14/74=19% for clean wounds and
3/15=20% for contaminated wounds, giving an overall rate of
17/89=19%. Although the rates are sufficiently alike for
standardisation to be unnecessary, we use these data to
illustrate the methods of indirect and direct
standardisation.

Indirect standardisation.

With 74 clean wounds and an expected rate of 5% there

would have been 3.7 infections and for the 15 contaminated
wounds with an expected infection rate of 10% there would
have been 1.5 infections. Thus in the 89 operations there
would have been 5.2 infections giving an expected rate (E)
of 5.8%.

In most cases, where the expected rates are less than

about 10%, one can use analysis methods for an SMR (O/E)
based on the Poisson distribution described in Section 2 of
this chapter. The R function rate() can be used for this
purpose. However, ICU and CCU mortality rates may be too
large to use these methods.

When the expected rates exceed about 10%, one can

calculate confidence limits for the overall observed
proportion in this case 17/89=19.1% and compare them with
the overall expected rate, here 5.2/89=5.8% using
proportion(). However, it is important to realise that this
is only an approximation as the combined rate will not be
binomially distributed. It is probably wise to avoid the

20
less conservative mid-p intervals. An alternative is to use
direct standardisation described in the following section.

Using R –

> load("[Link]")
> proportion()
Enter numerator 17
Enter denominator 89
Is a reference proportion available? y
Enter reference proportion .058
Proportion = 0.191, lower 95% limit = 0.115, upper limit = 0.288.
Mid-P 95% limits are 0.12 and 0.282.
P = 0, LR = 1/7006, mid-p value = 0, LR = 1/11033.

> load("[Link]")
> rate()
Enter numerator 17
Is there a denominator (y/n) y
Enter the denominator 89
Is a reference rate available? y
Enter reference rate .058
Rate = 0.19101, lower 95% limit = 0.11127, upper limit = 0.30583
Mid-P 95% limits are 0.11557 and 0.29889.
P = 0, LR = 1/3198.
Mid-p value = 0, LR = 1/4878
# Although the rate is above 10%, the result is still similar to that
obtained with proportion(). However, the latter should be used with
these data.

Direct standardisation.

This requires knowledge of the proportions of the

patients expected to be in each stratum in a reference
population such as a large group of hospitals, in this case
for clean (Class 1) and contaminated (Class 2) wounds for
lower limb vascular surgery. Suppose, for example, they are
known to be 0.9 and 0.1. The formula for the directly
standardised rate D is then D=Σ(Wi×Xi/Ni), where Wi are
weights representing the expected proportions of patients in
each of the strata (in this case 0.9 and 0.1). Also, Xi and
Ni are the number of surgical site infections and operations
in stratum i. The Wi must be scaled to sum to one, although
this scaling is performed automatically in the direct
standardisation R function stratifiedproportions(). In the
present case i=1..2 and X1=14 and X2=3, N1 is 74 and N2 is
15. The value of D with these data is
D=0.9*14/74+0.1*3/15=0.19 or 19%.

Confidence interval.

The method relies on re-scaling (Waller, Addy, Jackson

and Garrison 1994). For each stratum (i=1..I) calculate PiU,L

21
using the mid-p method. Then XU,L=ΣWiPiU,L, P=ΣWiPi, a re-
scaling factor R=(ΣWi2/N1)0.5/Σ(Wi/Ni0.5), and PU=P+R×(XU-P) and
PL=R×(P-XL). Newcombe uses an alternative method that
employs the square-and-add procedure he has described
([Link]/medicine/epidemiology_statistics/
research/statistics/[Link]). Some informal
comparisons suggest that the two methods give similar
results. Using the re-scaling method, the confidence
interval is 12% to 29%. If the expected rates were known to
be 5% and 10% for the two strata, the expected rate could be
calculated by Σ(Wi×Ei/Ni)=.9*.05+.1*.1=.055 and this could be
compared with the confidence limits; the lower limit of 12%
is much higher.

Using R -

> load(“[Link]”)
> stratifiedproportions()
Enter number of strata 2
ENTER to continue
Enter numerator in stratum 1 14
Enter denominator in stratum 1 74
Enter weight for stratum 1 .9
Press ENTER for the next stratum
Enter numerator in stratum 2 3
Enter denominator in stratum 2 15
Enter weight for stratum 2 .1
ENTER to continue
Do you wish to review the data (y/n)? y
Weighted average = 0.19.
Lower 95% limit = 0.122, upper limit = 0.285.

If you wish to review the data, for example because of

a possible data entry error, select y or Y when asked if you
wish to review the data. R’s simple spreadsheet appears with
the data that may be changed if necessary; when they are
satisfactory click on the X in the top right corner of the
spreadsheet to return to the R screen.

With these data the three approaches give very similar

confidence limits with the Poisson based limits being
slightly wider. However, the rates in the two strata are
very similar. Suppose that instead there were 15 SSIs in 150
class 1 operations and 3 SSIs in 15 class 2 procedures. The
confidence intervals, using .9 and .1 as the weights for the
direct standardisation, are shown below. Once again the
results for the three methods are similar.

Using R –

> stratifiedproportions()
Enter number of strata 2

22
ENTER to continue
Enter numerator in stratum 1 15
Enter denominator in stratum 1 150
Enter weight for stratum 1 .9
Press ENTER for the next stratum
Enter numerator in stratum 2 3
Enter denominator in stratum 2 15
Enter weight for stratum 2 .1
ENTER to continue
Do you wish to review the data (y/n)? n
Weighted average = 0.11.
Lower 95% limit = 0.071, upper limit = 0.168.
> proportion()
Enter numerator 18
Enter denominator 165
Is a reference proportion available? n
Proportion = 0.109, lower 95% limit = 0.066, upper limit = 0.167.
Mid-P 95% limits are 0.068 and 0.163.
> rate()
Enter numerator 18
Is there a denominator (y/n) y
Enter the denominator 165
Is a reference rate available? n
Rate = 0.10909, lower 95% limit = 0.06465, upper limit = 0.17241
Mid-P 95% limits are 0.06699 and 0.16869.

Stratified proportion data, differences between rates.

Indirect standardisation.

An approximate test to compare two indirectly

standardised rates has been described for SSIs (Horan and
Culver 1996). Since most expected adverse outcome rates are
less than 10%, the method employs the Poisson distribution,
and it is described in Section 2 of this chapter. It may not
be suitable for use with ICU and CCU mortality rates that
exceed 10%. However, the Mantel-Haenszel (MH) method
described below can be used to compare two groups of
stratified proportion data. It is possible to compare the
indirectly standardised rates, as described for control
charts by Hart, Lee, Hart and Robertson (2003), but these
rates are not binomially distributed so the method is only
approximate. Table 2 shows some stratified SSI data. For the
Class 1 operations a 5% rate can be expected and for the
Class 2 procedures the expected rate is 10%. The number of
operations in group 1 is 89 and the expected number of SSIs
is 74*.05+15*.1=5.2. For group 2 the expected number of SSIs
is 114*.05+36*.1=9.3. Thus in 89+150=239 procedures 14.5 SSI
would be expected so that the average expected SSI rate is
14.5/239=.06067. There were 17 SSIs in group 1 and 11 in
group 2. The indirectly standardised rate for group 1 is
then .06067*17/5.2=.1983. For the second group it is
.0607*11/9.3=.0718.

23
Using R –

> 89*(14.5/239)*17/5.2
[1] 17.65248 # Adjusted “number” in first group
> 150*(14.5/239)*11/9.3
[1] 10.76394 # Adjusted “number” in second group.
> x<-c(17.65248,10.76394)
> n<-c(89,150)
> [Link](x,n)
2-sample test for equality of proportions with continuity
correction
data: x out of n
X-squared = 7.3779, df = 1, p-value = 0.006603
alternative hypothesis: [Link]
95 percent confidence interval:
0.02506366 0.22810209
sample estimates:
prop 1 prop 2
0.1983425 0.0717596.

This result is very similar to the Mantel-Haenszel

result shown below (p-value = 0.006) and the rate difference
result from Section 2 of this chapter (p-value = 0.009).

Using R –

> 5.2/14.5 # Expected proportion of SSIs in first group

[1] 0.3586207
> [Link](17,28,.3586207) # 17 SSIs in first group
Exact binomial test # 28 total SSIs
data: 17 and 28
number of successes = 17, number of trials = 28, p-value = 0.009231
alternative hypothesis: true probability of success is not equal to
0.3586207
95 percent confidence interval:
0.4057682 0.7849572
sample estimates:
probability of success
0.6071429

Direct standardisation.

Rates may be compared using the directly standardised

difference. The Mantel-Haenszel (MH) method (Kirkwood and
Sterne 2003, Greenland and Robins 1985) is usually suitable
for performing the analysis although strictly speaking it
does not employ direct standardisation. Weights chosen for
each stratum are proportional to 1/(1/N1i+1/N2i) where N1i and
N2i are the numbers in each group in stratum i; the weights
must be scaled so that their sum is one. The sum of the
reciprocals of N1i and N2i is proportional to the variance of
their difference. Weights that are inversely proportional to
the variance are generally desirable. The formula for
calculating the directly standardised difference is

24
D=ΣWi(X1i/N1i-X2i/N2i). The MH procedure is implemented in the
R function mhci().

Recall that, for the clean surgery group, the before

and after review results were 14/74 and 8/114 respectively
and for the contaminated surgery group the before result was
3/15. The number of contaminated operations following the
process change was 36 and 3 of the wounds became infected.
These data are shown in Table 2. For these data the MH
difference is 12%.

Confidence interval.

Risk Difference.

The MH risk difference formula including a correction

for small numbers in the strata (Greenland and Robins 1985,
Sato 1989) is shown below. The 95% confidence interval
(Z=1.96) is 2.7% to 21%.

DU,L=D±Z×√[D×(∑Pi+∑Qi)/A], where
Pi={N1i2X2i-N2i2X1i+N1iN2i[(N2i-N1i)/2]}/(N1i+N2i)2,
Qi={[X1i(N2i-X2i)/(N1i+N2i)]+[X2i(N1i-X1i)/(N1i+N2i)]}/2 and
A=[∑N1i×N2i/(N1i+N2i)]2.

Risk ratio.

The MH risk ratio is

RR=Σ[X1iN2i/(N1i+N2i)]/Σ[X2iN1i/(N1i+N2i)]. Its confidence limits
are RRU,L=exp(loge(RR)±Z×[A/(B×C)]0.5), where
A=∑{[(X1i+X2i)N1iN2i-X1iX2i(N1i+N2i)]/(N1i+N2i)2},
B=∑[X1iN2i/(N1i+N2i)] and C=∑[X2iN1i/(N1i+N2i)]. For these data
the risk ratio is 2.63 with a 95% confidence interval from
1.28 to 5.39.

Significance test.

The formula for the large sample significance test is

Z=(ΣWiX1i/N1i-ΣWiX2i/N2i)/√V and χ21=Z2, where
V=Σ[WiPiQiNi/(NI-1)], Ni=N1i+N2i, Pi=(X1i+X2i)/(N1i+N2i) and
Qi=(1-Pi) (Miettinen 1985). For these data Z=2.73, χ21=7.48,
p-val=0.006.

From this analysis we conclude that the 12% fall in the

overall infection rate following the intervention is
unlikely to have been due to random variation. The NNT is 8
approximately so that the change noted for the clean (Class

25
1) surgery also occurred when both clean (Class 1) and
contaminated (Class 2) operations were combined.

When comparing stratified rates it is important to

ensure that the rate differences in each stratum are
approximately similar. If they differ markedly, it is
inappropriate to use a summary measure such as the
Mantel-Haenszel risk difference. For example, if the risk
difference in one stratum is positive and in another
negative, they would be likely to cancel out thus suggesting
that there was no difference when in fact there were
opposing differences in each stratum. This is called
interaction or, when it has a definite biological basis,
effect modification (Kirkwood and Sterne 2003).

There are formal tests called homogeneity tests for

determining whether the data within the strata are suitable
for amalgamation but they can be inefficient, especially
when the data are sparse. In addition, different tests are
required for difference and ratio estimators. We do not
include homogeneity tests. These are available in other
software, for example WINPEPI. However, the user should
examine the data and if different strata suggest markedly
differing results, within stratum analyses only should be
performed and the overall tests and confidence intervals
ignored. For the data in the above example the rate
difference in each stratum was 12%, so that there is no
evidence of inhomogeniety. When a summary measure is
required in the presence of inhomogeniety, for example for a
meta-analysis, the DerSimonian-Laird procedure may be
employed (Kirkwood and Sterne 2003). Meta-analysis is beyond
the scope of this work. However, the R meta library or
WINPEPI may be used.

Using R-

> load("[Link]")
> mhci()
Enter number of strata 2
ENTER to continue
Enter first numerator in stratum 1 14
Enter first denominator in stratum 1 74
Enter second numerator in stratum 1 8
Enter second denominator in stratum 1 114
Press ENTER for the next stratum
Enter first numerator in stratum 2 3
Enter first denominator in stratum 2 15
Enter second numerator in stratum 2 3
Enter second denominator in stratum 2 36
ENTER to continue
Do you wish to review the data (y/n)? n
Stratified risk ratio = 2.631.

26
Lower 95% limit = 1.284, upper limit = 5.39.
Stratified risk difference = 0.119.
Lower 95% limit = 0.027, upper limit = 0.21.
Chisq = 7.478, LR = 1/42, P = 0.006.

If you wish to review the data, for example because of

a possible data entry error, select y or Y when asked if you
wish to review the data. R’s simple spreadsheet appears with
the data that may be changed if necessary; when they are
satisfactory click on the X in the top right corner of the
spreadsheet to return to the R screen.

Very small samples.

The MH method is suitable for large samples. There may

be many strata with sparse data in each stratum or few
strata with large large counts in each stratum. However, in
IM and QI work in hospitals there are frequently few strata
with sparse data in each. The MH procedures over-estimate
significance and can give confidence intervals that are
inaccurate with these data.

Newcombe ([Link]/medicine/epidemiology
_statistics/research/statistics/[Link]) has
described a method for weighted average differences in
stratified data that employs the square-and-add procedure he
has described. The R function [Link] implements the method
using MH weights. The approximate p-value is calculated from
the value of Z employed when the difference confidence
interval just includes zero and the risk ratio confidence
interval is test based (Miettinen 1985). The test based
method is only valid when the p-value is not very small and
with these data that will usually be the case. Table 3 shows
some hypothetical data that illustrate the method.

Using R –

> load("[Link]")
> nci()
Enter number of strata 3
ENTER to continue
Enter first numerator in stratum 1 1
Enter first denominator in stratum 1 32
Enter second numerator in stratum 1 0
Enter second denominator in stratum 1 97
Press ENTER for the next stratum
Enter first numerator in stratum 2 4
Enter first denominator in stratum 2 43
Enter second numerator in stratum 2 3
Enter second denominator in stratum 2 142
Press ENTER for the next stratum
Enter first numerator in stratum 3 2
Enter first denominator in stratum 3 21

27
Enter second numerator in stratum 3 2
Enter second denominator in stratum 3 49
ENTER to continue
Do you wish to review the data (y/n)? n

Newcombe's method using exact binomial 95% confidence limits

MH weighted average in first group 0.073

First group confidence limits 0.036 to 0.158
MH weighted average in second group 0.018
Second group confidence limits 0.007 to 0.048
MH weighted difference 0.055
Difference confidence limits 0.007 to 0.14
MH weighted ratio 4.025
Approximate 95% ratio confidence limits 1.27 to 12.762
Approximate P=0.018, approximate LR=1/16

> d1<-.073 # Confidence limits using Newcombe’s method

> d2<-.018
> u1<-.158
> l1<-.036
> u2<-.048
> l2<-.007
> dl<-(d1-d2)-((d1-l1)^2+(u2-d2)^2)^.5
> dl
[1] 0.007
> du<-(d1-d2)+((d2-l2)^2+(u1-d1)^2)^.5
> du
[1] 0.14

> load("[Link]")
> mhci()
Enter number of strata 3
ENTER to continue
Enter first numerator in stratum 1 1
Enter first denominator in stratum 1 32
Enter second numerator in stratum 1 0
Enter second denominator in stratum 1 97
Press ENTER for the next stratum
Enter first numerator in stratum 2 4
Enter first denominator in stratum 2 43
Enter second numerator in stratum 2 3
Enter second denominator in stratum 2 142
Press ENTER for the next stratum
Enter first numerator in stratum 3 2
Enter first denominator in stratum 3 21
Enter second numerator in stratum 3 2
Enter second denominator in stratum 3 49
ENTER to continue
Do you wish to review the data (y/n)? n
Stratified risk ratio = 4.025.
Lower 95% limit = 1.333, upper limit = 12.154.
Stratified risk difference = 0.055.
Lower 95% limit = 0.001, upper limit = 0.109.
Chisq = 7.008, LR = 1/33, P = 0.008.

The MH test gives a p-value of .008. The Fisher exact

mid-p value implemented in WINPEPI is .02 for these data.
The MH procedure should be employed when the data not very

28
sparse. However, when there are few strata and few data
values within them, the Newcombe difference confidence
interval may be more reliable.

Using Graphics.

Graphical display of the results of these analyses is

often valuable. Figure 3 is an illustration of proportions
of patients admitted during 1994 and 1995 with heart failure
or shock (ANDRG 252) whose lengths of stay exceeded 5 days
in two hospitals that have been compared. Five days was the
median length of stay (LOS) in a large number of hospitals
that included the two being compared. These data are
described in more detail in Section 4 of this chapter. Also
the difference between the proportions for the two hospitals
is shown. In each case 95% confidence intervals have been
included in the graphical display. The proportions for the
two hospitals were 59% (95% confidence interval 53% to 65%)
and 43% (95% confidence interval 35% to 51%) for Hospital 1
and Hospital 2 respectively, and the difference was 16% (95%
confidence interval 6% to 25%).

The magnitude of the difference between the two

hospitals is clearly apparent in Figure 3 and, since the
lower confidence limit for the difference exceeds zero, it
is unlikely that it was due to random variation. However,
this does not necessarily mean that one hospital is more
efficient than the other as the characteristics of the
patients in the two institutions or the complexity of the
treatment available in each institution may differ. Further
study would be required to determine whether the populations
of the two hospitals were sufficiently alike for the
difference to be attributed to differing performance of the
institutions. Nevertheless, the value of a graphical display
for presenting the results of these investigations is clear.

Using R -

Figure 3 was produced by the following R commands –

> h<-c(.59,.43,.16)
> u<-c(.65,.51,.25)
> l<-c(.53,.35,.06)
> x<-c(1,2,3)
> k<-c("Hosp 1","","Hosp 2","","Diff")
> plot(h,ylim=c(0,.75),axes=F,main="Proportion exceeding reference data
median for ANDRG 252 and\ndifference between hospitals 1 and 2 with 95%
confidence intervals.",xlab="",ylab="")
> box()
> axis(side=1,tick=F,labels=k)
> axis(side=2)
> arrows(x,u,x,l,angle=90,code=3,col="red")

29
 In this section we have described methods for analysing
data in the form of proportions as follows -

1. Methods for single proportions

2. Methods for two or more independent proportions
3. The 1:1 Matched data method
3. Methods for stratified data
4. Graphical methods.

In the next section we deal with count and rate data

methods.

Section 2
Count and rate data
Introduction

We deal with count and rate data in this section. These

data occur frequently in QI and IM work. For example, most
nosocomial infections, excluding postoperative surgical site
infections, can be analysed as monthly counts of events.
Other examples include such adverse occurrences as patient
falls, medication errors, readmissions, pressure ulcers and
needlestick injuries.

When analysing nosocomial infection data, the number of

patients and objects that are spreading the organism,
commonly called the colonisation pressure or multiple
antibiotic resistant organism (MRO) burden, is very
important (Merrer, Santoli, Vecchi, Tran, Jonghe and Outin
2000). Other important factors include hygiene measures,
particularly hand washing, antibiotic usage, prevalence in
the community surrounding the hospital and the adequacy of
its isolation facilities and discharge planning.

First, we begin with methods for single samples,

followed by methods for two independent samples. Next we
describe methods for stratified data and a function for
analysing data in a large contingency table. Following this,
we discuss the common problem of breakdown on independence
and the likely ensuing increase in variation. We then refer
briefly to the use of graphical methods for the analysis of
these data; this aspect is dealt with more fully in the
following chapter that is devoted to control charts.

30
Statistical Analysis of rate and count data.

The statistical methods for count data analysis are

similar to those for proportions and it is sometimes
possible to use methods developed for analysing proportion
data on count and rate data. To do this the denominator of
the rate must be multiplied by a large constant so that its
variance is reduced to near zero (Rothman and Boice 1982).
Following the analysis, the result must then be multiplied
by the constant to return to the original scale.

When proportions are less than 10% and simple

approximate methods are being used, it is frequently
preferable to employ count and rate data methods rather than
methods for proportions.

Example.

During a 3-month surveillance period there were 12

unplanned readmissions for 32,190 occupied bed-days in a
small teaching hospital and during the same months in the
succeeding year there were 28 similar unplanned readmissions
for 37,440 occupied bed-days.

The number of readmissions was greater in the second

period but the hospital’s bed occupancy was also higher.
Also the numbers of readmissions could have been greater in
the second period because the surveillance may have been
more efficient, as discussed in Chapter 2. However, for the
present purpose, we shall use these data as if ascertainment
were similar in the two periods. The rates for the two
periods were respectively 0.37 and 0.75 per 1000 occupied
bed-days.

Certain patients are much more prone to requiring

readmission than others and some may need repeated
admissions. This may lead to increased variability and may
cause the methods of this section that are based on the
Poisson distribution to give misleading results. It is
therefore advisable whenever possible to check that the mean
and variance of these and similar count data are alike
before using the methods of this section.

During an 11 month period of relatively constant bed

occupancy, the following monthly counts of unplanned
readmissions were observed - 4, 5, 5, 2, 9, 6, 4, 6, 3, 7 &
7. The mean of these counts is 5.3 and the variance is 4.
Although one would hope to base these calculations on more
data, for example 20 to 24 monthly counts, there is no
evidence of excessive variation.

31
Single count or rate.

An estimate of each rate has already been calculated.

If X is the number of events and N is the person-time of
exposure, then the observed rate is R=X/N. A measure of the
precision of the rate, its confidence interval, is also
required. In likelihood terms, this tells us the range of
possible values of the rate that are supported by the data.

To obtain an accurate confidence interval, we employ

the inverse gamma function. For the lower 95% limit
[qgamma(0.025,X,1)]/N and for the upper value
[qgamma(0.975,X+1,1)]/N are employed in R. To get mid-p values
change X to X+1 for the lower limit and average this with
the result for X. For the upper limit substitute X for X+1
and average the two values. The resulting intervals are 0.19
to 0.65 and 0.20 to 0.63 for the mid-p value. An upper mid-p
value cannot be obtained if X=0.

Significance test.

It is often desirable to compare the observed rate with

a reference rate. For example, the unplanned readmission
rate for a surveillance period may be compared with a rate,
obtained for a large group of hospitals that is regarded as
a reference rate.

To illustrate the method, suppose that the reference

rate is 0.3 per 1000 occupied bed-days and it is to be
compared with the rate of 28 readmissions in 37,440 occupied
bed-days that occurred in the second surveillance period.
For these data, the expected count is E=3×37440/10000=11.232

We employ the gamma distribution function in R. When

the expected count E is greater than X, employ 2×
[1-pgamma(E,X+1,1)] and when X is greater than E, use
2×pgamma(E,X,1) to obtain two tailed p-values. With these data
the two tailed p-value is 0.000036. The value of Z, the
standard normal deviate, can be obtained using the qnorm
function and qnorm(0.0000362/2)=4.13. For mid-p values,
substitute X for X+1 when E is greater than X and X+1 for X
when X is greater than E and average the resulting 2 values.
With these data, the p-value is 0.0000138 when X+1 is
substituted for X and the average of 0.0000362 and 0.0000138
is 0.0000250. z=qnorm(0.0000250/2)=4.21.

In addition, the likelihood ratio is calculated using

LR<-exp(-z^2/2); it is 1 to 7202 for the mid-p value. This

32
means that the support from the data for the expected rate
of 0.3 per 1000 is very small.

Using R -

> load(“[Link]”)
> rate()
Enter numerator 28
Is there a denominator (y/n) y
Enter the denominator 37440
Is a reference rate available? y
Enter reference rate .0003
Rate = 0.00075, lower 95% limit = 5e-04, upper limit = 0.00108.
Mid-P 95% limits are 0.00051 and 0.00107.
P = 0, LR = 1/5061.
Mid-p value = 0, LR = 1/7202.

Count data, two independent rates.

The need to compare two independent rates arises

frequently in IM and QI studies. Usually this is
accomplished by using the rate ratio as an estimator of the
magnitude of the difference between the two rates (Kirkwood
and Sterne 2003). In a previous section we described the
number needed to treat (NNT). Although we are now dealing
with count and rate data rather than proportion data, we
believe that the same approach is worth following. The
inverse of the rate difference can be thought of as the
person-time needed to treat to prevent one complication.

Example.

If there have been 6 cases of ventilator associated

pneumonia in 1000 ventilator-days and a new method is able
to reduce this to 3, then the rate difference is 3 per 1000
ventilator-days. One case can be prevented each 330
ventilator-days. However, as we have remarked earlier, even
if this estimator is not used, it has been our experience
that staff working in IM and QI departments tend to have a
better intuitive understanding of the magnitude of
difference estimators than they do of ratio estimators. For
these reasons we prefer the rate difference to the rate
ratio when dealing with count data, although we include both
estimators and their confidence intervals.

Rate difference confidence interval.

Newcombe’s method for proportions can easily be adapted

for count and rate data and it is relatively simple to use.
First, one must obtain the confidence limits for each rate,
for example by using the method described above. Then the

33
formulas are DU=D+√{(R2-L2)2+(U1-R1)2} and
DL=D-√{(R1-L1)2+(U2-R2)2} where L1, L2 and U1, U2 are the lower
and upper confidence limits for the two rates (R1 and R2)
and D, DL and DU are the rate difference and its lower and
upper confidence limits respectively.

Rate ratio (RR) confidence interval.

The score method is employed (Graham, Mengersen and

Morton 2003). The formulas are –

Let A=2X1X2, B=Z2(X1+X2), C=√{Z2(X1+X2)[4X1X2+Z2(X1+X2)]} and

D=2X22. Then L=(N1/N2)×(A+B-C)/D and U=(N1/N2)×(A+B+C)/D.

Example (continued).

Comparing readmission rates of 12 in 32,190 and 28 in

37,440 occupied bed-days, the rate difference is 0.38 per
1000 bed-days with a 95% confidence interval of 0.02 to 0.73
per 1000 bed-days. The rate ratio is 2.01 and its 95%
confidence interval is 1.03 to 3.9.

Significance test.

A very conservative method is to regard X1/(X1+X2) as a

proportion with expected value N1/(N1+N2) and use the methods
for a single proportion as described above. For the above
data the p-value=0.055. Due to the conservative nature of
this method, mid-p values are preferred; the mid-p value is
0.039.

Using R -

> load(“[Link]”)
> tworates()
Are denominators available? y
Enter first numerator 28
Enter first denominator 37440
Enter second numerator 12
Enter second denominator 32190
Rate ratio = 2.01.
Lower 95% limit = 1.03, upper limit = 3.9.
Difference between rates = 0.00038.
Lower 95% limit = 2e-05, upper limit = 0.00073,
p-value = 0.03941, LR = 1/8.

Funnel plots for count data from several institutions.

34
 Funnel plots have been described in Section 1 of this
chapter. Here we illustrate their use with count data.
Annual bacteraemia data for a group of hospitals are shown
in Table 4. There are 3 hospital levels with Level 1 being
the large referral institutions. These data are analysed
with the R function groupfunnelc(). Figure 4 is the
funnelplot for the combined data in [Link]; there
are many outliers and the weighted variance for these data
is much larger than the mean. Level 1 and Level 2 & 3
hospitals need to be examined separately (Figure 5 & Figure
6). There are no high outliers in Figure 5
([Link]). The low outlier represents an obstetric
hospital. Figure 6 shows one high outlier ([Link]).
This is a Level 2 hospital that should review its system of
patient care; for example, if these bacteraemias were
related to intravenous devices, it is possible that they may
be left in place for too long. However, the explanation may
be that this hospital does work typical of a Level 1
institution, for example it may have a specialist oncology-
haematology or renal unit.

Using R –

> load(“[Link]”)
> groupfunnelc()
Do you have risk-adjusted data (Y/N) n
Loading data. # Data for level 1 hospitals
Data from clipboard (C) or file (F) f
Have you selected the correct directory (Y/N) y
Enter the name of the file [Link]
Do data column(s) have heading(s) (Y/N) y
Are data columns separated by spaces (S) or commas (C) c
The mean count is 53.833.
Do you want to change the mean value (Y/N) n
The mean denominator is 61701.
Do you want to change the mean denominator value (Y/N) n
Denominators.
Enter name for denominator (eg Bed-days) Bed-days
Denominators.
Per thousand (1000), hundred (100), ten (10), unit (1) Bed-days 1000
Chart heading.
Enter heading for chart eg Bacteraemias by hospital Bacteraemias level 1
hospitals # Figure 5

> groupfunnelc()
Do you have risk-adjusted data (Y/N) n
Loading data. # Data for level 2 & 3 hospitals
Data from clipboard (C) or file (F) f
Have you selected the correct directory (Y/N) y
Enter the name of the file [Link]
Do data column(s) have heading(s) (Y/N) y
Are data columns separated by spaces (S) or commas (C) c
The mean count is 9.333.
Do you want to change the mean value (Y/N) n
The mean denominator is 30438.

35
Do you want to change the mean denominator value (Y/N) n
Denominators.
Enter name for denominator (eg Bed-days) Bed-days
Denominators.
Per thousand (1000), hundred (100), ten (10), unit (1) Bed-days 1000
Chart heading.
Enter heading for chart eg Bacteraemias by hospital Bacteraemias level 1
hospitals # Figure 6

Chi-squared and trend tests for count data from several

institutions.

We illustrate the chi-squared test for comparing rates

using the data for the Level 2 & 3 hospitals.

> load("[Link]")
> [Link]()
Loading data.
Data from clipboard (C) or file (F) f
Have you selected the correct directory (Y/N) y
Enter the name of the file [Link]
Do data column(s) have heading(s) (Y/N) y
Are data columns separated by spaces (S) or commas (C) c
> datain # Data for level 2 & 3 hospitals
HospitalID Bacteraemias [Link]
1 1 23 40587
2 2 22 76995
3 3 3 34303
4 18 8 33722
5 19 7 17016
6 6 8 40829
7 23 14 29390
8 9 7 21893
9 12 12 31713
10 15 5 19716
11 17 2 10815
12 20 1 8271

> x<-sum(datain[,2])
> y<-sum(datain[,3])
> x # Total bacteraemias
[1] 112
> y # Total bed-days
[1] 365250
> z<-x/y # Average rate
> z
[1] 0.0003066393
> Exp<-z*datain[,3] # Expected bacteraemias for hospitals
> dd<-cbind(datain,Exp) # Adds expecteds to [Link]
> dd
HospitalID Bacteraemias [Link] Exp
1 1 23 40587 12.445569
2 2 22 76995 23.609692
3 3 3 34303 10.518648
4 18 8 33722 10.340490
5 19 7 17016 5.217774
6 6 8 40829 12.519775
7 23 14 29390 9.012129

36
8 9 7 21893 6.713254
9 12 12 31713 9.724452
10 15 5 19716 6.045700
11 17 2 10815 3.316304
12 20 1 8271 2.536214
# More than expected for hospital 1 and less for hospital 3 as shown in
Figure 6
> ch<-(datain[,2]-a)/a^.5 # Calculates residuals
> ch
[1] 2.9917649 -0.3312818 -2.3182474 -0.7278405 0.7802254 -1.2773755
1.6615046 0.1106703 0.7297154 -0.4252887 -0.7228181 -0.9646254
# Large residuals for hospitals 1 & 3
> chsq<-sum((datain[,2]-a)^2/a) # Chi-squared
> chsq
[1] 22.14404
> def<-length(ch)-1 # Degrees of freedom
> def
[1] 11
> 1-pchisq(chsq,def) # p-value
[1] 0.02328082
# Note 2 expecteds below 5. However, less that 1 in 5 below 5 and none
below 1 so chi-squared result should be reliable.

In addition, a trend test may be indicated to see if

the number of bacteraemias is related to the size of the
hospital. The easiest way to do this is to employ the built
in R function [Link](). However, because we are
dealing with count data, it is necessary to multiply the
denominators by a large constant, Rothman and Boice (1982)
recommend 1000000. Suppose we believed that the bacteraemia
rates were related to the size of the Level 2 & 3 hospitals.
The test shows that there is no evidence of a trend. When a
trend is present, the departure-from-trend should be
assessed as described in Section 1 of this chapter.

Using R –

> datain # From previous example

HospitalID Bacteraemias [Link]
1 1 23 40587
2 2 22 76995
3 3 3 34303
4 18 8 33722
5 19 7 17016
6 6 8 40829
7 23 14 29390
8 9 7 21893
9 12 12 31713
10 15 5 19716
11 17 2 10815
12 20 1 8271
> o<-order(datain[,3])
> o # Order from smallest to largest
[1] 12 11 5 10 8 7 9 4 3 1 6 2
> d<-datain[o,]
> d
HospitalID Bacteraemias [Link]

37
12 20 1 8271
11 17 2 10815
5 19 7 17016
10 15 5 19716
8 9 7 21893
7 23 14 29390
9 12 12 31713
4 18 8 33722
3 3 3 34303
1 1 23 40587
6 6 8 40829
2 2 22 76995
> [Link](datain[,2],datain[,3]*1000000)

Chi-squared Test for Trend in Proportions

data: a out of b * 1e+06 ,

using scores: 1 2 3 4 5 6 7 8 9 10 11 12
X-squared = 0.028, df = 1, p-value = 0.8671.

Note that these data frequently display overdispersion

and the above methods then can give misleading results.
Poisson regression is a safer alternative; when there is
evidence of overdispersion the standard errors and
confidence intervals can be corrected by using quasipoisson
in the glm() function. Poisson regression is described in
Chapter 5.

Using R –

> d1<-datain[,2];d2<-datain[,3]
> x<-1:length(d1)
> y<-glm(d1~x,poisson,offset=log(d2))
> summary(y)

Call:
glm(formula = d1 ~ x, family = poisson, offset = log(d2))

Deviance Residuals:
Min 1Q Median 3Q Max
-2.8780 -0.8375 -0.3280 0.8192 2.2660

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.96271 0.17383 -45.807 <2e-16 ***
x -0.02584 0.03043 -0.849 0.396
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 22.640 on 11 degrees of freedom

Residual deviance: 21.909 on 10 degrees of freedom
AIC: 71.229

Number of Fisher Scoring iterations: 4

# The residual deviance is about double the degrees of freedom

38
> y1<-glm(d1~x,quasipoisson,offset=log(d2))
> summary(y1)

Call:
glm(formula = d1 ~ x, family = quasipoisson, offset = log(d2))

Deviance Residuals:
Min 1Q Median 3Q Max
-2.8780 -0.8375 -0.3280 0.8192 2.2660

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.96271 0.25100 -31.724 2.28e-11 ***
x -0.02584 0.04394 -0.588 0.57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 2.084995)

Null deviance: 22.640 on 11 degrees of freedom

Residual deviance: 21.909 on 10 degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

We repeat of the above analysis using the count data

trend test {chisq=[sum(score×(o-e))]^2/[sum(score^2×e)-
(sum(score×e))^2/sum(o)]}.

Using R –

> [Link]()
Loading data. # Data for level 2 and level 3 hospitals
Data from clipboard (C) or file (F) f
Have you selected the correct directory (Y/N) y
Enter the name of the file g:/examples/[Link]
Do data column(s) have heading(s) (Y/N) y
Are data columns separated by spaces (S) or commas (C) c
> datain
HospitalID Bacteraemias [Link]
1 1 23 40587
2 2 22 76995
3 3 3 34303
4 18 8 33722
5 19 7 17016
6 6 8 40829
7 23 14 29390
8 9 7 21893
9 12 12 31713
10 15 5 19716
11 17 2 10815
12 20 1 8271
> a<-1:length(datain[,1]) # Scores for trend
> o<-datain[,2] # Observed bacteraemias
> py<-datain[,3] # Bed days
> e<-py*sum(o)/sum(py) # Expected bacteraemias
> ch<-sum((o-e)^2/e) # global chi-squared test
> ch
[1] 22.14404

39
> df<-length(o)-1 # degrees of freedom
> 1-pchisq(ch,df) # p-value
[1] 0.02328082
> res<-(o-e)/e^.5 # residuals, hospitals 1 & 3 are outliers
> res
[1] 2.9917649 -0.3312818 -2.3182474 -0.7278405 0.7802254 -1.2773755
[7] 1.6615046 0.1106703 0.7297154 -0.4252887 -0.7228181 -0.9646254
> tr<-(sum(a*(o-e)))^2/(sum(e*a^2)-(sum(e*a))^2/sum(o))
> tr # trend test
[1] 0.7217855
> 1-pchisq(tr,1) # p-value
[1] 0.3955589

> gch<-glm(o~a,poisson,offset=log(py)) # Poisson regression

> summary(gch)

Call:
glm(formula = o ~ a, family = poisson, offset = log(py))

Deviance Residuals:
Min 1Q Median 3Q Max
-2.8780 -0.8375 -0.3280 0.8192 2.2660

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.96271 0.17383 -45.807 <2e-16 ***
a -0.02584 0.03043 -0.849 0.396 # same as trend test
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 22.640 on 11 degrees of freedom # similar to global

# chi-squared
Residual deviance: 21.909 on 10 degrees of freedom
AIC: 71.229

Number of Fisher Scoring iterations: 4

> gch1<-glm(o~a,quasipoisson,offset=log(py)) # Regression repeated

# to deal with overdispersion (residual deviance 21.91 on 10 degrees
# of freedom
> summary(gch1)
Call:
glm(formula = o ~ a, family = quasipoisson, offset = log(py))

Deviance Residuals:
Min 1Q Median 3Q Max
-2.8780 -0.8375 -0.3280 0.8192 2.2660

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.96271 0.25100 -31.724 2.28e-11 ***
a -0.02584 0.04394 -0.588 0.57
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 2.084995)

Null deviance: 22.640 on 11 degrees of freedom

40
Residual deviance: 21.909 on 10 degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

Stratified count data, single group rate.

Example.

Surgical site infection (SSI) data for two surveillance

periods have been stratified by NNIS risk index (Table 3).
Although these are proportion data rather than counts or
rates, each proportion is less than 10% so the methods of
this section are appropriate for their analysis.

Stratified analysis for count data rates is similar to

that for proportion data using standardisation that may be
direct or indirect.

Indirect standardisation.

In indirect standardisation the observed count of

outcomes is compared to an expected count obtained by
calculating the number that would have occurred if some
reference rate were applied to the observed group
denominators.

Suppose that the expected SSI rate for risk index 0

operations is 1.5%, risk index 1 operations 3%, and risk
index 2 operations 5%. Applying these rates to the Group 1
data in Table 3, the expected number of SSIs would have been
E=32×0.015+43×0.03+21×0.05=2.82. The observed number was
O=7. The methods for a single rate may now be applied giving
a mid-p 95% confidence interval for the observed count of
3.1 to 13.7. The p-value (mid-p) is 0.034.

Direct standardisation.

This requires knowledge of suitable standardising

weights for each stratum. For example, it might be known
that, for a large group of hospitals, 33.5% of the patients
requiring the operation are in risk stratum 0, 48.2% are in
risk stratum 1, and 18.3% are in risk stratum 2 (the weights
should sum to one).

Some hospitals may treat more patients in one of the

higher risk groups than others and this will bias
comparisons. However, this problem can be overcome if the
data for each hospital are weighted using the same set of

41
standardising weights.

The directly standardised rate is r=ΣWi(Xi/Ni), where

the Wi are the weights for the strata and Xi and Ni are the
numerator and denominator for stratum i. For the year 1
data, i=1, 2 and 3 and X1 is 1, X2 is 4 and X3 is 2. Also N1
is 32, N2 is 43 and N3 is 21. The value of r with these data
is r=0.335×1/32+0.482×4/43+0.183×2/21=0.073.

Confidence intervals for directly standardised rates.

A suitable estimate can be obtained using the

re-scaling method shown in Section 1 for standardised
proportions modified by employing the gamma distribution
functions as described for a single count data rate to
obtain the limits for the individual strata. Using this
method, the mid-p confidence interval is 0.042 to 0.155. For
the year 2 data in Table 3, the rate is 0.018 with the mid-p
95% confidence interval from 0.010 to 0.043.

Using R -

> load(“[Link]”
> stratifiedrates()
Enter number of strata 3
ENTER to continue
Enter numerator in stratum 1 1
Enter denominator in stratum 1 32
Enter weight for stratum 1 .335
Press ENTER for the next stratum
Enter numerator in stratum 2 4
Enter denominator in stratum 2 43
Enter weight for stratum 2 .482
Press ENTER for the next stratum
Enter numerator in stratum 3 2
Enter denominator in stratum 3 21
Enter weight for stratum 3 .183
ENTER to continue
Do you wish to review the data (y/n)? y
Weighted average = 0.07273.
Lower 95% limit = 0.04164, upper limit = 0.15467.

The weighted average of the expected SSI rates is

.335×.015+.482×.03+.183×.05=.028635. This is below the lower
confidence limit of .042.

Note that the data may be reviewed in R’s simple

spreadsheet (if necessary, answer y or Y to Do you wish to
review the data so that any errors may be corrected). The
spreadsheet is closed by clicking on the X at the top right
corner of the screen.

42
Stratified count data, differences between rates.

Indirect standardisation.

Horan and Culver (1996) describe a large sample method.

However, we employ the binomial method used for the analysis
of two independent sample rate data without stratification.
O1 is binomial with denominator O1+O2 and expected value
E1/(E1+E2). For the above data, O1=7, E1=2.82, O2=5 and
E2=8.165, E1/(E1+E2)=0.2567 and the p-value=0.033. The mid-p
value is 0.02 with the LR equal to 1/15. This indicates that
a zero difference receives 1/15 times the support of the
observed difference.

Using R -

> load(“[Link]”)
> proportion()
Enter numerator 7
Enter denominator 12
Is a reference proportion available? y
Enter reference proportion .2567
Proportion = 0.583, lower 95% limit = 0.277, upper limit = 0.848.
Mid-P 95% limits are 0.313 and 0.819.
P = 0.033, LR = 1/10, mid-p value = 0.02, LR = 1/15.

Direct standardisation.

As with proportion data, the Mantel-Haenszel (MH)

difference (Greenland and Robins 1985) can be used to
compare directly standardised rates although the method is
not strictly one of direct standardisation. The weights
chosen for each stratum are proportional to 1/(1/N1i+1/N2i)
where N1i and N2i are the person-time denominators in each
group in stratum i and the weights are scaled so that their
sum is one.

The formulas for calculating this difference, its

confidence interval and statistical significance are
D=ΣWi(X1i/N1i-X2i/N2i), DU,L=D±z√[ΣWi2(X1i/N1i2+X2i/N2i2)] and
χ2=(∑WiX1i/N1i-∑WiX2i/N2i)2/[∑Wi2(X1i+X2i)/(N1i+N2i)].

For the data in Table 3, the Mantel-Haenszel rate

difference is 0.055 with a 95% confidence interval of minus
0.002 to 0.111. The significance test result is χ2=6.77
(z=√6.77=2.6), p-value=0.009.

The rate ratio may also be calculated by the

Mantel-Haenszel method and the formula is

43
RR=[∑X1iN2i/(N1I+N2i)]/[∑X2iN1i/(N1i+N2i)]. Its confidence
limits are RRU,L=exp[loge(RR)±Z×√A], where A=B/(C×D),
B=∑[(X1i+X2i)N1iN2i/(N1i+N2i)2], C=∑X1iN2i/(N1i+N2i) and
D=∑X2iN1i/(N1i+N2i). For these data the rate ratio is 4 and
its 95% confidence interval is 1.3 to 12.5.

Although the significance test shows a small p-value,

the rate difference does not differ from zero. This
anomalous result is probably due to the very small samples
in the strata as the Mantel-Haenszel rate difference is not
a small sample estimate. If the directly standardised rates
for each year using the same weights are employed, a 95%
confidence interval for their difference can be calculated
using Newcombe’s formulas in nid() described below. This
method may be preferable when the data are sparse but it
requires further evaluation.

Using R.

> load(“[Link]”
> mhid()
Enter number of strata 3
ENTER to continue
Enter first numerator in stratum 1 1
Enter first denominator in stratum 1 32
Enter second numerator in stratum 1 0
Enter second denominator in stratum 1 97
Press ENTER for the next stratum
Enter first numerator in stratum 2 4
Enter first denominator in stratum 2 43
Enter second numerator in stratum 2 3
Enter second denominator in stratum 2 142
Press ENTER for the next stratum
Enter first numerator in stratum 3 2
Enter first denominator in stratum 3 21
Enter second numerator in stratum 3 2
Enter second denominator in stratum 3 49
ENTER to continue
Do you wish to review the data (y/n)? n
Stratified risk ratio = 4.03.
Lower 95% limit = 1.29, upper limit = 12.54.
Stratified risk difference 0.05469.
Lower 95% limit -0.00157, upper limit 0.11095.
Chisq 6.77, LR = 1/30, P = 0.009.

Note that the data may be reviewed in R’s simple

spreadsheet (answer y or Y to Do you wish to review the
data). Any errors may be corrected. The spreadsheet is
closed by clicking on the X at the top right corner of the
screen.

Very small samples.

44
 When there are few strata and the data within each
stratum are sparse, the MH procedures may give inaccurate
results. This is not uncommon in IM and QI work. The R
function nid() employs Newcombe’s square-and-add procedure
as described in Section 1 of this chapter.

Using R –

> load("[Link]")
> nid()
Enter number of strata 3
ENTER to continue
Enter first numerator in stratum 1 1
Enter first denominator in stratum 1 32
Enter second numerator in stratum 1 0
Enter second denominator in stratum 1 97
Press ENTER for the next stratum
Enter first numerator in stratum 2 4
Enter first denominator in stratum 2 43
Enter second numerator in stratum 2 3
Enter second denominator in stratum 2 142
Press ENTER for the next stratum
Enter first numerator in stratum 3 2
Enter first denominator in stratum 3 21
Enter second numerator in stratum 3 2
Enter second denominator in stratum 3 49
ENTER to continue
Do you wish to review the data (y/n)? n

Newcombe's method using exact Poisson 95% confidence limits

MH weighted average in first group 0.073

First group confidence limits 0.036 to 0.169
MH weighted average in second group 0.018
Second group confidence limits 0.007 to 0.05
MH weighted difference 0.055
Difference confidence limits 0.006 to 0.152
MH weighted ratio 4.025
Approximate 95% ratio confidence limits 1.212 to 13.373
Approximate P=0.023, approximate LR=1/13

> d1<-.073 # Confidence limits using Newcombe’s method

> d2<-.018
> u1<-.169
> l1<-.036
> u2<-.05
> l2<-.007
> dl<-(d1-d2)-((d1-l1)^2+(u2-d2)^2)^.5
> dl
[1] 0.006
> du<-(d1-d2)+((d2-l2)^2+(u1-d1)^2)^.5
> du
[1] 0.152

It is important when using a stratified analysis that

the results in each of the strata show approximately similar
effects. If there is marked inhomogeniety, for example the

45
rate difference differing from stratum to stratum, the
summary result should be ignored and the within sample
differences only should be reported. When a summary measure
is required in the presence of inhomogeniety, for example
for a meta-analysis, the DerSimonian-Laird procedure may be
employed (Kirkwood and Sterne 2003). Meta-analysis is beyond
the scope of this work. However, the R meta library or
WINPEPI may be used.

The importance of count data variation.

The formulas described above are appropriate for

performing significance tests and calculating confidence
intervals for Poisson distributed count data. The variance
of a series of counts must be approximately equal to their
mean for these methods to be appropriate, and for this to
occur, individual observations must be independent.
Bacteraemias, device related infections and medication
errors usually occur independently.

With some IM and QI data, independence cannot be

assumed. For example, the number of readmissions, pathology
tests or hospital acquired infections may be counted rather
than the number of patients being readmitted, or having
pathology tests or infections. It will usually be correct to
count the number of patients, not events.

Nevertheless, it is sometimes necessary to count the

number of events, and when this happens, methods are needed
that take into account the clustering that occurs due to
some patients having repeated admissions, tests or
infections. The statistical methods then become more complex
(Glynn and Buring 1996) and expert advice may need to be
sought.

When a patient carrying a contagious organism such as

MRSA is admitted to a ward, potential spread from this
patient will ensure that further infections in the ward do
not occur independently. In addition, many adverse outcomes
such as patient falls, readmissions, pressure ulcers and
needlestick injuries exhibit this problem because some staff
and patients are more susceptible than others.

A further difficulty is that some hospital areas have

much higher rates than others for some adverse outcomes.
When a number of different processes produce Poisson
distributed count data with different mean rates and they
are considered together, the resulting amalgamated data may
have greater variation than would be expected for a Poisson
process. For example, patients in an ICU may be more exposed

46
to and more susceptible to infection by a multiple
antibiotic resistant organism than patients in general
wards.

When this occurs, one of two different processes may

arise. First, when the data mix randomly, the result is a
linear function and it is also approximately Poisson
distributed. Second, if mixing is non-random, the result is
a weighted sum of Poissons and this is not Poisson
distributed.

As a result of the above breakdowns in independence,

variation will often be increased and employing Poisson
distribution methods will then result in low p-values and
narrow confidence limits that are incorrect. These problems
are a particular source of difficulty when control charts
are being used, and this will be dealt with further in
Chapter 4.

Frequently the negative binomial distribution will

describe these data better than the Poisson distribution.
The negative binomial distribution is also a compound
Poisson distribution with the mean count having a gamma
distribution. This property makes it valuable for modeling
count data that have large variability.

Unfortunately, there is as yet little guidance for

dealing with these data in the hospital epidemiology
literature although there is a very recent paper by
Spiegelhalter (2005). We illustrate the problem with several
trivial contrived examples. It is hoped that this discussion
will stimulate interest in this important problem among
hospital epidemiologists, and that statisticians will begin
to describe accurate simple approximate methods for
analysing these and similar data.

Suppose that a department has F=4 needlestick injuries

in one month and its staff wish to calculate a confidence
interval for this count. Using the Poisson methods described
earlier in this chapter, the mid-p 95% confidence interval
is calculated to be 1.35 to 9.51. However, suppose it is
known that when the mean monthly count for needlestick
injuries is M=6, its variance is V=12 and not approximately
6 as would be expected for Poisson distributed data.

The negative binomial distribution has 2 parameters

that we label S and P. S represents a fixed number of
successes in a series of Bernoulli trials and P the
probability of success in each trial. We label the observed
value F for failures. Unlike the binomial distribution where

47
S is a variable and S+F=N is fixed, the negative binomial
variable is F and S is fixed.

When data follow a negative binomial distribution, we

can calculate S=M2/(V-M) and P=M/V. For the above data,
S=62/(12-6)=6, and F=4. We can now use the beta distribution
to obtain mid-p 95% confidence limits for F. Let
A=qbeta(0.025,S,F) and B=qbeta(0.025,S,F+1). Then the upper mid-p
95% confidence interval for F is U=[(S-A×S)/A+(S-B×S)/B]/2.
The lower limit L is found by substituting 0.975 for 0.025
in the above formulas. The mid-p 95% confidence limits are
1.167 to 15.457. These are much wider than the corresponding
Poisson limits.

It is possible to test whether an observed monthly

count differs from a mean value. For example, suppose that
the count for the current month is 12. Does this differ
significantly from 6? Using the conventional count data
method that employs the Poisson distribution, the two-tailed
p-value is 0.04 so that a conventional test would suggest a
statistically significant difference.

When the negative binomial distribution is employed the

p-value is 0.12. When M=6 and V=12, P=0.5.
For F greater than M, pval={[1-pbeta(P,S,F)]+[1-pbeta(P,S,F+1)]}/2
and, for F less than M, p-val=[pbeta(P,S,F)+pbeta(P,S,F+1)]/2.

The R commands and output for F=4 and F=12 are as

follows –

> load(“[Link]”)
> overdispersed1()
Single sample overdispersed count data.
Enter observed count 4
Enter reference mean 6
Enter reference variance 12
Lower 95% confidence limit = 0.952, upper confidence limit = 16.868.
Lower mid-p value = 1.167, upper mid-p value = 15.457.
LR = 1/1, z = 0.313, p-value = 0.754.
Mid-P LR = 1/1, mid-p z = 0.481, p-value = 0.631.

> overdispersed1()
Single sample overdispersed count data.
Enter observed count 12
Enter reference mean 6
Enter reference variance 12
Lower 95% confidence limit = 4.722, upper confidence limit = 38.968.
Lower mid-p value = 4.971, upper mid-p value = 37.596.
LR = 1/3, z = 1.463, p-value = 0.143.
Mid-P LR = 1/3, mid-p z = 1.555, p-value = 0.12.

Suppose now that it is of interest to determine whether

the counts of 4 and 12 differ. In addition, it would be

48
desirable to be able to calculate a confidence interval for
this difference. When we employ a Bayesian method using
diffuse beta(1,1) priors and beta distributed random numbers
(Antelman 1997), the 95% credible interval is –3.94 to
26.47, p-val=0.186. We can conclude that the difference is
not significant. The corresponding Poisson derived interval
is 0.16 to 15.84, p-val=0.046. Thus it can be seen that the
Poisson methods underestimate p-values and produce
confidence intervals that are too narrow when the variances
of count data exceed their mean values.

The R commands and output are shown below. When the

function is called, the R spreadsheet appears. The two
observed counts are entered on the first line. The reference
means are then entered on the next line and the variances on
the following one. The reference means and variances may
differ. Finally, the exposure is entered on the final line.
For these data, the exposures were the same and no
person-time or bed-days denominators were available so the
default values of 1 were left unchanged. When the data have
been entered, click on the X at the top right corner of the
screen to return to R.

load(“[Link]”)
overdispersed2()
Difference = -8.
95% confidence limits are 3.93855 and -26.46584.
LR = 1/2, z = 1.324, p-value = 0.186.

There is currently much misguided interest in comparing

outcomes and some managers may be tempted to compare the
work of two departments such as ICUs using surveillance
data, for example by analysing weekly or monthly counts of
MROs. Usually such comparisons would entail employing the
common Poisson count data methods that may be incorrect, and
may suggest that there is a difference when one does not
exist.

Such departments may have markedly differing patient

populations and risk adjustment may be difficult or
impossible. Moreover the distributions for the two
departments may differ so that the parameter S may not be
the same for each of them. It may be difficult to obtain
sufficient monthly or weekly data during periods marked only
by endemic activity to obtain stable estimates of mean
values and variances. In addition, the departments may be of
differing sizes so that obtaining comparable denominators
may be difficult. Methods have been described for dealing
with negative binomial data having denominators of varying
sizes (Bissell 1972). However, as described in Chapter 2, it

49
is often difficult to determine the correct denominators.
The denominator employed is often occupied bed-days. There
are 2 difficulties with this. First, more severely ill
patients are usually more susceptible and their bed-days may
fluctuate less than overall bed-days. Secondly, LOS which
influences occupied bed-days is also often a risk factor.

If comparisons of this nature must be made, it is much

more likely that sensible results will be obtained if system
and process surveillance is performed and the resulting data
used for making comparisons. For example, is hand washing
performed with equal frequency by the staff of the two
departments? Are there in place sensible guidelines for the
use of antibiotics and are these guidelines implemented? Is
there an effective program for isolation of patients
carrying potentially dangerous MROs? Are current guidelines
for the management and prevention of nosocomial infections
up to date and are they being followed correctly?

To use negative binomial methods, it is necessary to

have data that enable mean outcome rates and their variances
to be calculated. For monthly counts it would seem desirable
to have at least 20 available from a period known to be free
of any epidemic activity.

We must gain much more experience at understanding how

variation occurs with count data adverse events especially
colonisations and infections due to MROs. It is possible to
see variances that are similar to mean values with new
isolates of some MROs when isolation procedures, antibiotic
use, and hand washing and other hygiene measures are
excellent and the MRO burden is low. However, if there is
less than excellent IM, for example there may be
insufficient numbers of isolation beds, variability can
increase even in the absence of an epidemic. This is a
difficult area and much work is needed to improve
understanding. In some cases, it may be useful to perform a
sensitivity analysis.

Count data nosocomial infections are often better

analysed in control charts that display monthly occurrence
rates. Good departments examine their systems and processes
regularly and constantly strive to improve them, using for
example the Deming cycle and QI tools described in Chapter
1. When this is done, the sequential analysis of outcome
data using control charts that when appropriate employ
negative binomial control limits, provides reliable
information for making decisions. These methods are
described in Chapter 4.

50
 Management frequently tries to employ statistical
methods to prevent substandard performance, for example
using league tables and star ratings. However, to prevent
substandard performance it is necessary for management to
have a firm understanding of systems and the courage,
resources and authority to act when there is an environment
that is known to be associated with unsatisfactory outcomes.

Consider the situation where one death is expected. It

takes at least four to occur in a similar time period for
statistical significance to be attained. Reliance on
statistical analysis will result in unacceptable delay if it
is possible, employing systems analysis, to determine that
an environment exists in which excess deaths may be expected
to occur. However, when an institution first analyses its
systems and institutes measures to improve, then monitoring
and analysis of relevant data using statistical methods such
as control charts can be an invaluable adjunct to systems
analysis and improvement.

Funnel plots for count data from several institutions when

there is overdispersion.

Funnel plots for summarising count data have been

described for data that are not overdispersed earlier in
this section. There we referred to the file [Link]
that contains bacteraemia data for a number of hospitals,
some of which are level 1 and others are level 2 or 3
institutions. The funnel plot for the combined data (Figure
4) shows marked overdispersion. The plots for Level 1
Hospitals (Figure 5) and for Level 2 & 3 Hospitals combined
(Figure 6) do not have this problem. However, to illustrate
the method for dealing with overdispersion using the
negative binomial distribution, we employ the R function
[Link] to analyse these data and the result is
shown in Figure 7. We do not recommend mixing distributions
in this manner and use these data only to illustrate the
method. When displaying data from a number of hospitals, the
data should be independent and overdispersion is likely to
be due to heterogeneity of the type described above. We will
want to explore this heterogeneity rather than compensate
for it statistically.

Using R –
> groupfunneln()

Do you have risk-adjusted data (Y/N) n

Loading data.
Data from clipboard (C) or file (F) f
Have you selected the correct directory (Y/N) y

51
Enter the name of the file [Link]
Do data column(s) have heading(s) (Y/N) y
Are data columns separated by spaces (S) or commas (C) c
The mean count is 24.167.
Do you want to change the mean value (Y/N) n
The mean denominator is 40859.
Do you want to change the mean denominator value (Y/N) n
The weighted variance (V) is 199.663.
Do you want to change V (Y/N) n
Denominators.
Enter name for denominator (eg Bed-days) Bed-days
Denominators.
Per thousand (1000), hundred (100), ten (10), unit (1) Bed-days 1000
Chart heading.
Enter heading for chart eg MRSA by hospital Bacteraemias all hospitals

# To calculate weighted mean and weighted variance. The vector d1 is for

# the counts of the outcome and d2 is for the corresponding
# denominators. In the function [Link].
> load("[Link]")
> [Link]()
Loading data.
Data from clipboard (C) or file (F) f
Have you selected the correct directory (Y/N) y
Enter the name of the file g:/examples/[Link]
Do data column(s) have heading(s) (Y/N) y
Are data columns separated by spaces (S) or commas (C) c
> datain
HospitalID Bacteraemias [Link]
1 15 5 19716
2 1 23 40587
3 2 22 76995
4 17 2 10815
5 3 3 34303
6 18 8 33722
7 19 7 17016
8 20 1 8271
9 5 28 22822
10 4 1 13409
11 6 8 40829
12 7 101 105308
13 23 14 29390
14 9 7 21893
15 10 116 118939
16 8 37 59239
17 13 40 50489
18 12 12 31713
> d1<-datain[,2];d2<-datain[,3]
> m<-mean(d1)
> dbar<-mean(d2)
> dw<-d2/dbar
> v<-sum(dw*(d1/dw-m)^2/(length(dw)-1))
> a1<-paste("The mean count is",[Link](round(m,3)),".\n")
> a2<-paste("The mean denominator is",[Link](round(dbar,0)),".\n")
> a3<-paste("The weighted variance (V)
is",[Link](round(v,3)),".\n")
> cat(a1)
The mean count is 24.167.
> cat(a2)
The mean denominator is 40859.

52
> cat(a3)
The weighted variance (V) is 199.663.

Using graphics.

As described in the section on proportion data,

graphical display of rates and their differences is an
important method of data presentation. The methods for
proportions, illustrated in the previous section of this
chapter, may also be used for count data rates. For example,
Figure 8 shows the readmission rate data per 1000 occupied
bed-days for each of the years (12 in 32190 and 28 in 37440
bed-days), and for their difference, together with 95%
confidence limits.

Figure 8 was produced with the following R commands –

> a<-c(.37,.75,.38)
> b<-c(.65,1.08,.73)
> d<-c(.19,.5,.02)
> x<-c(1,2,3)
> l<-c("Year 1","","Year 2","","Diff")
> plot(a,ylim=c(0,1.2),type="p",axes=F,xlab="",ylab="Per 1000 occupied
bed-days")
> box()
> axis(side=1,tick=F,labels=l)
> axis(side=2)
> arrows(x,b,x,d,angle=90,code=3,col="red")
> title(main="Unplanned readmission rates, rate difference and\n95%
confidence intervals.")

Chi-squared and trend tests for larger tables.

Sometimes data will be in a large contingency table.

For example, an audit of hand-washing may have been
conducted in 4 wards, classifying the staff as –

1. Usually washing hands between patients,

2. Washing hands between patients about half the time, and
3. Infrequently washing hands between patients.

Suppose that, for Ward 1 the numbers of staff for

Usually, Average, and Infrequently were 12, 8, and 3
respectively. For the remaining 3 wards the respective
numbers were 5, 7, and 9 for Ward 2; 4, 17, and 11 for Ward
3; and for Ward 4 they were 6, 15, and 10 (hypothetical
data).

Using R –

> load("[Link]")
> largetable()

53
Enter the number of columns in the table 4
Enter the number of rows in the table 3 # The R spreadsheet appears
Observed values. # Enter data and close
A B C D
1 12 5 4 6
2 8 7 17 15
3 3 9 11 10

Expected values.
A B C D
1 5.80 5.30 8.07 7.82
2 10.10 9.22 14.06 13.62
3 7.09 6.48 9.87 9.56

Residuals
A B C D
1 3.36 -0.17 -1.98 -0.89
2 -1.00 -1.09 1.25 0.59
3 -2.08 1.33 0.52 0.20

Unordered data 2. Ordered rows 3. Ordered rows & columns Enter 1, 2 or 3

2 # The data in the rows are ordered
Pearson's Chi-squared test.
Chisq = 14.339, DF = 6, P-value = 0.026.
Kruskal-Wallis rank sum test .
Chisq = 10.44, DF = 3, P-value = 0.015.
Group ranks 36.91 59.33 60.12 56.74.

The Fisher Exact test is an alternative to the

chi-squared test when samples are small.

Using R –

a<-c(12,8,30);b<-c(5,7,9);d<-c(4,17,11);e<-c(6,15,10)
> F<-[Link](a,b,d,e)
> F
a b d e
1 12 5 4 6
2 8 7 17 15
3 30 9 11 10
> [Link](F)

Fisher's Exact Test for Count Data

data: F
p-value = 0.01179
alternative hypothesis: [Link]

When the data are in ordered categories, the

Kruskal-Wallis test often has more power to detect a
difference than the global Pearson’s chi-squared test.
Examination of the average ranks, expected values and the
standardised residuals suggests that Ward 1 differs from the
others. This can be examined further by amalgamating the
data for the other 3 wards and comparing the amalgamated
data with that of Ward 1.

54
 For the differences among the 3 wards (B, C and D) the
chi-squared result using largetable() is 2.5 with 4 DF and a
p-value of 0.647, so that there is no evidence of a
difference among them. Next the data for these 3 wards can
be amalgamated and compared with the Ward 1 data also using
largetable() - the resulting chi-squared value is 11.99 with
2 DF and a p-value of 0.0025. This suggests that the
difference would have been between Ward 1 and the other
wards.

However, suppose that, in conducting the audit, you had

classified the wards by the quality of their handwashing
facilities and found that Ward 1 had the best facilities,
followed by Wards 2, 3, and 4 respectively. A trend test
might then be warranted to examine the relationship of hand
washing to hand washing facilities. Running the function
again and selecting option 3 for trends in both variables
gives the following result –

1. Unordered data 2. Ordered rows 3. Ordered rows & columns Enter 1, 2

or 3 3
Pearson's Chi-squared test.
Chisq = 14.339, DF = 6, p-value = 0.026. Trend test chisq = 4.24, p-
value = 0.04. Kendall's Tau = 0.17.
Departure from trend chisq = 10.101, DF = 5, p-value = 0.072.

The evidence for a trend is not strong. Had these data

been genuine, the analysis would have suggested that Ward 1
would have been superior rather than the differences among
the wards being due to a trend associated with the quality
of their handwashing facilities. An alternative explanation
would be that there was better leadership and teamwork
concerning infection prevention in Ward 1.

In some cases, one may want to look further, for

example, for differences among the wards while controlling
for differences in the handwashing facilities. Although it
is possible to conduct a stratified analysis, it will
usually be best to use an advanced statistical method such
as Poisson regression that is described in Chapter 5.

In summary, this section has dealt with statistical

methods for the analysis of count data. The following
aspects have been addressed –

1. Statistical methods for counts and rates.

2. Dealing with excessive variation
3. Graphical methods for counts and rates.
4. Large contingency tables.

55
 The next section of the chapter deals with data in
ordered categories.

Section 3
Categorical data
Introduction

In this section we deal with categorical data,

especially data in ordered categories. For an estimator for
these data, we describe the RIDIT, also called the C
statistic and AUC. Somers’ delta and the RIDIT are
equivalent and we describe their relationship. In addition,
we describe a method for calculating the RIDIT confidence
interval. The significance test for these estimators is the
rank sum test and we describe a method for its calculation.
Next, we extend these methods to stratified data using
standardisation. Finally, we show that graphical methods are
useful for summarising these data.

The importance of these methods stems from their

usefulness in analysing adverse events that have “near
misses”, patient and staff satisfaction surveys and LOS
data. The latter can be thought of as ordered categorical
data because of the many tied values for shorter lengths of
stay within most DRGs. Severe adverse patient outcomes are
often preceded or accompanied by more minor adverse events
that can be thought of as near misses. Incorporating the
additional information of these near misses will often
increase the power of an investigation, permitting earlier
detection of problems and reducing the possibility of false
negative states due to undetected quality difficulties.

Categorical data methods are used infrequently in QI

studies at present. However, there has been considerable
interest in the C statistic that is used to estimate AUC,
the area under Receiver Operating Characteristic (ROC)
curves. This estimator is equivalent to the mean RIDIT, the
estimator for ordered categorical data. We shall return to
these estimators in Chapter 5 that deals with logistic
regression.

In most cases where categorical QI and IM data are

analysed, they will be in ordered categories and this
section will concentrate on methods for ordered categorical
data.

56
 As described above, it frequently happens that an
increase in the rate of some complication is preceded or
accompanied by an increase in the number of near misses. For
example, surgical site infections are frequently classified
as superficial or deep/organ space, the latter being the
more serious, and an increase in deep infections may be
preceded by or accompanied by an increase in superficial
infections. Similarly, an increase in postoperative deaths
following a major operation such as colectomy,
pancreatectomy, oesophagectomy, or abdominal aneurysm repair
is likely to be preceded and accompanied by an increase in
major postoperative complications such as sepsis,
haemorrhage, anastomotic breakdown, or organ injury or
failure, or increases in returns to the operating theater or
prolonged ICU stay.

These ideas lead naturally to the use of ordered

categorical methods. For example, uncomplicated wounds,
superficial SSIs, and deep/organ space SSIs are ordered in
severity but no definite numerical classification is
possible. Making use of this ordered category information
will usually increase the power of an investigation and
improve our ability to detect statistically significant
changes.

The term RIDIT, coined by Bross (1958), stands for

“relative to an identified distribution”. For ordered
categorical data in two independent samples, the mean RIDIT
is equivalent to the C statistic, although they are usually
calculated differently. Both the mean RIDIT and the
equivalent C statistic are estimators for the Wilcoxon rank
sum test (Selvin 1977). They are also equivalent to another
estimator in this group called Somers’ delta (Reynolds 1977)
except that the coefficients for the former two estimators
lie between zero and one while those for the latter lie
between minus one and plus one. Thus, when the mean RIDIT or
C is zero, Somers’ delta is minus one, when they are 0.5
Somers’ delta is zero, and when they are one, Somers’ delta
is also one. Somers’ delta may be thought of as a
generalisation of the difference between two independent
proportions to data in ordered categories.

RIDIT≡C≡AUC=Delta/2+0.5.

The mean RIDIT is the estimated probability that a

randomly chosen individual from a comparison group will be
more affected than a person selected at random from a
reference group (Selvin 1977). The C statistic equivalently
represents the probability that a randomly chosen diseased
subject is correctly ranked with greater suspicion than a

57
randomly chosen non-diseased subject (Hanley and McNeil
1982). These estimators also measure AUC the area under the
ROC curve that is frequently used to assess discrimination
in logistic regression. In the latter, the cases are
diseased and the non-cases or controls non-diseased and the
data in the ordered categories are the probabilities from
the logistic regression equation. Used in this way, C and
the ROC curve are important for studying sensitivity and
specificity, and discrimination (Altman, Machin, Bryant and
Gardner 2000); this latter use will be discussed further in
Chapter 5 which deals with logistic regression.

Another way of thinking of a RIDIT is as the score for

each ordered category that is equivalent to a percentile
rank (Kantor, Winkelstein, and Ibrahim 1968). This aspect is
of particular importance when describing LOS and similar
data. The distributions of these data are markedly skewed so
that using means and standard deviations as summary measures
may be misleading (Davies 1998). In addition, LOS data
within most DRGs typically have many values tied at lower
lengths of stay so that these data can be regarded as
categorical. Data that are markedly skewed or categorical
are often better summarised by medians and percentiles
(Davies 1998), or by RIDITS. Usually, with these data, tests
of significance are best performed by the rank sum test
provided the data in each study group have similar shaped
distributions.

There is a substantial literature on estimators for the

rank sum test. The Hodges-Lehmann estimators (Altman,
Machin, Bryant and Gardner 2000) are useful for data that
are numerical. However, they are limited by the fact that
the one sample version is sensitive to asymmetry in the
data. Data such as LOS and cost of stay are usually strongly
positively skewed with occasional marked outlier values
causing severe asymmetry. In addition, the presence of
numerous tied values often makes Hodges-Lehmann estimators
less attractive. For these reasons, we do not describe the
Hodges-Lehmann estimators further.

Comparison of a single sample with a reference value.

Example (hypothetical data).

The RIDIT method is best illustrated by example.

Suppose a hospital surveillance program has detected 2
deaths and 5 major postoperative complications among 50 ASA
class 1 and class 2 patients undergoing a major elective
surgical procedure such as abdominal aortic aneurysm repair.
An examination of 200 consecutive similar operations

58
performed when complication rates were considered acceptable
might have revealed 2 deaths and 4 major postoperative
complications. The number of deaths from the large series
would suggest that a mortality rate of 1% might be
acceptable for this procedure among ASA class 1 and class 2
patients undergoing elective surgery; with this expected
rate, the 2 observed deaths in 50 operations would not be
statistically significant (p–value=0.14). However, this
analysis neglects to use the information in the intermediate
group (postoperative complications) and making use of these
extra data increases the power of the test and thus enables
earlier detection of adverse changes.

Proceed as shown in Table 5 to calculate the reference

data RIDIT scores (Fleiss 1981). The numbers in column 1 of
Table 5 are divided by 2 to give those in column 2. The
accumulated entries in column 1 are in column 3 one row
down. Column 4 is calculated by summing columns 2 and 3.
Finally the RIDIT score for each category is obtained by
dividing the result in column 4 by N, the sample size which
in this case is 200.

Now the corresponding RIDIT values for the observed

data are calculated. There were 2 deaths and 5 postoperative
complications so that there were 43 uncomplicated outcomes.
The mean RIDIT for the observed data is
(43×0.485+5×0.98+2×0.995)/50=0.555. The mean RIDIT for the
reference group is always 0.5. The best way to test the
significance of the mean RIDIT for the observed data is to
use the rank sum test for which the p-value=0.002. Its
formula is Z=(|Nc-Nd|-cc)/√V, where
V=N1N2(N3-N-Σ(Ui3-Ui))/[3N(N-1)], N1 and N2 are the numbers of
observations in each group and N=N1+N2. Nc and Nd are the
numbers of concordant and discordant pairs, Ui is the number
of values in both groups tied in category i,
cc=(2N-U1-Uk)/[2(k-1)] is a continuity correction and k is
the number of categories.

Using R –

> load(“[Link]”)
> categoryranksum()
Enter the number of data categories 3

Group1 Group2
1 194 43
2 4 5
3 2 2

Wilcoxon Rank Sum Test p-value = 0.002.

RIDIT Score = 0.555, 95% confidence limits are 0.505 to 0.605.

59
Rank sum variance = 123810.1, RIDIT Variance = 0.0006498741.

When the function categoryranksum() is called, the R

spreadsheet appears for the data to be entered. When this is
complete, the spreadsheet is closed by clicking on the X at
the top right corner of the screen. The rank sum and RIDIT
variances are provided as then may be required for analysing
stratified data (see below).

A difficulty with these data is the small numbers in

some of the categories. However, a chi-squared test using
simulation gives some reassurance, although it fails to take
the ordered nature of the data into account. Also two
proportions can be obtained by amalgamating the near miss
and major adverse event categories to give 6/200 and 7/50
and these differ.

Using R -
> a<-c(194,43) # Test with simulated p-values
> b<-c(4,5)
> d<-c(2,2)
> x<-[Link](a,b,d)
> [Link](x,[Link]=T,B=10000)

Pearson's Chi-squared test with simulated p-value (based on

10000 replicates)

data: x
X-squared = 9.8717, df = NA, p-value = 0.0063.

The Fisher Exact test may be used as an alternative to

the chi-squared test.

> a<-c(194,43) #Fisher Exact test

> b<-c(4,5)
> d<-c(2,2)
> x<-[Link](a,b,d)
> [Link](x)

Fisher's Exact Test for Count Data

data: x
p-value = 0.00775
alternative hypothesis: [Link])

>load(“[Link]”) # near miss and death combined 7/50 v 6/200

> twoproportions()
Enter first numerator 6
Enter first denominator 200
Enter second numerator 7
Enter second denominator 50
Difference between proportions 0.11.
Lower 95% limit 0.029, upper limit 0.227.
Z = 2.72, P = 0.006, LR = 1/41.

60
Ratio 4.67.
Ratio 95% confidence limits are 1.69 and 12.63.

A pair of values is concordant if the one in the group

with the smaller values, that is the one with the larger
numbers in the lower categories, is less than the one in the
other group. They are discordant if the reverse is true, and
they are tied if equal. We illustrate these concepts below.

Altman, Machin, Bryant and Gardner (2000) describe a

method for calculating AUC or the mean RIDIT and it is
employed in the R function categoryranksum(). For each value
(j=1..J) in the first group f calculate Xj=Σ(fj>s++.5fj=s+)
and for each value (j=1..J) in the second group s calculate
Yj=Σ(sj>f++.5sj=f+), where s+ and f+ stand for all the values
in the first and second groups respectively. Then
AUC=ΣXj/N1N2 and its non-null variance is V=S2x/N1+S2y/N2,
where S2x=Σ((Xj/N2-AUC)2)/(N1-1) and S2y=Σ((Yj/N1-AUC)2)/(N2-1).

An alternative approach that may have more appeal to

some workers because of its relationship to the rank sum
test is now shown. Here, AUC=(M2-M1)/(N1+N2), where M1 and M2
are the mean ranks of the two groups and N1 and N2 the
numbers of observations.

Using R –

> x<-c(rep(1,194),rep(2,4),rep(3,2))
> y<-c(rep(1,43),rep(2,5),rep(3,2))
> x1<-rep(1,length(x))
> y1<-rep(2,length(y))
> xy<-c(x,y)
> xy1<-c(x1,y1)
> z<-[Link](xy1,xy)
> zz<-rank(z[,2])
> s<-[Link](z[,1],zz)
> m1<-mean(s[,2][s[,1]==1])
> m2<-mean(s[,2][s[,1]==2])
> auc<-(m2-m1)/length(s[,1])+.5
> auc
[1] 0.5549
> v1<-var(s[,2][s[,1]==1])/length(s[,1][s[,1]==1])
> v2<-var(s[,2][s[,1]==2])/length(s[,1][s[,1]==2])
> v<-(v1+v2)/length(s[,1])^2
> v # variance of auc
[1] 0.0006499546 # 0.0006498741 by previous method

Calculation of the variance of the mean RIDIT is

complicated by the fact that two formulas are required – one
for a significance test called a null variance and one for
calculation of confidence intervals called a non-null
variance. Since the rank sum test is used to test whether

61
the mean RIDIT differs from zero, the null variance of the
mean RIDIT is of less importance for hypothesis testing.
However, it is required when these methods are employed in
control charts. A simple approximate null and non-null
variance formula for a mean RIDIT (Fleiss 1981) is
V=1/12N1+1/12N2, where N1 and N2 are the sizes of the
observed samples (Fleiss also describes a more complicated,
and slightly more accurate formula). Sometimes the reference
sample is very large and then the formula will be V≈1/12N
approximately where N is the size of the observed sample. We
shall employ the latter formula in control charts for
ordered categorical data that are described in Chapter 4.
The non-null variance is required to calculate confidence
limits for the mean RIDIT and when the data are stratified
(see below). For these data the non-null variance is
0.00065. This gives a 95% confidence interval of 0.505 to
0.605.

Comparison of two independent samples.

The RIDIT methods described above can be used for

studying two independent samples of categorical data.
Another method similar to that described for calculating the
C statistic is also of interest. The methods give identical
results and derive from the fact that the rank sum test can
be performed either by ranking the data or by considering
whether the individual data values in the two samples are
concordant or discordant. As we have described, a pair of
values is concordant if the one in the group with the
smaller values, that is the one with the larger numbers in
the lower categories, is less than the one in the other
group. They are discordant if the reverse is true, and they
are tied if equal. The alternative formula for the mean
RIDIT is (Nc-Nd)/[2×(Nc+Nd+Nt)]+0.5, where Nc is the number of
concordant pairs, Nd is the number of discordant pairs and
Nt is the number of tied pairs. When comparing two groups,
it is usually convenient to employ the larger group as the
reference sample. If the two groups contain roughly equal
numbers, use the one with the smaller values, that is the
one that has the larger numbers in the lower categories, as
the reference group.

Example.

Suppose that following the 50 operations described

above, when there was evidence of an impending quality
problem, a thorough reassessment of the relevant system and
processes of patient care was undertaken. Suppose also that
surveillance of the next 50 operations disclosed 1

62
postoperative complication and 0 deaths, there being 49
normal outcomes (Table 6).

We now describe how the numbers of concordant and

discordant pairs are calculated. We use the second group
with the bigger numbers of data values in the lower
categories as the reference group. There are
Nc=49×(5+2)+1×2=345 concordant pairs. There are 43 values in
the first group which are smaller (in the sense that normal
is lower than intermediate and intermediate that represents
a major postoperative complication without death is lower
than death) than one in the second group so that Nd=43 and
Nc-Nd=345-43=302. In addition, we require Nt, the number of
tied pairs, Nt=43×49+5×1+2×0=2112. Therefore
Nc+Nd+Nt=345+49+2112=2506 so that the mean
RIDIT=302/(2×2506)+0.5=0.56. The rank sum test p-value is
0.027 showing that the change, if it really were to have
occurred, would be unlikely to have arisen by chance. Had
consideration been given solely to the deaths, the
difference between 2 in 50 operations and 0 in 50 is not
statistically significant (P = 0.25).

Using R -

> load (“[Link]”)

> categoryranksum()
Enter the number of data categories 3

Group1 Group2
1 49 43
2 1 5
3 0 2

Wilcoxon Rank Sum Test p-value = 0.027 .

RIDIT Score = 0.56 , 95% confidence limits are 0.508 to 0.613 .
Rank sum variance = 18610.1 , RIDIT Variance = 0.0007135445 .

With the above data the numbers in several of the

categories are very small so the statistical significance
result should be viewed with caution although practically
important improvement certainly would appear to have
occurred, had these data been real. Unfortunately,
[Link]() in R cannot provide exact results when there
are ties in the data. The chi-squared test with simulation
and the Fisher Exact test just fail to reach conventional
statistical significance; however, they are likely to be
inefficient when the data categories are ordered. The
difference between the proportions 7/50 and 1/50 (the near
miss and death cateories have been amalgamated) suggests a
difference that is sufficiently large for random variation
to be an unlikely explanation for its occurrence.

63
Using R -

> a<-c(43,49) # chi-squared test

> b<-c(5,1)
> d<-c(2,0)
> x<-[Link](a,b,d)
> x
a b d
1 43 5 2
2 49 1 0
> [Link](x,[Link]=T,B=10000)

Pearson's Chi-squared test with simulated p-value (based on

10000 replicates)

data: x
X-squared = 5.058, df = NA, p-value = 0.06729.

> a<-c(43,49) # Fisher Exact test

> b<-c(5,1)
> d<-c(2,0)
> x<-[Link](a,b,d)
> [Link](x)

Fisher's Exact Test for Count Data

data: x
p-value = 0.06692
alternative hypothesis: [Link]

> twoproportions()
Enter first numerator 1
Enter first denominator 50
Enter second numerator 7
Enter second denominator 50
Difference between proportions 0.12.
Lower 95% limit 0.018, upper limit 0.237.
Z = 2.21, P = 0.027, LR = 1/12.
Fisher Exact p-value = 0.059. # mid-p value 0.033
Ratio 7.
Ratio 95% confidence limits are 1.19 and 42.9.

The mean RIDIT and its variance may be calculated from

the ranks as described above.

Using R -

> x<-c(rep(1,49),rep(2,1))
> y<-c(rep(1,43),rep(2,5),rep(3,2))
> x1<-rep(1,length(x))
> y1<-rep(2,length(y))
> xy<-c(x,y)
> xy1<-c(x1,y1)
> z<-[Link](xy1,xy)
> zz<-rank(z[,2])
> s<-[Link](z[,1],zz)
> m1<-mean(s[,2][s[,1]==1])

64
> m2<-mean(s[,2][s[,1]==2])
> auc<-(m2-m1)/length(s[,1])+.5
> auc
[1] 0.5604
> v1<-var(s[,2][s[,1]==1])/length(s[,1][s[,1]==1])
> v2<-var(s[,2][s[,1]==2])/length(s[,1][s[,1]==2])
> v<-(v1+v2)/length(s[,1])^2
> v
[1] 0.0007147739

Stratified data.

The above methods may be used with stratified data and

direct standardisation. For example, patients undergoing
surgery may be stratified by ASA class. ASA classes 1 and 2
and ASA classes 3, 4, and 5 may be amalgamated giving rise
to two strata.

The rank sum test for stratified data is

Z=Σ(Nci-Ndi)/√ΣVi, where Vi is the variance for stratum i.
The continuity correction is not used with stratified data.

The mean RIDIT can be similarly generalised to

ΣWiRbari, where the weights (Wi) are scaled to sum to one and
Rbar stands for the mean RIDIT. Also its variance is V=ΣWi2Vi
so that the confidence limits for a directly standardised
mean RIDIT are Rbar±Z√V.

Note that, when calculated using reference database

data, RIDIT scores have an expected value of 0.5 and an
approximate variance of 1/12. Therefore, RIDIT scores from
stratified data can be combined by adding them and obtaining
a mean RIDIT value by dividing by the number added. The
approximate variance value is then 1/12N. The expected value
remains 0.5. This is analogous to using indirect
standardisation employing RIDIT scores and it is useful for
combining data for LOS described in the Section 4 of this
chapter and, in relation to control charts, in Chapter 4.

The hypothetical data in Table 7 might be those for

similar surgery on ASA class 3-5 patients performed at the
same time as that described above for ASA class 1-2 patients
(Table 6). The weights would be 0.6 and 0.4 for ASA class
1-2 and ASA class 3-5 respectively if there were, in a
reference series of operations, 6 of the former operations
to every 4 of the latter.

Using R -

> load(“[Link]”) # Table 7 data

> categoryranksum()

65
Enter the number of data categories 3

Group1 Group2
1 23 15
2 2 4
3 1 2

Wilcoxon Rank Sum Test p-value = 0.151.

RIDIT Score = 0.585, 95% confidence limits are 0.468 to 0.703.
Rank sum variance = 4100.303, RIDIT Variance = 0.003603228.

> x<-c(rep(1,23),rep(2,2),rep(3,1)) # alternative method using ranks

> y<-c(rep(1,15),rep(2,4),rep(3,2))
> x1<-rep(1,length(x))
> y1<-rep(2,length(y))
> xy<-c(x,y)
> xy1<-c(x1,y1)
> z<-[Link](xy1,xy)
> zz<-rank(z[,2])
> s<-[Link](z[,1],zz)
> m1<-mean(s[,2][s[,1]==1])
> m2<-mean(s[,2][s[,1]==2])
> auc<-(m2-m1)/length(s[,1])+.5
> auc
[1] 0.5851648
> v1<-var(s[,2][s[,1]==1])/length(s[,1][s[,1]==1])
> v2<-var(s[,2][s[,1]==2])/length(s[,1][s[,1]==2])
> v<-(v1+v2)/length(s[,1])^2
> v
[1] 0.003610234

The calculations for obtaining a directly standardised

mean RIDIT and the rank sum test values are summarised as
follows.

Table 6 data (ASA class 1 and 2).

Nc-Nd=302 as already described, the continuity correction

cc=26.5, where N=100 is the size of the combined sample,
U1=92 and Uk=2 are the values tied in the first and last
categories, and k=3 is the number of categories in the
table. In addition, the rank sum variance V=18610.1, where
N1=50 and N2=50 are the size of the respective samples and Ui
is the number tied at category i (92, 6 and 2 respectively),
and the p-value=0.027. In addition the mean RIDIT=0.56 and
its non-null variance V=0.000714, so that the 95% confidence
interval is 0.508 to 0.613.

Table 7 data (ASA class 3, 4, and 5).

Nc-Nd=93, the continuity correction cc=13.25, the rank

sum variance V=4100.3, p-value=0.151. In addition, the mean
RIDIT=0.585, and its non-null variance V=0.003603, so that
its 95% confidence interval is 0.468 to 0.703. The
calculation of the rank sum variance in R is as follows -

66
> n1<-21
> n2<-26
> n<-n1+n2
> u<-c(38,6,3)# the number of tied values in each row
> v<-(n1*n2*(n^3-sum(u^3)))/(3*n*(n-1))
> v
[1] 4100.303

For the combined data.

Z=(302+93)/√(196101.1+4100.3)=2.62, p-value=0.009. In
addition, the mean RIDIT Rbar=0.6*0.56+0.4*0.585=0.57 and
its variance V=0.62×0.000714+0.42×0.003603=0.000834. The 95%
confidence interval is Rbar±Z√V=0.514 to 0.625, where
Z=1.96. Thus, even though the result for the second stratum
is not significant, by using the data in both strata a more
powerful analysis is possible.

Using graphics.

Graphical display of RIDITS and their confidence

intervals is an important aspect of their presentation.
Figure 9 is an illustration of the application of these
methods to LOS data. These data have already been discussed
in the first section of this chapter where they were used to
illustrate graphical methods for proportions. LOS data are
very suitable for analysis using RIDITS because their
distributions are highly skewed and, within a DRG, there are
many tied values at most shorter lengths of stay so that
methods for ordered categorical data are therefore
appropriate. The use of RIDITs for analysing LOS data will
be dealt with further in the next section of this chapter.

Figure 9 shows the mean RIDIT and its confidence limits

for LOS for DRG 252 (Heart failure and shock) for two
hospitals in 1994-95. The first line in the chart summarises
the result for Hospital 1, using the data from a large
number of hospitals as the reference group.

Hospital 1 has a mean RIDIT of 0.56 with a 95%

confidence interval from 0.53 to 0.6; this indicates that
the distribution of LOS for this hospital is higher than for
the reference group. Also in the figure, Hospital 2 is
compared with the reference data. Here the mean RIDIT is
0.49 with a 95% confidence interval from 0.44 to 0.53. Thus
LOS for Hospital 2 has a similar distribution to the
reference data.

The third line in the figure shows a comparison of

67
Hospital 1 and Hospital 2, with Hospital 2 being the
reference group. Relative to Hospital 2, Hospital 1 has a
mean RIDIT of 0.58, with a 95% confidence interval of 0.53
to 0.63. Thus Hospital 1 has longer lengths of stay than
Hospital 2, as already demonstrated for proportions
exceeding the median value, and shown in Figure 3. However,
the RIDIT method is likely to be more powerful for detecting
such differences because it takes into account the entire
LOS distributions. The method in Figure 3 only examines the
proportion of cases exceeding 5 days, which was the median
LOS for this DRG in the reference database.

The R commands to produce Figure 9 are as follows.

> h<-c(.56,.48,.58)
> u<-c(.6,.53,.63)
> l<-c(.53,.44,.53)
> x<-c(1,2,3)
> k<-c("Hosp 1 v Ref","","Hosp 2 v Ref","","Hosp1 v Hosp2")
> plot(h,ylim=c(.3,.7),axes=F,main="Mean RIDITs and confidence limits
for\ndifferences between reference values and 2 hospitals for
ANDRG252.",xlab="",ylab="")
> box()
> axis(side=1,tick=F,labels=k)
> axis(side=2)
> arrows(x,u,x,l,angle=90,code=3,col="red").

In the next section of this chapter we describe another

method for studying LOS and similar data that invloves the
Kaplan-Meier Life Table and the log rank test. These methods
are a powerful alternative for analysing time-to-event data

In summary, this section has addressed methods for

dealing with categorical data, especially those in ordered
categories. The following areas have been covered -
(1) RIDITS, the C statistic, AUC and Somers’delta.
(2) The rank sum test
(3) Stratification and direct standardisation
(4) Graphical methods.

The following section deals with numerical data.

Section 4
Numerical Data
Introduction.
In this section we examine methods for numerical data.
Most numerical QI data have severely skewed distributions
and numerous tied values making conventional methods of
analysis difficult. We concentrate on the use of medians,

68
percentiles, and graphical methods, and for time-to-event
data like LOS, Life Table methods. To obtain percentile
values, we use the RIDIT methods described in the previous
section of this chapter.
It is important to distinguish between methods that
provide management with information for planning services,
and methods that are useful to help workers improve the
quality of their systems and processes of care. Poor
management planning will have an important effect on the
quality of care and will increase common cause variation.
However, our chief interest in this book is with methods to
aid staff in improving systems and processes of care. For
example, it is important for management to have accurate
information about the average number of patients treated
within DRG’s so that resources can be deployed properly.
However, there is nothing that staff can (or should) do
about these numbers as it is their job to treat the patients
who present for care.

It is essential that any analysis convert data into

usable information and present it in a clear and simple way
to those who are striving to improve their work. Complex
methods that these workers lack a feel for are unlikely to
succeed.

LOS demonstrates the difficulties that can be

encountered and we shall use it to illustrate a general
approach that we believe should be satisfactory for QI work.
LOS is generally stratified by DRG with possible further
grouping by severity and co-morbidity. These data have
severe positive skew in almost all DRG’s with very long
right tails in their distributions and there are often large
outlier values present. Moreover they differ in their
distribution from DRG to DRG.

Thus, while in a few DRG’s an exponential or log-normal

distribution is approximated, for the majority of DRG’s this
is not the case. In addition, in some DRG’s with appreciable
mortality, the lengths of stay for survivors and those who
die may differ. For all these reasons, the usual methods for
analysing numerical data, which rely on approximations to
the normal distribution, are of limited value.

With very large samples, confidence intervals can be

calculated for the mean LOS using large sample normal theory
that should be approximately correct due to the Central
Limit Theorem. However, when a t -test is used to compare
LOS data from different hospitals, its power is often low
compared with the corresponding rank test. Nevertheless, it

69
is necessary to be aware that rank tests can give incorrect
results if the distributions in the groups being compared
have different shapes (Hart 2001).

Thompson and Barber (2000) show that rank tests can

mislead with some economic data. This occurs when the shapes
of the distributions of the data in the samples differ. When
this occurs, it is often more important to find out why than
to do significance tests. For example, there may be patients
in one of the groups that differ in their responses and
outcomes from others in that group and it is then more
important to identify the characteristics of the differing
subgroups than it is to do significance tests. These authors
recommend t-tests and bootstrapping. It is essential to
examine the distributions of the data in the two groups
being compared by using simple histograms and boxplots.

We believe that it is useful to employ percentiles to

analyse numerical surveillance data such as LOS. Approximate
percentile values are obtained using RIDITS. As well as the
data for an individual hospital, there will usually be
available data for a number of hospitals in a region or
state that can be used as a reference database. There are
two ways in which percentiles can be used and it is
necessary to have a clear idea of the differences in the two
approaches.

First, percentile values can be calculated for a DRG

from the reference database using RIDITs. These values can
then be applied to the data for an individual hospital and a
mean percentile (RIDIT) value for that institution can be
obtained. If the mean RIDIT value (Rbar) for a hospital is
in excess of 0.5 there may be a quality problem and further
systems analysis may be required.

It may also be useful to look at other percentile

values such as the 90th. If the 90th percentile value for the
reference database is 10 days, we would expect that about
10% of the patients in a hospital would have LOS values in
excess of 10 days for the relevant DRG. If the proportion
were much larger a quality problem might exist.

Secondly, the percentile values for the individual

institution can be obtained. For example, we may employ the
90th percentile for LOS for a particular DRG. In this case
we seek the numerical value at the 90th percentile for that
institution. If the 90th percentile for LOS is say 20 days
and that for the data in the reference database or from
another institution is only 15 days, there may be a quality
problem and an analysis of the relevant underlying systems

70
may be required.

When applying percentiles calculated from reference

database values to individual institutions, we are
interested in mean percentiles (RIDITs) or the proportions
of patients with numerical values above the relevant
percentiles. When using percentiles calculated from data
from individual institutions, we seek the corresponding
numerical values. In the latter case, the numerical values
from the reference database for those percentiles may also
be useful.

LOS data are time-to-event data and Life Table methods

may be useful for their analysis. This approach is
illustrated at the end of this section.

Using means and trimmed means.

Means and trimmed means have been used to describe LOS

data. The mean of a sample of LOS data for a particular DRG
may be useful for determining needed resources but its value
for QI studies is limited because of the distribution of the
data. This has led to the use of the trimmed mean, which
measures the mean of the data with outlier values excluded.
One way to achieve this is to exclude values greater than
the upper inner fence.

The upper inner fence is calculated by adding 1.5 times

the difference between the value of the 25th and 75th
percentiles (the interquartile range) to the 75th percentile
value using the reference data. The inner fence is thought
to be a good cutoff point for detecting outlier values.
However, an alternative is to use, for example, the
reference data 90th percentile in QI studies.

Reports of the mean and trimmed mean are usually

accompanied by the coefficient of variation, which is the
standard deviation divided by the mean. The coefficient of
variation serves as a measure of the relative variability of
the data. However, as Davies (1998) has pointed out, means
and standard deviations can be misleading with markedly
skewed data. An estimator that is derived from two
misleading estimators is itself unlikely to aid real
understanding. Thus the above estimates, while useful for
some administrative purposes, are limited in their ability
to describe data for QI.

Using medians and percentiles.

An alternative estimator is the median (Davies 1998),

71
which is the value at the 50th percentile. On its own it too
is of limited value when used with LOS data. One reason for
this is that there are often many tied values at the median
so that it can be very insensitive, and two sets of data can
differ substantially and still share a common median. To see
this, suppose that two hospitals each treated 100 patients
from a DRG with an expected median value of 5 days. In each
case suppose that 15% of the LOS values were tied at 5 days
but that the first hospital’s 5 days LOS values were between
its 36th and 51st lengths of stay, whereas the second
hospital’s values were between the 48th and 63rd lengths of
stay. Although the median value is the same for these
hospitals, they have LOS distributions with differing
locations.

We propose that it is reasonable to assume that the 15

values could be thought of as being spread uniformly over
the 5th day, that is between 4.5 and 5.5 days. Although this
is somewhat artificial, it allows the median to be estimated
within that period by linear interpolation. By doing this,
it is possible to use medians to distinguish between
hospitals like those in the above hypothetical example.

A convenient way to determine the position of the

percentiles of a distribution of ordered categorical data is
to calculate RIDIT scores that have been described in the
previous section of this chapter.

Example.

We use los data from the 1994-95 Australian national

diagnosis related group (ANDRG 252 heart failure and shock)
in a large central or reference database for a large number
of hospitals, and for two of its component hospitals. These
data for Hospital 1, Hospital 2 and all the hospitals in
tabular format are in [Link], [Link] and
[Link]. They are also shown in a single column in
[Link], [Link] and [Link]. We have
already referred to these data in section 3 of this chapter.

Using R.

The following can be used to obtain the LOS values in a

single column if these data are available with numbers of
patients grouped by LOS –

> load(“[Link]”) # Function for performing [Link] from a text

file or (in Windows) the clipboard. The clipboard has been used here;
the columns LOS and NUMBER have been copied from [Link].
> [Link]()

72
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
# The data are imported into a [Link] called datain.
> a<-0;b<-0
> for (i in 1:length(datain[,1])){a<-rep(datain[i,1],datain[i,2]);b<-
c(b,a)}
> b<-b[-1] # Removes the first data value in b that is 0.
> [Link](b,file="[Link]",sep=",",
[Link]=F,[Link]="LOSH2")

These data may be available as a single column and it

may then be required to group the numbers of patients by
LOS. The following R commands will achieve this –

> bb<-table(b) # The column of LOS values for Hospital 2 are in the
vector b.
> bc<-dimnames(bb)
> bc<-[Link](bc$b)
> ba<-[Link](bb)
> a<-[Link](bc,ba)
> [Link](a,file="[Link]",sep=",",[Link]=F,
[Link]=c("LOS","NUMBER"))

Proceed as follows to find the median value, the values

for the 25th (lower quartile), 75th (upper quartile), and
preferably also the 90th percentile for the data from both
the reference database and the individual hospitals.

Using R –

> load("[Link]")
> [Link]() # [Link]
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> median(datain[,1])
[1] 5
> quantile(datain[,1],.25)
25%
3
> quantile(datain[,1],.75)
75%
9
> quantile(datain[,1],.9)
90%
16

A convenient way to determine the position of the

percentiles of a distribution of ordered categorical data is
to calculate RIDIT scores that have been described in the
previous section of this chapter.

Using R –

> [Link]() # [Link]

73
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> b<-datain[,1]
> bb<-table(b)
> bc<-dimnames(bb)
> LOS<-[Link](bc$b)
> NUMBER<-[Link](bb)
> a<-[Link](LOS,NUMBER)
> A<-a[,2]/2
> B<-0
> B[2]<-a[1,2]
> for(i in 2:(length(a[,2])-1)){B[i+1]<-B[i]+a[i,2]}
> D<-A+B
> RIDIT<-D/sum(a[,2])
> DRG252RIDIT<-[Link](a,A,B,D,RIDIT)
> [Link](DRG252RIDIT,file="[Link]", sep=",",[Link]=F)
> DRG252RIDIT[DRG252RIDIT[,1]==5,]
bc ba A B D RIDIT
5 5 379 189.5 1772 1961.5 0.4692584

> length(datain[,1][datain[,1]==5])
[1] 379
> length(datain[,1][datain[,1]<5])
[1] 1772
> length(datain[,1][datain[,1]<6])
[1] 2151
> (1772+.5)/4180
[1] 0.4240431
> (2150+.5)/4180
[1] 0.5144737
> 4.5+(5.5-4.5)*(0.5-0.424)/(0.5145-0.424)
[1] 5.339779

It can be seen that, for the reference data, the median

value for the reference database is 5. From [Link]
the corresponding RIDIT score is seen to be 0.4693. There
were 379 patients who had an LOS of 5 days and because of
these large numbers of tied values, there is no RIDIT score
precisely equal to 0.5. However, the number of patients with
LOS values less than 5 days was 1772. Thus if the first
value at 5 days were considered to be unique, its RIDIT
score would be (1772+.5)/4180=0.424, the RIDIT score being
the number less than the value in question plus half of
those at that value, all divided by the total number of
values. The last value for 5 days was case number 2151. The
RIDIT score for this case would be (2150+.5)/4180=0.5145
(assuming it also to be unique). Thus, the median value lies
in this interval, and it is therefore 5 days. However,
assuming that the 379 lengths of stay of 5 days are spread
uniformly in the interval 4.5 to 5.5 days, we can by linear
interpolation obtain a value which would enable it to be
differentiated from other hospital LOS data which also have
a median of 5 days. The interpolation formula is
LOSM=LOS1+(LOSN-LOS1)(RM-R1)/( RN-R1), where the subscripts N,

74
M, and 1 refer to the last of the tied values for 5 days,
the median value, and the first of the tied values
respectively, and RM=0.5 is the RIDIT score for the median.
In addition, R1 and RN refer to the RIDIT scores for the
first and last category. Thus
LOSM=4.5+(5.5-4.5)(0.5-0.424)/(0.5145-0.424)=5.3 days.

Table 8 shows the relevant values for the median and

selected percentiles, first assuming all values for a day
are tied, and second assuming that an interpolated figure is
appropriate (shown in brackets). Note that, although the
method has been illustrated by LOS data, it should be
equally useful for other similar data such as that for cost
of stay and waiting times. Note also that, for the
interpolation method, one uses the percentiles calculated
from the data for the institution concerned. One does not
use RIDIT scores calculated from a reference database that
are used to calculate mean RIDITS for the data from the
individual hospitals. The latter are discussed in the
section on rank tests and mean RIDITs below.

Using boxplots.

A useful way to compare LOS and similarly distributed

data graphically is to employ boxplots that are
conventionally summary graphs of the minimum, maximum and
25th, 50th and 75th percentile values for the distribution
being displayed. The boxplots in Figure 10A, from
[Link] are difficult to read because of the large
outlier values. The upper inner fences are all less than 25
so Figure 4B shows the data for LOS up to 25. It can be seen
in Figure 4B that Hospital 1 has higher values and Hospital
2 similar values to the reference data for LOS in ANDRG 252.

Using R –

> [Link]() # los252refa

Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> hosp<-rep(0,length(datain[,1]))
> hospbox<-cbind(hosp,datain)

> [Link]() # los252h1a

Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> hosp<-rep(1,length(datain[,1]))
> hospbox1<-cbind(hosp,datain)

> [Link]() # los252h2a

Loading data.

75
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> hosp<-rep(2,length(datain[,1]))
> hospbox2<-cbind(hosp,datain)

> a<-c(hospbox[,1],hospbox1[,1],hospbox2[,1])
> b<-c(hospbox[,2],hospbox1[,2],hospbox2[,2])
> hospbox0<-[Link](a,b)
> hospbox0[1,]
a b
1 0 3
> boxplot(hospbox0)
> boxplot(hospbox0[,2]~hospbox0[,1])
> [Link](hospbox0,file="[Link]", sep=",",[Link]=F)

# Data in hospbox0 are in 2 columns, the first is 0 for the reference

data, 1 for Hospital 1 and 2 for Hospital 2.
> hospbox0[,1]<-hospbox0[,1]+1 # Add 1 so the first column values are 1,
2 and 3
> l<-c("","Reference","","Hospital 1","","Hospital 2","")
> boxplot(hospbox0[,2]~hospbox0[,1],ylab="Days",main="Boxplots of LOS\n
for Reference data and Hospitals 1 & 2",axes=F)
> box()
> axis(side=c(2,3,4))
> axis(side=1,tick=F,labels=l) # Figure 10(a)

> a<-hospbox0[,2]<=25
> data1<-hospbox0[a,]
> boxplot(data1[,2]~data1[,1],ylab="Days",main="Boxplots of LOS\n for
Reference data and Hospitals 1 & 2",axes=F)
> box()
> axis(side=c(2,3,4))
> axis(side=1,tick=F,labels=l) # Figure 10(b)

Example (continued).

Although the inner fence is conventionally used, the

90th percentile for the data for each hospital is a more
valuable alternative. It is a useful guide to which data
values are outliers and, in addition, it excludes a known
proportion of that hospital’s data. Figure 11 is a boxplot
employing the 90th percentile. It can be seen that, for
hospital 1, 10% of lengths of stay exceed 21 days, and for
hospital 2, 10% of lengths of stay are above 18 days.

Using R –

> All<-c(1,3,5,9,16) # Figure 11

> H1<-c(1,4,7,12,21)
> H2<-c(1,3,5,8,18)
> l<-c("","Reference","","Hospital 1","","Hospital 2","")
> boxplot(All,H1,H2,ylab="Days",main="Boxplots of LOS to 90th
percentiles\n for Reference data and Hospitals 1 & 2",axes=F)
> box()
> axis(side=c(2,3,4))
> axis(side=1,tick=F,labels=l)

76
Using relative frequencies.

For QI work, it is desirable that methods that

summarise the entire data distribution be available, and
these will frequently be graphical.

Example (continued).

Figure 12 shows the LOS data for ANDRG 252. It is a

comparison of the reference database values with those for
Hospital 1 and Hospital 2. The LOS values are on the
horizontal axis and the proportion of patients for each LOS
value are on the vertical axis. Hospital 1 has less short
lengths of stay (1-7 days) and more in the range from 8 to
30 days than the reference values while Hospital 2 has
slightly more lengths of stay around five days and less
around 8 to 14 days than the reference values. However,
Hospital 1 has relatively less lengths of stay in the 2 to 6
days range and more at higher values than Hospital 2. The
statistical significance of these differences and their
practical importance can be assessed by the rank sum test
and mean RIDITS, and these are described below.

Using R –

> load(“[Link]”)
> [Link]()# [Link]
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> d0<-datain # Reference data
> [Link]() # [Link]
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> d1<-datain # Hospital 1 data
> [Link]() # [Link]
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> d2<-datain # Hospital 2 data
> d0[,2]<-d0[,2]/sum(d0[,2])
> d1[,2]<-d1[,2]/sum(d1[,2])
> d2[,2]<-d2[,2]/sum(d2[,2])
> d0[,1][d0[,1]>30]<-30
> d1[,1][d1[,1]>30]<-30
> d2[,1][d2[,1]>30]<-30
> x<-sum(d0[,2][d0[,1]==30])
> e0<-d0[d0[,1]<30,]
> e0[length(e0[,1])+1,1]<-30
> e0[length(e0[,1]),2]<-x
> x<-sum(d1[,2][d1[,1]==30])
> e1<-d1[d1[,1]<30,]
> e1[length(e1[,1])+1,1]<-30
> e1[length(e1[,1]),2]<-x

77
> x<-sum(d2[,2][d2[,1]==30])
> e2<-d2[d2[,1]<30,]
> e2[length(e2[,1])+1,1]<-29 # No 29; therefore need a zero before the
30
> e2[length(e2[,1]),2]<-0
> e2[length(e2[,1])+1,1]<-30
> e2[length(e2[,1]),2]<-x
> plot(e0[,2]~e0[,1],type="l",lwd=2,col="blue",lty=1, ylim=c(0,
max(c(e0[,2],e1[,2],e2[,2]))),ylab="Proportion",xlab="LOS, Reference
solid blue, Hospital 1 dashed red, Hospital 2 dotted
green.",main="Distribution of LOS for Reference data and\n for Hospital
1 and Hospital 2 data, DRG 252.")
> lines(e1[,2]~e1[,1],lwd=2,col="red",lty=2)
> lines(e2[,2]~e2[,1],lwd=2,col="dark green",lty=3)
# Figure 12.

Although plotting Figure 12 is not difficult, the

re-arrangement of the data in preparation is fiddling and
time-consuming. Others wishing to use this chart to
illustrate the distributions of LOS values may prefer to
perform the data preparation in a spreadsheet (experienced R
users may know of simpler ways to produce similar charts).

Although there are substantial differences between

Hospital 1 and Hospital 2, it should be noted that factors
other than the efficiency of the hospital are important for
this DRG. For example, it has a large number of elderly
patients admitted in acute pulmonary oedema, many of whom
come from nursing homes. Such patients may settle quickly
with appropriate treatment and be ready for discharge to
their nursing homes in a couple of days. However, when the
patient is ready to leave hospital, a nursing home bed may
be unavailable or home help may need to be arranged, and
this may delay discharge. Thus the ability of a hospital to
discharge these patients may be determined largely by
factors such as the capacity of community and family
resources to receive them, although careful discharge
planning is likely to make a considerable difference. In
addition, hospitals with specialist services and teaching
roles may have different lengths of stay from community
hospitals for reasons that are not related to quality or
efficiency. Another factor is that patients in a particular
DRG may not be homogeneous. For example, patients in ANDRG
252 in the early 1990’s at a town with a large retirement
population often had chronic pulmonary disease as well as
cardiac disease whereas in a city hospital there were many
more patients with hypertensive renal disease.

Using transformations.

Another approach is to use a transformation. The marked

positive skew of these data suggests that a logarithmic

78
transformation may be appropriate. Unfortunately, the
majority of QI and other staff working in hospitals do not
find that the logarithm of days conveys to them very useful
information.

An alternative transformation that we have already

referred to and that has great potential is the percentile
transformation. LOS data are aggregated for large hospital
systems, for example the data for all the hospitals in a
state or country. As we have shown, it is easy to rank all
the data for a particular DRG for a particular year and to
assign an approximate percentile value to each LOS using
RIDITS. The only limitation is that there must be sufficient
data available in the reference database for the DRG in
question; it is not possible to get stable values for the
upper percentiles with less than 200 subjects and 500 is
preferable. For example, with a database containing 200
records, the 97.5th percentile would have only 5 values
above it so that it may be unduly influenced by outlier
values and therefore may be unreliable. However, this
difficulty is not a serious one as those DRG’s that must be
excluded represent relatively uncommon disorders. The RIDIT
scores for ANDRG 252 in 1994-95 for the reference data and
for two hospitals are shown in the files [Link],
[Link] and [Link].

Using RIDITS and percentiles.

Using R –
> hosps<-c(.36,.45,.5,.56,.68) # To produce Figure 13
> p90<-.62
> HospX<-.6
> boxplot(hosps,main="RIDITs for all hospitals ANDRG252\nand hospital
X",ylab="Mean RIDITs",xlab="Boxplot with 90th percentile.")
> points(HospX,pch="+",font=2,cex=1.25)
> points(p90,pch="-",font=2,cex=2)
> legend(locator(n=1),legend=c("Study Hospital","90th
percentile"),pch=c("+","-"),cex=1.25)

Having calculated percentile values using RIDIT scores

for lengths of stay for a particular DRG from a reference
database, these values can then be allocated to the lengths
of stay data for the hospital being studied. The mean RIDIT
can then be calculated for that hospital’s data to give an
estimate of its mean percentile value. If the hospital from
which these data come is inefficient or if it has problem
patients referred to it, for example cardiac patients
referred to a specialist cardiology center, the mean RIDIT
may be larger than for the reference data, for which the
mean RIDIT is always 0.5. An advantage of this approach is

79
that percentiles are widely used and understood by hospital
staff. These rank percentile values are also valuable for
constructing control charts, which are the subject of
Chapter 4. In addition, as is shown in Figure 7, the mean
RIDITs for a DRG for a group of hospitals can be displayed
in a boxplot with the value for an individual hospital of
interest marked separately.

Example (continued).

[Link] shows LOS values for Hospital 1 and 2.

The R function twosample() can be used to study the two
groups.

Using R –

> load(“[Link]”)
> twosample()
Loading data. # Data from clipboard using [Link]
# Column 1 has 1 for hospital 1 and 2 for hospital 2
# Column 2 has LOS values
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y

[Link] [Link] tests [Link] [Link] [Link] 1

First sample
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 4.00 7.00 10.09 12.00 158.00

Second sample
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 3.000 5.000 8.377 8.000 97.000

[Link] [Link] tests [Link] [Link] [Link] 2

Welch Two Sample t-test, t = 1.41, DF = 366, P-value = 0.16, LR = 3.

First mean = 10.094, second mean = 8.377, difference = 1.718.
95% confidence interval for the difference is -0.679 to 4.114.
95% confidence limits for first mean are 8.551 to 11.637.
95% confidence limits for second mean are 6.534 to 10.219.

Rank sum test Z = 2.78, p-value = 0.00536, LR = 48.

Variance = 6935508, concordances = 18052, discordances = 25396.
MEAN RIDIT = 0.57898, 95% confidence limits are 0.52455 to 0.63341.
RIDIT variance 0.0007711566.

Median test p-value 0.00416, LR = 61.

Fligner-Killeen dispersion test p-value = 0.0018, LR = 130.

The t-test suggests that there is no difference between

the 2 distributions and the rank sum test suggests a
substantial difference. The dispersion test shows that the
two distributions have differing dispersions. However, their

80
distributions are otherwise similar (Figure 14). The median
test suggests that there is a difference between their
medians of 5 and 7. The RIDIT difference is 0.579 with a 95%
confidence interval from 0.525 to 0.633.

[Link] [Link] tests [Link] [Link] [Link] 3

# Figure 14(a)
Enter a name for the boxplot charts eg Boxplots for group 1 and group 2
Figure 8A. Boxplots for hospital 1 and hospital 2. # Boxplot is Figure
14(a).

[Link] [Link] tests [Link] [Link] [Link] 4

# Figure 14(b)
Enter a heading for the first histogram Figure 8B. Hospital 1 LOS in
days.

Enter a heading for the second histogram Hospital 2 LOS in days.

Enter the number of breaks 12

[Link] [Link] tests [Link] [Link] [Link] 5

Using proportions within quartiles.

Example (continued).

Figure 15 shows the proportions of LOS values that are

in the 4 quartiles for the 1994-95 ANDRG 252 data using the
quartile values from the reference data (1st quartile data
to the 25th percentile, 2nd from above 25th to 50th
percentile, 3rd from above 50th to 75th and 4th from above 75th
to 100th percentile). It is clear that Hospital 1 has a
large proportion of values in the 4th quartile and Hospital
2 is similar to the reference data.

Using R –
> descriptive()
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y

1. Summary 2. Boxplot 3. Histogram 4. QQ plot 5. [Link] 6. Quit 1

Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 3.000 5.000 8.349 9.000 297.000
Number = 4180, SD = 12.946.
95% Confidence interval for the mean is 7.955982 to 8.741147.

1. Summary 2. Boxplot 3. Histogram 4. QQ plot 5. [Link] 6. Quit 6

Using R continued -

> load("[Link]") # Figure 15

> [Link]()
Loading data. # Loading from [Link] via clipboard
Data from clipboard (C) or file (F) c

81
Do data column(s) have heading(s) (Y/N) y
> a<-datain[,1]
> [Link]()
Loading data. # Loading from [Link] via clipboard
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> b<-datain[,1]
> x<-c(length(a[a<=3]),length(a[a>3 & a<=5]),length(a[a>5 &
a<=9]),length(a[a>9]))
> y<-c(length(b[b<=3]),length(b[b>3 & b<=5]),length(b[b>5 &
b<=9]),length(b[b>9]))
> x<-x/sum(x)
> y<-y/sum(y)
> l<-c("1st","2nd","3rd","4th")
> m1<-"Hospital 1 data DRG 252\nby quarters of reference distribution"
> m2<-"Hospital 2 data DRG 252\nby quarters of reference distribution"
> mat<-matrix(1:2,1,2)
> layout(mat)
> barplot(x,main=m1,col=c("red","orange","yellow","brown"))
> axis(side=1,labels=l)
> barplot(y,main=m2,col=c("red","orange","yellow","brown"))
> axis(side=1,labels=l)
> mat<-matrix(1:1,1,1)
> layout(mat)

Using the proportion exceeding median and 90th percentile.

A very useful approach is to determine the proportion

of values for a DRG and hospital that exceed the median
value calculated from the reference database. In this case
the grouped reference median value is used, for example 5
days for ANDRG 252, and not the interpolated value for the
individual hospitals. Although this approach appears to be
inefficient because it discards the information about the
individual LOS values, it is in fact surprisingly effective.

Example (continued).

For the ANDRG 252 data and Hospital 1, the proportion

exceeding the median is 0.59 with 95% confidence interval
0.53 to 0.65 and the p-value for its difference from the
reference value of 0.4854 is 0.0004. The reference value is
0.4854 rather than 0.5 because of tied values at the median
in the reference database. Figure 3 shows the values for
Hospital 1 and Hospital 2, and the difference between them.

Using R –

> [Link]()
Loading data. # [Link]
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> q<-length(datain[,1][datain[,1]>5]) # 5 is the reference median
> q
[1] 170

82
> q1<-length(datain[,1])
> q1
[1] 287
> [Link]()
Loading data. # [Link]
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> d<-length(datain[,1][datain[,1]>5])
> d
[1] 2029
> d1<-length(datain[,1])
> d1
[1] 4180
> d/d1
[1] 0.4854067
> load("[Link]")
> proportion()
Enter numerator 170
Enter denominator 287
Is a reference proportion available? y
Enter reference proportion .4854
Proportion = 0.592, lower 95% limit = 0.533, upper limit = 0.65.
Mid-P 95% limits are 0.535 and 0.648.
P = 0.00036, LR = 588, mid-p value = 0.00029 , LR = 712.

In addition to examining the proportions exceeding the

median, the proportions exceeding a higher percentile such
as the 90th of the reference data can be valuable for
studying outlier values. This approach is likely to be most
useful in a control chart and will be discussed further in
Chapter 4.

Using Rank tests and Mean RIDITS.

The frequency of tied values at many lengths of stay,

together with the skewness of their distributions within
DRGs, make methods for ordered categorical data, including
the rank sum test, attractive for analysing these data.
However, these methods are not distribution free. The one
sample or signed rank test usually cannot be employed
because it is sensitive to severe asymmetry in the data, a
prominent characteristic of LOS data. If a one sample test
is required with these data, it will usually be necessary to
use the binomial test, the single sample test for
proportions described in Section 1 of this chapter.

In addition, the two sample or rank sum test, also

equivalently called Kendall’s S test and the Mann Whitney
test requires that the shapes of the distributions of the
samples being compared are similar. Two samples from the
same DRG usually will have similar shaped distributions.

In spite of its interpretation as an average percentile

rank, the mean RIDIT or equivalent C statistic may appear

83
inscrutable to some QI staff. However, suppose a mean RIDIT
of say 0.6 is obtained in a comparison between two groups of
LOS values. This means that if a single LOS value is sampled
randomly from each group there is a 2×(0.6-0.5)=20% greater
chance that the larger value comes from the group with the
distribution of values shifted towards the higher lengths of
stay than the other group. This version is called Somers’
delta. It was described in Section 3 of this chapter. As
already noted, it is a generalisation of the difference
between two proportions to data in ordered categories. In
addition, when the mean RIDIT is 0.6 and one value is
sampled at random from each group, the odds are 0.6 divided
by 0.4 or 3 to 2 that the larger value of the two comes from
the group with the distribution of values shifted towards
the higher lengths of stay.

Another interpretation of the mean RIDIT as the C

statistic is commonly used in logistic regression. Imagine
two hospitals where one had all its lengths of stay larger
than the other. In this case the lengths of stay would
“discriminate” perfectly between two hospitals and LOS could
be used to determine which hospital an individual patient
came from; the mean RIDIT or C statistic would be one. Now
suppose that the lengths of stay for the two hospitals were
identical, so that LOS would be useless for discriminating
between the hospitals; in this case the mean RIDIT would be
0.5.

In practice, there is usually some separation with there

being many LOS values common to each hospital but with one
having a number of values higher than the other. The mean
RIDIT is a measure of the strength of this separation. When
the mean RIDIT is thought of in this way, it can be seen
that it provides useful information about the magnitude of
the difference between two groups.

As we have already described, the C statistic is often

used in conjunction with a Receiver Operating Characteristic
(ROC) chart; it represents the area under the ROC curve
(AUC). ROC charts do not seen to be as useful for
summarising LOS data as the mean RIDIT. ROC charts will be
described in Chapter 6 that deals with logistic regression.

Figure 13 shows how a mean RIDIT value for a particular

hospital can be compared with the mean RIDIT values for a
group of hospitals using a boxplot. Suppose that mean RIDITs
for a particular DRG are calculated for each hospital
supplying data to a large central database and that, for a
particular DRG, they range from 0.36 to 0.68 with median

84
value of 0.5, upper and lower quartiles of 0.56 and 0.45
respectively, and 90th percentile 0.62. Figure 13 shows the
result for Hospital X whose mean RIDIT=0.6. Note that, in
this case, we employ the mean RIDIT for each hospital using
RIDIT scores derived from the reference database values.

Using R –

# Calculates RIDIT scores for Hospital 1 using reference data values

> [Link]() # [Link]
Loading data.
Data from clipboard (C) or file (F) c
> d<-datain[,1]
> [Link]() # [Link]
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> a<-datain[,1]
> for (i in 1:length(a)){if (length(which(d==a[i]))==0){d<-c(d,a[i])}}
> dd<-table(d);dc<-dimnames(dd);dc<-[Link](dc$d);da<-[Link](dd)
> a0<-da/2;b0<-0;b0[2]<-da[1]
> for(i in 2:(length(da)-1)){b0[i+1]<-b0[i]+da[i]}
> d0<-a0+b0;ri<-d0/sum(da)
> b<-NA
> for (i in 1:length(a)){q<-which(dc==a[i]);b[i]<-ri[q]}
> datain<-[Link](a,b)
> m<-mean(datain[,2])
> for (i in 1:length(a)){if (length(which(d==a[i]))==0){d<-c(d,a[i])}}
> dd<-table(d);dc<-dimnames(dd);dc<-[Link](dc$d);da<-[Link](dd)
> a0<-da/2;b0<-0;b0[2]<-da[1]
> for(i in 2:(length(da)-1)){b0[i+1]<-b0[i]+da[i]}
> d0<-a0+b0;ri<-d0/sum(da)
> b<-NA
> for (i in 1:length(a)){q<-which(dc==a[i]);b[i]<-ri[q]}
> o<-order(a);a<-a[o];b<-b[o]
> ridlos252h1a<-[Link](a,b)
> x<-[Link](table(a))
> y<-[Link](table(b))
> RIDIT<-[Link]([Link](y[,1]))
> LOS<-x[,1];NUMBER<-x[,2]
> ridlos252h1b<-[Link](LOS,RIDIT,NUMBER)
> m
[1] 0.5662692
# Mean RIDIT is in m, RIDIT values for Hospital 1 using reference data
# are in the variables riditlos252h1a and riditlos252h1b. The former has
# data sorted, the latter has them in a table.

Example (continued).

Using the ordered categorical data methods, a

comparison of the reference database LOS data for ANDRG 252
with that from Hospital 1 gives a mean RIDIT for the latter
of 0.57. The 95% confidence interval is 0.53 to 0.60. The
mean RIDIT value shows that a randomly selected value from
the Hospital 1 data has a 2×(0.57-0.5)=14% greater chance of
being larger than a randomly selected value from the

85
reference data. This is Somer’s delta described in Section
3; it is a generalisation of the difference between two
proportions to ordered categorical data. The corresponding
odds are 4 to 3. The rank sum test value is Z=3.77,
p-value<0.0001. For the Hospital 2 and reference data
comparison using the ANDRG 252 data, the mean RIDIT is 0.49
(95% confidence interval 0.44 to 0.53, Z=0.53, p-value>0.2).
For the comparison of Hospital 1 and Hospital 2, and using
the Hospital 2 data as the reference, the mean RIDIT=0.58
(95% confidence interval 0.53 to 0.63, Z=2.7,
p-value=0.007). Somers’ delta is 2×(0.58-0.5)=16%, so that
there is a 16% greater chance that a randomly selected
patient from Hospital 1 would have a longer LOS than a
similarly selected patient from Hospital 2. These
differences are summarised in Figure 3. The corresponding
odds are approximately 4 to 3.

Combining data.

Frequently it will be desirable to combine data from

several DRG’s. For example, the staff of a medical unit may
wish to obtain an overall view of LOS for patients in the
major medical groups or a surgical unit may have a similar
need for their major abdominal of other surgical groups. In
addition, as severity and co-morbidity groupings appear
within DRG’s, or more DRG’s are created to cater for
severity and co-morbidity, this need may increase.

Since each LOS RIDIT score is an LOS value transformed

to a percentile rank, the staff of a medical unit could add
RIDIT scores, derived from reference database values, for
patients in major medical DRGs. They could then calculate a
mean RIDIT value by dividing by the number N of LOS values
used as described in the Section 3 of this chapter. Since
the reference mean RIDIT value is always 0.5 and the
approximate variance is 1/12N, this is equivalent to using
indirect standardisation. However, combining data from DRGs
having very different data distributions should be avoided.
For example, eye surgery patients often have very short LOS
values and to attempt to combine them with a DRG like ANDRG
252 Heart failure and shock would not produce a sensible
result.

The method of direct standardisation can also be

applied to LOS exceeding the mean RIDIT. Although mean LOS
values can also be analysed in this way, using for example a
bootstrap method to obtain confidence limits, the resulting
summaries may not be very illuminating to QI staff. However,
overall or weighted average lengths of stay may be required

86
for resource planning tasks.

Suitable standardising weights for direct

standardisation are the numbers in the relevant categories
in the reference database, scaled so as to sum to one.
However, there is a major limitation. Some high volume DRG’s
such as cataract surgery normally have very short lengths of
stay and very large numbers are tied at the median. These
should be analysed separately and not combined with data
from other longer stay DRG’s usually representing more
complex medical problems.

Example.

The 1994-95 LOS data for ANDRG 249 (Circulatory

disorders with Acute Myocardial Infarction without invasive
cardiac procedure uncomplicated) for the reference data and
the two hospitals are in [Link], [Link] and
[Link]. The percentile values are shown in Table 9.
These data are summarised in Figure 16. It is similar to
Figure 10, the corresponding figure for the ANDRG 252 data.

Using R -

> ref<-c(11.4,8.6,6.5,4.5,1) # Figure 16

> h1<-c(14.4,10.6,8.3,5.5,1)
> h2<-c(9,6.9,5.7,4.6,1)
> l<-c("","Reference","","Hospital 1","","Hospital 2","")
> boxplot(ref,h1,h2,ylab="Days",main="Boxplots of LOS to 90th
percentiles\nfor reference data and hospitals 1 and 2.",axes=F)
> box()
> axis(side=c(2,3,4))
> axis(side=1,tick=F,labels=l)

Figure 17 corresponds for ANDRG 249 to Figure 12 for

the ANDRG 252 data; these charts display the proportion of
patients for each LOS. The figure shows a comparison of
Hospital 1 and Hospital 2 with the reference data. The
difference between the two hospitals is greater with these
data than with the ANDRG 252 data. The data are in 2 columns
the files [Link], [Link] and [Link],
column 1 has the LOS values and column 2 the count of
patients with that LOS.

Figure 18 shows the proportions of LOS values in the

four quartiles obtained from the reference data. Hospital 1
has more in the 4th quartile and Hospital 2 has very few LOS
values in either the third or fourth quartiles.

Figure 19 shows the proportions exceeding the median

value derived from the reference database for Hospitals 1

87
and 2, and for the difference between them. This difference
is more marked than with ANDRG 252. The corresponding chart
for the latter is Figure 8. Figure 20 shows the result when
the data for ANDRG 249 and ANDRG 252 are combined. There are
clearly large differences between the hospitals, and these
differences are similar for both DRG’s.

Using R –

> [Link]() # Figure 17

Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> d0<-datain # [Link]
> [Link]()
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> d1<-data # [Link]
> [Link]()
Loading data.
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> d2<-datain # [Link]
> d0[,2]<-d0[,2]/sum(d0[,2])
> d1[,2]<-d1[,2]/sum(d1[,2])
> d2[,2]<-d2[,2]/sum(d2[,2])
> d0[,1][d0[,1]>30]<-30
> d1[,1][d1[,1]>30]<-30
> d2[,1][d2[,1]>30]<-30
> x<-sum(d0[,2][d0[,1]==30])
> e0<-d0[d0[,1]<30,]
> e0[length(e0[,1])+1,1]<-30
> e0[length(e0[,1]),2]<-x
> x<-sum(d1[,2][d1[,1]==30])
> e1<-d1[d1[,1]<30,]
> e1[length(e1[,1])+1,1]<-30
> e1[length(e1[,1]),2]<-x
> x<-sum(d2[,2][d2[,1]==30])
> e2<-d2[d2[,1]<30,]
> e2[length(e2[,1])+1,1]<-29 # No 29; therefore need a zero before the
30
> e2[length(e2[,1]),2]<-0
> e2[length(e2[,1])+1,1]<-30
> e2[length(e2[,1]),2]<-x
> plot(e0[,2]~e0[,1],type="l",lwd=2,col="blue",lty=1, ylim=c(0,
max(c(e0[,2],e1[,2],e2[,2]))),ylab="Proportion",xlab="LOS, Reference
solid blue, Hospital 1 dashed red, Hospital 2 dotted
green.",main="Distribution of LOS for Reference data and\n for Hospital
1 and Hospital 2 data, DRG 252.")
> lines(e1[,2]~e1[,1],lwd=2,col="red",lty=2)
> lines(e2[,2]~e2[,1],lwd=2,col="dark green",lty=3)

In R (Figure 18) –

> descriptive() # Use function descriptive to get quartiles

Loading data. # [Link]
Data from clipboard (C) or file (F) c

88
Do data column(s) have heading(s) (Y/N) y

1. Summary 2. Boxplot 3. Histogram 4. QQ plot 5. [Link] 6. Quit 1

Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 4.000 6.000 7.129 9.000 96.000
Number = 2853, SD = 5.43.
95% Confidence interval for the mean is 6.929649 to 7.328325.

1. Summary 2. Boxplot 3. Histogram 4. QQ plot 5. [Link] 6. Quit 6

> load("[Link]")
> [Link]()
Loading data. # Loading from [Link] via clipboard
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> a<-datain[,1]
> [Link]()
Loading data. # Loading from [Link] via clipboard
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> b<-datain[,1]
> x<-c(length(a[a<=4]),length(a[a>4 & a<=6]),length(a[a>6 &
a<=9]),length(a[a>9]))
> y<-c(length(b[b<=4]),length(b[b>4 & b<=6]),length(b[b>6 &
b<=9]),length(b[b>9]))
> x<-x/sum(x)
> y<-y/sum(y)
> l<-c("1st","2nd","3rd","4th")
> m1<-"Hospital 1 data DRG 249\nby quarters of reference distribution"
> m2<-"Hospital 2 data DRG 249\nby quarters of reference distribution"
> mat<-matrix(1:2,1,2)
> layout(mat)
> barplot(x,main=m1,col=c("red","orange","yellow","brown"))
> axis(side=1,labels=l)
> barplot(y,main=m2,col=c("red","orange","yellow","brown"))
> axis(side=1,labels=l)
> mat<-matrix(1:1,1,1)
> layout(mat)

Using R cont –

> a<-c(.58,.66,.73) # Figure 19

> b<-c(.23,.3,.37)
> d<-c(.26,.36,.45)
> x<-c(a[1],b[1],d[1])
> y<-c(a[2],b[2],d[2])
> z<-c(a[3],b[3],d[3])
> x1<-c(1,1,1)
> y1<-c(2,2,2)
> z1<-c(3,3,3)
> l<-c("Hosp 1","","Hosp 2","","Diff")
> plot(x,ylim=c(0,1),type="p",pch="-",cex=2,axes=F,xlab="",
ylab="Proportion exceeding reference median",col="red")
> points(y,pch=16,col="red")
> points(z,pch="-",cex=2,col="red")
> lines(a~x1,col="red")
> lines(b~y1,col="red")
> lines(d~z1,col="red")
> box()
> axis(side=1,tick=F,labels=l)

89
> axis(side=2)
> title(main="Proportion exceeding reference median ANDRG 249\nand
difference between hospitals 1 & 2 with confidence intervals.")

> a<-c(.57,.62,.66) # Figure 20

> b<-c(.33,.38,.43)
> d<-c(.18,.25,.32)
> x<-c(a[1],b[1],d[1])
> y<-c(a[2],b[2],d[2])
> z<-c(a[3],b[3],d[3])
> x1<-c(1,1,1)
> y1<-c(2,2,2)
> z1<-c(3,3,3)
> l<-c("Hosp 1","","Hosp 2","","Diff")
> plot(x,ylim=c(0,1),type="p",pch="-",cex=2,axes=F,xlab="",
ylab="Proportion exceeding reference median",col="red")
> points(y,pch=16,col="red")
> points(z,pch="-",cex=2,col="red")
> lines(a~x1,col="red")
> lines(b~y1,col="red")
> lines(d~z1,col="red")
> box()
> axis(side=1,tick=F,labels=l)
> axis(side=2)
> title(main="Proportion exceeding reference median ANDRG 249 & ANDRG
252\ncombined and difference between hospitals 1 & 2 with confidence
intervals.")

> h<-c(.62,.43,.69) # Figure 21

> u<-c(.67,.46,.75)
> l<-c(.58,.4,.64)
> x<-c(1,2,3)
> k<-c("Hosp 1 v Ref","","Hosp 2 v Ref","","Hosp1 v Hosp2")
> plot(h,ylim=c(.3,.8),axes=F,main="Mean RIDITs and confidence limits
for\ndifferences between reference values and 2 hospitals for
ANDRG249.",xlab="",ylab="")
> box()
> axis(side=1,tick=F,labels=k)
> axis(side=2)
> arrows(x,u,x,l,angle=90,code=3,col="red")
# Figure 20 and Figure 21 are produced differently. Users could decide
which they prefer.

Figure 21 for ANDRG 249 corresponds to Figure 9 for the

ANDRG 252 data. The mean RIDIT values for the two hospitals
demonstrate that the average rank percentile value for the
Hospital 1 data is 0.62, and for hospital 2 it is 0.43. The
mean RIDIT for Hospital 1 using Hospital 2 as the reference
is 0.69. This means that, relative to Hospital 2 with an
adjusted average rank percentile value of 0.5, the Hospital
1 mean rank percentile value is 0.69. Although the latter
has been calculated directly, it may be found as follows -
Rbar=0.62-0.43+0.5=0.69. Finally, Figure 22 shows mean RIDIT
values when the data for ANDRG 249 and ANDRG 252 are
combined.

Table 10 summarises the results of the calculation of

90
mean RIDITs for both sets of data (ANDRG 249 and the
previously described ANDRG 252), first for each DRG
separately, and then for both combined. Tables 11 and 12
summarise the results of the calculations for the
proportions exceeding the median of the reference values,
using the methods of Section 1 of this chapter.

The standardising weights used for the combined results

are proportional to the numbers in the reference database in
each DRG scaled so that they sum to one. For ANDRG 249 the
number was 2853 and for ANDRG 252 it was 4180. Thus the
first weight is 0.406 and the second weight is 0.594. H1 and
H2 refer to the two hospitals, R is the reference database
and NNV is the non-null RIDIT variance.

In R –

> rrh1<-.406*.62+.594*.57
> vrh1<-.406^2*.000526+.594^2*.000318
> urh1<-rrh1+1.96*vrh1^.5
> lrh1<-rrh1-1.96*vrh1^.5

> rrh2<-.406*.43+.594*.49
> vrh2<-.406^2*.000297+.594^2*.0005
> urh2<-rrh2+1.96*vrh2^.5
> lrh2<-rrh2-1.96*vrh2^.5

> rh2h1<-.406*.69+.594*.58
> vh2h1<-.406^2*.000838+.594^2*.000771
> uh2h1<-rh2h1+1.96*vh2h1^.5
> lh2h1<-rh2h1-1.96*vh2h1^.5

> wrh1<-.406*(292875-170065)+.594*(641244-482292)
> vrh1<-.406^2*499349914+.594^2*1775149374
> wrh1/vrh1^.5

> wrh2<-.406*(202305-276722)+.594*(304326-321028)
> vrh2<-.406^2*538193403+.594^2*973706520
> wrh2/vrh2^.5
> 2*(1-pnorm(1.93))

> wh2h1<-.406*(21534-8813)+.594*(25396-18052)
> vh2h1<-.406^2*3939259+.594^2*6935508
> wh2h1/vh2h1^.5

The differences are calculated using the R function

twoproportions() described in Section 1 of this chapter. The
combined results have been calculated with the Mantel-
Haenszel method in the R function mhci() (Section 1). The
relevant proportions are for DRG 249 REF 1408/2853 H1
115/175 H2 56/188 and for DRG 252 REF 2029/4180 H1 170/287
H2 70/162.

The need to understand when to use reference database

91
RIDIT scores and scores derived from an individual
hospital’s data is a potential source of confusion. To
calculate mean RIDITs for each hospital, to determine the
proportions of LOS values for each hospital within reference
data percentile categories, or to use RIDIT scores for
indirect standardisation, use the reference database RIDIT
scores. To calculate percentile values for individual
hospitals so that the corresponding numerical values may be
found, use the data for the hospitals concerned.

Survival methods.

Survival analysis may be useful for analysing ime-to-

event data like LOS and wating times. We illustrate with
[Link] that has two columns, the first having a 1
for Hospital 1 and a 2 for Hospital 2. The second column
contains the LOS data. Figure 17(a) is the survival curve
fot the 2 hospitals; due to outlier values, it is difficult
to see any difference. Figure 17(b) shows the same data but
with LOS longer than 50 days truncated to 50 days. The log
rank test p-value with the data before truncation is 0.024,
and the median LOS values are 7 days for Hospital 1 and 5
days for Hospital 2. The difference is similar to but not as
marked as that using the rank sum test. Setting rho=1 in
survdiff() changes the test to the Peto & Peto modification
of the Gehan-Wilcoxon test that may give more weight to the
more frequent shorter LOS values. This gives a result
similar to the rank sum test.

Using R –

> library(survival)
> load("[Link]")
> [Link]()
Loading data. # [Link]
Data from clipboard (C) or file (F) c
Do data column(s) have heading(s) (Y/N) y
> su<-survfit(Surv(LOS)~HOSP,data=datain)
> su
Call: survfit(formula = Surv(LOS) ~ Gp, data = datain)

n events median 0.95LCL 0.95UCL

Gp=1 287 287 7 6 7
Gp=2 162 162 5 4 6
>
plot(survfit(Surv(LOS)~HOSP,data=datain),col=c("blue","green"),main="LOS
data for Hospital 1 (blue) and Hospital 2
(green)",xlab="LOS",ylab="Probability not discharged",lwd=2) # Figure
23(a)
> survdiff(Surv(LOS)~HOSP,data=datain)
Call:
survdiff(formula = Surv(LOS) ~ Gp, data = datain)

92
 N Observed Expected (O-E)^2/E (O-E)^2/V
Gp=1 287 287 308 1.41 5.11
Gp=2 162 162 141 3.08 5.11

Chisq= 5.1 on 1 degrees of freedom, p= 0.0238

> datain[,2][datain[,2]>50]<-50
>
plot(survfit(Surv(LOS)~HOSP,data=datain),col=c("blue","green"),main="LOS
data for Hospital 1 (blue) and Hospital 2 (green)",xlab="LOS truncated
at 50 days",ylab="Probability not discharged",lwd=2) # Figure 23(b),
LOS>50 days truncated to 50 days to make survival curve easier to view.
# Peto & Peto modification of the Gehan-Wilcoxon test.
> svd<-survdiff(Surv(LOS)~HOSP,data=datain,rho=1)
> svd
Call:
survdiff(formula = Surv(LOS) ~ HOSP, data = datain, rho = 1)

N Observed Expected (O-E)^2/E (O-E)^2/V

Gp=1 287 144.1 160.4 1.67 8.01
Gp=2 162 95.1 78.8 3.40 8.01

Chisq= 8 on 1 degrees of freedom, p= 0.00465 Peto & Peto modification

of the Gehan-Wilcoxon test.

This section has described methods for dealing with

numerical data that arise in hospital epidemiology work.
These data are frequently not normally distributed. Methods
based on medians, ranks and percentiles are often most
useful for analysing these data.

Section 5
Miscellaneous Topics
Making Decisions (Decision Analysis)

It is often important to determine whether some actions

such as employing screening tests should be incorporated
into a process. For example, should otherwise well patients
undergoing relatively minor surgery such as hernia repair or
excision of skin lesions be required to have a preoperative
chest X-ray? Some people may consider such investigations
important for medico-legal reasons.

Suppose that a preoperative test has a sensitivity of

95% and a specificity of 95%. Thus if the patient has a
lesion there would be 95 chances in 100 that it will be
detected. Similarly, if there is no lesion, there are 95
chances in 100 that the test will be negative. Suppose also
that, in this group of patients, the probability of a lesion
being present is 1 in 200.

93
 From these data we can say that the probability of the
lesion is 0.005 and the probability of a positive test in
its presence is 0.005×0.95=0.00475. Similarly the
probability of no lesion is 0.995 and the probability of a
positive test in the absence of a lesion is
0.995×(1-0.95)=0.04975. From these figures the probability
of a positive test is 0.04975+0.00475=0.05425 and the
probability of a lesion when the test is positive is
0.00475/0.05425=0.088. Thus only 0.088×100≈9 in every 100
positive tests (about 1 in 11) will signify the presence of
a definite lesion and the remaining 91 will ultimately be
shown to be false positive results.

These calculations, called Bayes’ Rule, illustrate the

importance of the prevalence of the condition in question in
determining the usefulness of the test. Thus if 200 tests
are taken there will be about 11 positive results
(0.05425×200), only 1 of which would represent the presence
of a lesion and the remainder would be false positives. In
addition, one negative result is likely to occur in the
presence of the lesion (false negative result) in every 20
groups of 200 tests performed
0.00025/(0.94525+0.00025)≈1/4000).

Although such testing is not necessarily wrong, it is

important to understand what it costs and also what further
costs are generated in investigating the false positives to
establish that no lesion exists. Often these investigations
can be quite invasive and there may be risk of harm,
including unfounded anxiety and inconvenience, as well as
added cost.

Furthermore, it should be determined what the result

would be if this money were spent differently. For example,
if it were used to improve IM systems and processes it might
be possible to prevent 3 or 4 hospital acquired infections
that might be more serious and costly than the one lesion
that the test finds.

Two further concepts of importance are the disutility

of missing the lesion and the sensitivity of the decision
analysis. If the lesion were to cause death in the
postoperative period it would be much more serious than if
it simply had a low probability of causing a minor
postoperative complication. In the former case its
disutility would be measured at 100% (or its utility zero),
and in the latter case its disutility might be estimated as
10%. Thornton and Lilford (1995) show how to incorporate
disutility, or its opposite utility, in a decision analysis.

94
 The probabilities of disease presence that are used in
the decision analysis may not be known precisely and it may
be necessary to estimate them. Because of this it is usual
to select a range of probabilities and to determine the
outcome of the analysis for each of them. By doing this it
is possible to judge which course of action would be
appropriate in a particular set of circumstances. For
example, the sensitivity analysis may indicate that the same
course of action is appropriate over a wide range of
probabilities of disease presence, and this would add to
confidence in the analysis. Alternatively, there may be a
cut off point where each course of action is equally
appropriate. Thus on one side of this threshold, for example
when the probability of disease is below a certain level,
the test would not be indicated, whereas if it were above
this level, the test should be performed.

Assessing agreement.

QI work often involves staff reading hospital records

to extract evidence of adverse outcomes. In addition, IM
staff must make judgements about whether or not there is a
surgical site infection or other nosocomial infection. It is
important that there are carefully standardised definitions
of the adverse occurrences and other events to be detected,
and that the people doing the assessment are able to agree
about their presence or absence.

Agreement has conventionally been assessed by the Kappa

statistic (Fleiss 1981). However, the proportion of
agreement for both normal and abnormal assessments may be
preferable (Grant 1991), although it is probable that the
exercise of performing a formal assessment of agreement is
much more important than the statistic used to estimate the
ability of observers to agree. We illustrate the method
using the hypothetical data in Grant’s paper shown in Table
13.

For these data the proportion of agreement for

abnormality is e/(e+f+g)=16/22=0.74, and the proportion of
agreement for normality is h/(h+f+g)=28/34=0.82. Approximate
confidence intervals for these proportions can be calculated
using proportion() as described in Section 1 of this
chapter. The 95% limits are 0.5 to 0.89 and 0.65 to 0.93
respectively.

Grant (1991) has emphasized the importance of obtaining

samples of sufficient size and, if the proportion of
abnormal results occurring naturally is small, selecting

95
samples for testing that contain a higher proportion of
abnormal findings. This would be easy to do when
documentation is to be assessed but more difficult with an
infrequently occurring nosocomial infection. Grant (1991)
has stated that the proportions of agreement that are
considered satisfactory will differ depending on the
subject, but that a confidence interval including 50%
agreement will generally be inadequate provided samples of
sufficient size are employed. The method can easily be
extended to the assessment of multiple observers as
described in Grant’s (1991) paper.

An important aspect is the assessment of observer bias,

which is illustrated in Table 14. In this case there are the
same number of disagreements as before but they all occur
when Observer A indicates a normal and Observer B an
abnormal result. The probability of getting such a result by
chance is 2/26=0.032 for a two tailed test. Since this is
unlikely, it is probable that the two observers are making
systematically different assessments of patients who may be
borderline, and that further training is required to improve
overall assessment capability.

Example.

Two medical students undertook as their Social and

Preventative Medicine project a study of postoperative
respiratory complications (McGrath and Morton 1986). This
involved searching patients’ files for evidence of such a
complication. In order to assure themselves that they were
doing this properly they first performed an independent
retrospective assessment of 82 patient files to determine
whether a respiratory complication existed. The result of
their study is shown in Table 15.

For these data the proportion of agreement for

normality was 0.94 (95% confidence interval 0.85 to 0.98)
and the proportion of agreement for abnormality was 0.82
(95% confidence interval 0.60 to 0.95). They concluded that
their level of agreement was satisfactory and thus were able
to proceeded with their study.

References.

Abramson J and Gahlinger P “PEPI Version 4" 2000

([Link]

Altman D, Machin D, Bryant T, and Gardner M “Statistics with

Confidence” 2nd edition British Medical Journal London 2000.

96
Antelman G “Elementary Bayesian Statistics” Cheltenham
Edward Elgar 1997.

Bissell A “A Negative Binomial Model for Varying Element

Sizes” Biometrika 1972;59:435.

Bross I “How to Use RIDIT Analysis” Biometrics 1958;14:19.

Clayton D and Hills M “Statistical Models in Epidemiology”

Oxford Science Publications Oxford 1993.

Dalgaard P “Introductory Statistics with R” New York

Springer 2002.

Davies H "Informative Presentation of Summary Data" Hospital

Medicine 1998;59:154.

Davies H, Crombie I, and Tavakoli M “When Can Odds Ratios

Mislead” British Medical Journal 1998;316:989.

Fleiss J “Statistical Methods for Rates and Proportions”

John Wiley and Sons New York 1981.

Gart J and Nam J “Approximate Interval Estimation of the

Difference in Binomial Parameters: Correction for Skewness
and Extension to Multiple Tables” Biometrics 1990;46:637.

Gibberd R, Pathmeswaran A, and Burtenshaw K “Using Clinical

Indicators to Identify Areas for QI” Journal of Quality in
Clinical Practice 2000;20:136.

Glynn R and Buring J “Ways of Measuring Rates of Recurrent

Events” British Medical Journal 1996;312:364.

Goodman S “Toward Evidence-Based Medical Statistics 2: The

Bayes Factor” Annals of Internal Medicine 1999;130:1005.

Graham P, Mengersen K and Morton A “Confidence Limits for

the Ratio of Two Rates Based on Likelihood Scores:
Non-iterative Method” Statistics in Medicine 2003;22:2085.

Grant J “The Fetal Heart Rate is Normal, Isn’t It? Observer

Agreement of Categorical Assessments” The Lancet
1991;337:215.

Greenland S and Robins J “Estimation of a Common Effect

Parameter from Sparse Follow-up Data” Biometrics 1985;41:55.

Hanley J and McNeil B “The Meaning and Use of the Area under
a Receiver Operating Characteristic (ROC) Curve” Diagnostic

97
Radiology 1982;143:29.

Hart A “Mann-Whitney Test is not just a Test of Medians:

Differences in Spread are Important” British Medical Journal
2001;323:391.

Hart M, Lee K, Hart R and Robertson J “Application of

Attribute Control Charts to Risk-Adjusted Data for
Monitoring and Improving Health Care Performance” Quality
Management in Health Care 2003;12:5.

Hofer T, Bernstein S, Hayward R, and DeMonner S “Validating

Quality Indicators for Hospital Care” The Joint Commission
Journal on Quality Improvement 1997;23:455.

Horan T and Culver D “Methods for Comparing Surgical Site

Infection Rates” in APIC Infection Control and Applied
Epidemiology edited by Olmsted R 1996 St Louis Mosby-Year
Book Inc.

Ihaka R and Gentleman R. “R: A Language for Data Analysis

and Graphics” Journal of Computational and Graphical
Statistics 1996;5:299.

Kantor S, Winkelstein W, and Ibrahim M “A Note on the

Interpretation of the RIDIT as a Quantile Rank” American
Journal of Epidemiology 1968;87:609.

Kirkwood B and Sterne J “Essential Medical Statistics” 2nd

ed Oxford Blackwell Science 2003.

McGrath G and Morton A “A Study of Postoperative Respiratory

Complications” Unpublished Social and Preventative Medicine
Project, University of Queensland 1986.

Maindonald J and Braun J “Data Analysis and Graphics”

Cambridge University Press 2003.

Merrer J, Santoli F, Vecchi A, Tran B, Jonghe B and Outin H

“Colonisation Pressure and Risk of Acquisition of
Metyhicillin Resistant Staphylococcus Aureus in a Medical
Intensive Care Unit” Infection Control and Hospital
Epidemiology 2000;21:718-723.

Miettinen O “Theoretical Epidemiology” John Wiley and Sons

New York 1985.

Myatt M “Open Source Solutions–R” ([Link]/epitools/).

Nam J “Confidence Limits for the Ratio of Two Binomial

98
Proportions Based on Likelihood Scores: Non-iterative
Method” Biometrics Journal 1995;37:375.

Newcombe R “Interval Estimation for the Difference Between

Independent Proportions” Statistics in Medicine 1998;17:873.

Newcombe R “Logit Confidence Intervals and the Inverse Sinh

Transformation” The American Statistician 2001;55:200.

Reynolds H “The Analysis of Cross-classifications” The Free

Press New York 1977.

Rothman K and Boice J “Epidemiologic Analysis with a

Programmable Calculator” 2nd ed. 1982 Newton MA:
Epidemiology Resources.

Sackett D, Strauss S, Richardson W, Rosenberg W, and Haynes

R "Evidence-based Medicine How to Practice and Teach EBM"
2nd edition New York Churchill Livingstone 2000.

Salemi C, Morgan J, Kelleghan S, and Hiebert-Crape B

“Severity of Illness Classification for Infection Control
Departments: A Study in Nosocomial Pneumonia” American
Journal of Infection Control 1993;21:117.

Selvin S “A Further Note on the Interpretation of RIDIT

Analysis” American Journal of Epidemiology 1977;105:16.

Sato T “On the Variance Estimator for the Mantel-Haenszel

Risk Difference” Biometrics 1989;45:1323.

Spiegelhalter D “Funnel Plots for Comparing Institutional

Performance” Statistics in Medicine 2005;24:1185.

Thompson S and Barber J “How Should Cost Data in Pragmatic

Randomised Controls be Analysed” British Medical Journal
2000;320:1197.

Thornton J and Lilford R “Decision Analysis for Medical

Managers” British Medical Journal 1995;310:791.

Woodworth G “Biostatistics A Bayesiian Introduction” New

York John Wiley and Sons 2004.

Verzani J “Using R for Introductory Statistics” Chapman &

Hall/CRC 2004.

Waller J, Addy C, Jackson K and Garrison C “Confidence

Intervals for Weighted Proportions” Statistics in Medicine

99
1994;13:1071.

Dean A, Arner T, Sangam S, Sunki G, Friedman R, Lantinga M,

Zubieta J, Sullivan K, and Smith D “Epi Info 2000, a
Database and Statistics Program for Public Health
Professionals for Usa on Windows 95, 98, NT, and 2000
Computers” Centers for Disease Control and Prevention,
Atlanta, Georgia, USA, 2000

van der Tweel I “Repeated Looks at Accumulating Data: To

Correct or not to Correct” European Journal of Epidemiology
2005;20:205.

100