funnel: Funnel plot for benchmarking health units

Description

Produces a variety of funnel plots comparing health units or ICUs (intensive care units) making easy to identify those units which deviate from the group. There is a function that calculates all the values required and returns the values for all units and the funnel, and there is a function that calls graphical parameters from the former values. The options of funnels available are the funnel for rate, for ratio of rates, for proportions, for difference of proportions and for ratio of proportions.

The funnel for rates are usually plots of either SMR or SRU at vertical axis. If the direct method is chosen, the horizontal axis will display the number of admissions. If the indirect method is chosen instead, the expected number of deaths will be displayed for SMR or the expected length of stay for SRU. As consequence of this differentiation, the interpretation regarding the clssification of the points displayed will be the same in every case.

The funnel for ratio of rates are usually plots of ratios of SMRs (or SRUs) within the same units. These two SMRs are, for example, from the same units in different time periods. Therefore, it expresses how the SMR changed over time. If the number of expected deaths is different in both periods, the plot will return at the horizontal axis a parametrization of the geometric mean of the expected number of deaths in both periods for each unit. If the number of expected deaths is identical in both periods, the plot will return at the horizontal axis the arithmetic mean of the observed number of deaths in both periods for each unit.

The funnel for proportions plots on the vertical axis the percentage of observed deaths of the units and on the horizontal axis the number (volume) of admissions. The funnel for ratio of proportions and for difference of proportions are usually used to express the fraction of deaths of the same units in different time period. Therefore, they express how the fraction of deaths changed over time in each unit. If one picks the difference of proportions, the horizontal axis will display a parametrization of the arithmetic mean of the number of admissions in both periods. If one picks the ratio of proportions, the horizontal axis will display a parametrization of the geometric mean of the number of admissions in both periods.

Usage

funnel(
  unit,
  y,
  n,
  n1,
  n2,
  o,
  o1,
  o2,
  e,
  e1,
  e2,
  lambda1 = sum(o1)/sum(n1),
  lambda2 = sum(o2)/sum(n2),
  pi1 = sum(o1)/sum(n1),
  pi2 = sum(o2)/sum(n2),
  y.type = c("SMR", "SRU"),
  p = c(0.95, 0.998),
  theta,
  method = c("normal", "exact"),
  direct = FALSE,
  myunits = rep(0, length(unit)),
  option = c("rate", "ratioRates", "prop", "diffProp", "ratioProp"),
  printUnits = TRUE,
  plot = TRUE,
  digits = 5,
  overdispersion = FALSE,
  ...
)
# S3 method for funnel
print(x, ...)
# S3 method for funnel
plot(
  x,
  ...,
  col = c("darkblue", "paleturquoise3", "gray26"),
  lwd = 2,
  lty = c(2, 6, 1),
  bty = "n",
  pch = 21,
  pt.col = "white",
  bg = "orange",
  pt.cex = 1.5,
  auto.legend = TRUE,
  text.cex = 0.7,
  text.pos = NULL,
  mypts.col = "darkblue",
  printUnits = x$printUnits,
  xlab = x$xlab,
  ylab = x$ylab,
  xlim = x$xlim,
  ylim = x$ylim
)
rateFunnel(
  unit,
  y,
  n,
  o,
  e,
  y.type,
  p,
  theta = 1,
  method = c("exact", "normal"),
  direct,
  ...,
  printUnits,
  auto.xlab = TRUE,
  xlab = c("Volume of cases", "Expected values"),
  ylab = y.type[1],
  xlim = c(0, max(rho)),
  ylim = c(min(lowerCI[[which(p == max(p))]]), max(upperCI[[which(p == max(p))]])),
  myunits,
  digits,
  overdispersion
)
changeRateFunnel(
  unit,
  n1,
  n2,
  o1,
  e1,
  o2,
  e2,
  lambda1,
  lambda2,
  y.type,
  p,
  ...,
  printUnits,
  auto.xlab = TRUE,
  xlab = c("Average observed count", "Expectation per period"),
  auto.ylab = TRUE,
  ylab = c(paste0(y.type[1], "'s Ratio"), paste0("Log(", y.type[1], "'s Ratio)")),
  ylim = c(max(lowerCI[[which(p == max(p))]]) - 1.5 * theta, min(upperCI[[which(p ==
    max(p))]]) + 1.5 * theta),
  xlim = c(0, max(rho)),
  myunits,
  digits,
  overdispersion
)
propFunnel(
  unit,
  o,
  n,
  theta,
  p,
  method = c("exact", "normal"),
  ...,
  printUnits,
  ylab = "%",
  xlab = "Volume",
  ylim = c(0, min(upperCI[[which(p == max(p))]]) + 2.5 * theta),
  xlim = c(0, max(n)),
  myunits,
  digits,
  overdispersion
)
changePropFunnel(
  unit,
  o1,
  o2,
  n1,
  n2,
  p,
  pi1,
  pi2,
  method = c("diff", "ratio"),
  ...,
  printUnits,
  xlab = "Sample size per period",
  auto.ylab = TRUE,
  ylab = c("Proportions difference", "Proportions ratio log"),
  ylim = c(max(lowerCI[[which(p == max(p))]]) - 6 * theta, min(upperCI[[which(p ==
    max(p))]]) + 6 * theta),
  xlim = c(0, max(rho)),
  myunits,
  digits,
  overdispersion
)

Arguments

unit

A factor vector representing the unit names.

A numeric vector representing the "Standardized rate" for each unit, usually the SMR (Standardized Mortality Ratio), or possibly the SRU (Standardized Resource Use), accordind to y.type . It's also called "indicator".

A numeric vector representing the case volume, or number of admissions, for each unit.

n1, n2

If one picks option = "ratioRates" or option = "diffProp" or option = "ratioProp", then n1 and n2 are numeric vectors representing the total of admissions at 1st and 2nd periods, respectively.

A numeric vector representing the observed death. Acceptable values are 0 (absence) or 1 (presence).

o1, o2

If one picks option = "ratioRates" or option = "diffProp" or option = "ratioProp", then o1 and o2 are numeric vectors representing the observed deaths at 1st and 2nd periods, respectively.

Used only when option = "rate" and direct = FALSE. This is a numeric vector representing the expected number of deaths.

e1, e2

If one picks option = "ratioRates", e1 and e2 are numeric vectors representing the expected number of deaths at 1st and 2nd periods, respectively.

lambda1, lambda2

Values correponding to the rate at which a death occurs in the institutions at the 1st and 2nd periods, respectively. It is assumed that the parameters o1 and o2 are distributed as oi ~ Poisson(lambdai) when option = "ratioRates". The default value for lambdai is lambdai = sum(oi)/sum(ni), where ni is the value of the parameter n1 or n2 when i equals 1 or 2. lambdai is the estimate for the mean of the poisson distribution.

pi1, pi2

Values correponding to the probability for the occurrence of a death in the institutions at the 1st and 2nd periods, respectively. Its assumed that the paramenters o1 and o2 are distributed as oi ~ Bin(pii,ni) when option = "diffProp" or option = "ratioProp". ni is the value of the parameter n1 or n2 when i equals 1 or 2. The default value for pii is sum(oi)/sum(ni), the estimate for the mean of the poisson distribution.

y.type

A character vector representing the indicator type. It is used to name the vertical axis if option = "rate" or option = "ratioRate" and ignored otherwise. It usually is 'SMR' or 'SRU'.

A confidence level numeric vector. The function will return a confidence interval for each value in p. The default is 2 and 3 standard deviations (p = c(.95, 998)).

theta

The target value which specifies the desired expectation for institutions considered "in control". Used when option = "prop" or option = "rate". Usually, this function internally estimates a theta to represent a central tendency of the group. However, one may want to set a pre-specified value for theta to indicate a "baseline" parameter for comparison (e.g 1 for option = "rate" or .20 for option = "prop"). If this is the case, the horizontal line representing theta may not be centralized in the funnel or may be even outside the funnel, making the plot look unusual.

method

There are two kinds of approximations for the CI, as mentioned in direct parameter. The one from the exact distribuition (binomial or poisson) and the one from the normal distribution. So, method is a character vector representing the kind of approximation desired, being "exact" (default) or "normal" the two options. It is used when option = "rate" or option = "prop". The original report makes no formal comparison of which method is best, however it is mentioned that the funnels from different methods should look identical or very similar if all units have 100 or more observations. If any unit has less, the funnel from the normal approximation may mislead the interpretation. See details.

direct

Logical (default = FALSE); Used when option = "rate". If TRUE, we assume the vector of rates "y" is being reported as a rate per (say) 1000 individuals, and that it has been transformed to a proportion between 0 and 1. The associated error (horizontal axis) will be measured accordingly to the size of the populations n. The CI - confidence interval - is calculated by a binomial distribution. If FALSE, the associated error will be measured accordingly to the expected number of deaths e. The CI is calculated by a poisson distribution instead. See details.

myunits

A numeric vector coded with 0 and 1 indicating which units one would like to benchmark among all units. These will be highlighted with dots of different collors in the plot.

option

A character specifying the type of funnel plot one wants to produce. It can assume "rate", "ratioRates", "prop", "diffProp" or "ratioProp". If option = "rate", funnel plots a standardized rate y versus the expected number of deaths or case volume (number of unit admissions) for all units. If option = "ratioRate", funnel can be used to compare units at two diferent periods. It plots a ratio of rates y versus a precision parameter rho. If option = "prop", funnel plots a proportion y versus its case volume (number of admissions). If option = "ratioProp" or option = "diffProp", funnel can be used to compare units at two diferent periods. It plots a ratio (or difference) of proportions y versus a precision parameter rho. See details.

printUnits

Logical (default = TRUE); If TRUE, the units are identified in the plot and printed in de console. The numbers plotted correspond to the row numbers printed in the console.

plot

Logical; If TRUE (default), the correspondent graphic is plotted with the standard options.

digits

Integer indicating the number of decimals to be used in the output.

overdispersion

Logical (default = FALSE); If TRUE, introduces an multiplicative over-dispersion factor phi that will inflate the CI null variance. See details.

...

Further arguments passed to plot.

An object of class 'funnel'.

col

A character vector representing the colors for the CI funnel lines. Must have same length of p + 1 with the target line color in the last position.

lwd

A positive number specifying the lines width. It's the same for all lines in the plot. See par.

lty

A numeric vector representing the CI lines types. See par.

bty

A character string which represents the type of box which is drawn around plots. See par.

pch

Either an integer or a single character specifying a symbol to be used as the default in plotting points. See points for possible values and their interpretation. Note that only integers and single-character strings can be set as a graphics parameter (and not NA nor NULL).

pt.col

A character specifying the points colors.

A character specifying the color to be used for the points background when pch = 21 (default). See par.

pt.cex

A numerical value giving the amount by which plotting points should be magnified relative to the default. See par.

auto.legend

Logical; If TRUE (default), prints a legend with default arguments.

text.cex

A numerical value giving the amount by which plotting text should be magnified relative to the default. See par.

text.pos

A position specifier for numbers that correspond to the units in the plot. Values of 1, 2, 3 and 4, respectively indicate positions below, to the left of, above and to the right of the points.

mypts.col

A character representing the color used to benchmark the units specified in myunits.

xlab, ylab

A title for the x and y axis. See title

xlim, ylim

Limits of horizontal and vertical axis. These limits are defined in the funnel plot and passed to plot.funnel. The user may redefine the limits in plot.funnel. Ultimately, these arguments are passed to plot.default.

auto.xlab, auto.ylab

Logical. If TRUE, one is not able to change x and y axis labels, respectively.

Value

A table with unit names, y, observed (Obs), expected (Exp) and admissions (N) for each unit, a binary column showing which units one would like to highlight in the plot (myunits) and final columns show which units are out of control.

Details

For every possible value of option, if overdispersion = TRUE, the CI can be inflated by a overdispersion parameter phi. There is a test for overdispersion which inflates the funnel if it's necessary. An "Winsorized" over-dispersion parameter is estimated and is used to inflate the funnel limits if it is significantly greater than 1. The parameter phi is returned as an funnel object.
If option = "rate", funnel plots a standardized rate y versus the expected number of deaths or volume value for several units.

To choose the direct argument, one should pay attention if one wants to use a Direct or Indirect Standardized Rate. If direct, we assume the rate is reported as a rate per (say) 1000 individuals, then it is treated as a proportion. If indirect, we assume it is a cross-sectional data that leads to a standardized event ratio.

In many circumstances we can assume an exact or approximate normal distribution for the data. Using the method argument, one could choose between "exact" or "normal". For direct standardized rates, the exact distribuition is binomial and for indirect standardized rates, the exact distribuition is poisson. Assume rho is the precision parameter (volume, for direct rates; expected value, for indirect rates). The original report claims that, for rho > 100, the normal and exact curves almost coincide. So, one could perfectly use normal approximation if ones data parameter precision is greater than 100, in general.

The console warns if there are units with volume/expected value less than 100.

phi = (1/total) * sum((y - theta) ^ 2 * rho)/g(theta)

var(y|theta,rho) = (phi * g(theta))/rho
If option = "ratioRate", funnel can be used to compare units at two diferent periods. It plots a ratio of rates y versus a precision parameter rho.

Suppose we have two measures for each institution: O1; E1 in a baseline period and O2; E2 in a subsequent period, and we wish to assess the change in the underlying rate (SMR or SRU). We shall only consider the ratio of rates option. The exact method will automatically be applied if E1 = E2, and the indirect method, of normal approximations, otherwise. On this second method, for low (especially zero) counts the funnel function adds 0.5 to all parameters O and E in order to stabilize the estimates.

Y = (O1/E1)/(O2/E2) and the target theta = lambda2/lambda1.
When E1 = E2, y is plotted versus the average observed count (rho).

When E1 is different of E2, i.e., it is used normal approximation. It is convenient to work on a logarithmic scale so that log(theta) is a target for log(Y). Y is plotted versus a different rho depending on the chosen rate.
If option = "prop", funnel plots a proportion y versus its volume. It is used for cross-sectional data. Suppose in each institution that O events are observed out of a sample size of N:

The indicator is the observed proportion y = O/N

Assume N is the precision parameter (volume). Similarly to when option = "rate", for N > 100 the normal and exact curves almost coincide. So, one could perfectly use normal approximation on the parameter method if ones data parameter precision is greater than 100, in general.

phi = (1/total) * sum((y - theta) ^ 2 * N)/g(theta)

var(y|theta,N) = (phi * g(theta))/N
If option = "ratioProp" or option = "diffProp", funnel can be used to compare units at two diferent periods. It plots a ratio (or difference) of proportions y versus a precision parameter rho to assess the change in the underlying proportion from pi1 to pi2. Normal approximations are used throughout, and for low (especially zero) counts, the function adds 0.5 to all arguments r and 1 to all arguments n in order to stabilize the estimates.

In the case option = "diffProp", the indicator is Y = (O2/N2 - O1/N1) and theta = pi2 - pi1. If option = "ratioProp", the indicator is Y = (O2/N2)/(O1/N1) and theta = pi2/pi1. It is convenient to work on a logarithmic scale, so that log(theta) is a target for log(Y) in this case as well.

For these two parameter options, the precision parameter (plotted at horizontal axis) can be interpreted as approximately the sample size per period.

References

Spiegelhalter, David J. "Funnel plots for comparing institutional performance." Statistics in medicine 24.8 (2005): 1185-1202.

Examples

Run this code

# NOT RUN {
# Loading data
data(icu)

# Some edition
icu$Saps3DeathProbabilityStandardEquation <- icu$Saps3DeathProbabilityStandardEquation / 100
icu <- icu[-which(icu$Unit == "F"),]
icu$myunits <- ifelse(icu$Unit == "A",1,0) #my units
icu <- droplevels(icu)

# Getting the cross-sectional arguments to use in funnel
x <- SMR.table(data = icu, group.var = "Unit",
              obs.var = "UnitDischargeName", pred.var = "Saps3DeathProbabilityStandardEquation")
myunit_names <- unique(icu$Unit[which(icu$myunits == 1)])
x$myunits <- ifelse(x$Levels %in% myunit_names, 1,0)

# Analysis of proportions
f1 <- funnel(unit = x$Levels[-1], o = x[-1,]$Observed, theta = x$Observed[1] / x$N[1],
            n = x[-1,]$N, method = "exact", myunits = x$myunits[-1], option = "prop", plot = FALSE)
f1
plot(f1, main = "Cross-sectional proportions")

# To analyze rates (SMR)
f2 <- funnel(unit = x$Levels[-1], y = x[-1,]$SMR, method = "exact", direct = TRUE,
            theta = x$SMR[1], e = x[-1,]$Expected, n = x[-1,]$N, o = x[-1,]$Observed,
            option = "rate", plot = FALSE)
f2
plot(f2, main = "Cross-sectional rate (SMR)")

# Creating a variable containing month information about each admission
icu$month <- as.numeric(format(as.Date(icu$UnitAdmissionDateTime),"%m"))

# First quarter
dt1 <- icu[which(icu$month %in% c(1,2,3)),]

# Second quarter
dt2 <- icu[which(icu$month %in% c(4,5,6)),]

# Getting the two period arguments to use in funnel
z <- SMR.table(data = dt1, group.var = "Unit", obs.var = "UnitDischargeName",
              pred.var = "Saps3DeathProbabilityStandardEquation")
w <- SMR.table(data = dt2, group.var = "Unit", obs.var = "UnitDischargeName",
              pred.var = "Saps3DeathProbabilityStandardEquation")

z$myunits <- ifelse(z$Levels %in% myunit_names, 1,0)
w$myunits <- ifelse(w$Levels %in% myunit_names, 1,0)
# To analyze periods using ratio rates with e1 = e1
f3 <- funnel(unit = z$Levels[-1], n1 = z$N[-1], o1 = z$Observed[-1],
            e1 = z$Expected[-1],
            n2 = w$N[-1], o2 = w$Observed[-1], e2 = z$Expected[-1],
            myunits = z$myunits[-1], option = "ratioRates", plot = FALSE)
f3
plot(f3, main = "Ratio of SMRs of periods with same expectation of death")


# To analyze periods using ratio rates with e1 =! e1
f4 <- funnel(unit <- z$Levels[-1], n1 = z$N[-1], o1 = z$Observed[-1],
            e1 = z$Expected[-1], n2 = w$N[-1], o2 = w$Observed[-1], e2 = w$Expected[-1],
            option = "ratioRates", plot = FALSE)
f4
plot(f4, main = "Ratio of SMRs of periods with different expectation of death",
    ylim = c(-1.5,1.5), xlim = c(0,200))

# To analyze periods by difference in proportions
f5 <- funnel(unit <- z$Levels[-1], n1 = z$N[-1], o1 = z$Observed[-1],
            n2 = w$N[-1], o2 = w$Observed[-1], option = "diffProp", plot = FALSE)
f5
plot(f5, main = "Difference in proportions of death for two periods")

# To analyze periods by ratio of proportions
f6 <- funnel(unit <- z$Levels[-1], n1 = z$N[-1], o1 = z$Observed[-1],
            n2 = w$N[-1], o2 = w$Observed[-1], option = "ratioProp", plot = FALSE)
f6
plot(f6, main = "Ratio of proportions of death for two periods")

rm(icu, x, z, w, dt1, dt2, unit, f1, f2, f3, f4, f5, f6)

# }

Run the code above in your browser using DataLab