*! version 1.0.2 12Feb2024

*! version 1.0.1 23Jan2024

*! version 1.0.0 28Oct2023

Title

funnelinst -- Funnel plot for comparing institutional performance

Syntax

funnelinst for a proportion:

funnelinst prop
num1
den1
[num2
den2]
[if] [in]
[, options]

funnelinst for a rate:

funnelinst rate
num1
den1
[num2
den2]
[if] [in]
[, options]

funnelinst for an indirect standardized rate:

funnelinst smr
num1
den1
[num2
den2]
[if] [in]
[, options]

funnelinst for a continuous variable:

funnelinst mean
num1
den1
[num2
den2]
[if] [in]
[, options]

In the syntax for funnelinst prop and rate, the numerator is the number of events and the denominator is the population size.

In the syntax for funnelinst smr, the numerator is the observed number of events and the denominator is the predicted number of events.

In the syntax for funnelinst mean, the numerator is the mean and the denominator is the standard error.

For all subcommands, when only two variables are specified (numerator1, denominator1), funnelinst analyzes the performance cross-sectionally.
When four variables are specified (numerator1, denominator1, numerator2, denominator2), funnelinst analyzes the change in the performance indicator
between two timepoints (longitudinal).

target(numlist) target rate with a maximum of two values; default is to use the computed overall average rate

chtype(string) change type for prop and rate subcommands. Available options are "diff" (difference), "ratio" (ratio),
and "or" (odds ratio)

npoints(#) number of points to use for constructing the funnel plot denominator; default is npoints(200)

xrangehigh(#) highest value to use in constructing the funnel plot denominator; default is 1% higher than the largest denominator value

logtrans log transformation of the performance measure Y (numerator/denominator); available for subcommands prop, rate and smr

arcsintrans arc sin transformation of the performance measure Y (numerator/denominator); available for subcommands prop and rate

sqrttrans square root tranformation of the performance measure Y (numerator/denominator); available for subcommand smr

pval(numlist) p-value(s) used as the basis for generating confidence intervals (up to a maximum of three); default is pval(0.001 0.025)

madjust adjusts the comparisonwise error rate to account for multiple testing; madjust overrides p-values specified in pval

exact perform exact estimation; available for subcommands prop, rate and smr when two variables are specified, and smr when four
variables are specified

winsor(#) percent of observations to be winsorized; default is no winsorizing winsor(0)

overdisp(string) adjust for overdispersion; options are additive overdisp(add) and multiplicative overdisp(mult)

ypercent multiplies performance measure by 100; available for subcommands prop and smr

ratedenom multiplies y-axis values by specified amount (e.g. 1000) and divides x-axis values by that specified amount; available for subcommand rate

ymin(#) presents observations and CI values in the figure above specified minimum

ymax(#) presents observations and CI values in the figure below specified maximum

saving(filename) specify the filename where the CIs, generated by funnelinst, are to be saved.
The suboption "replace" is allowed

figure[(twoway_options)] allows all twoway options to enhance the basic figure

Description

funnelinst generates funnel plots for comparing units (e.g. physicians, hospitals, etc.) on a given performance measure (e.g. readmission, mortality, etc.)
either cross-sectionally or longitudinally (i.e. the change between two timepoints), as described in Spiegelhalter (2005a, 2005b, 2012). The funnel plot may
be thought of as a control chart in which units that fall outside of the confidence limits be considered "out of statistical control". A particularly salient
feature of funnelinst is an option to account for over-dispersion of the performance measure due to unmeasured risk factors, which could lead to a large
number of units being inappropriately classified as out of control. funnelinst further offers adjustment of the comparisonwise error rate to account for
multiple testing using both the Bonferroni and False Discovery Rate methods (Jones et al. 2008).

Options

target(numlist) target rate for the performance indicator Y (numerator/denominator). When two values are specified (the maximum allowed),
funnelinst treats the target values as a range. The default is to use the computed overall average rate of the performance indicator.

chtype(string) type of analysis to evaluate change for prop and rate subcommands. Available options are "diff" (difference), "ratio" (ratio),
and "or" (odds ratio). Only one change type is available for subcommands smr and mean, therefore chtype(string) is not required
for those subcommands. chtype(string) can only be specified when two numerators and two denominators are specified.

npoints(#) specifies the number of points to use for constructing the funnel plot denominator (x-axis); default is npoints(200).

xrangehigh(#) specifies the highest value to use in constructing the funnel plot denominator (x-axis); default is 1% higher than the
largest denominator value.

logtrans log transforms the performance measure Y for computing control limits. However, the funnel plot axes are presented on the original scale.
logtrans is available for subcommands prop, rate and smr.

arcsintrans arc sin transforms the performance measure Y for computing control limits. However, the funnel plot axes are presented on the original scale.
arcsintrans is available for subcommands prop and rate.

sqrttrans square-root transforms the performance measure Y for computing control limits. However, the funnel plot axes are presented on the original scale.
sqrttrans is available for the subcommand smr.

pval(#) p-values to be used as the basis for generating confidence intervals. Up to three p-values may be specified. The default is pval(0.001 0.025)
which are roughly equivalent to 2 and 3 standard deviations.

madjust uses both the Bonferroni and False Discovery Rate (FDR) methods to adjust the comparisonwise error rate to account for multiple testing. The Bonferroni method
is extremely conservative and becomes even more so as sample size increases (given that the comparisonwise error rate = 0.05 / N). Conversely, the FDR method uses a step-down
approach which is likely to identify more outlier institutions (see Jones et al. 2008). An FDR of 5% means that among all institutions called significant, 5% of these are
truly null on average (Storey and Tibshirani, 2003). madjust will override any pval(#) that the user may have specified. The graph will present
Bonferroni and FDR adjusted 95% confidence intervals as well as a one-sided unadjusted 95% CI.

exact performs exact estimation, rather than the default normal approximation. Available for subcommands prop, rate and smr for cross-sectional analyses,
and smr for longitudinal analyses.

winsor(#) indicates the percent of the observations to be winsorized in each tail. More specifically, all values above (below) the percent specified are
replaced by the value immediately below (above) the percent specified. winsor(#) must be less than 50. The default is no winsorizing.

overdisp(string) adjusts for overdispersion of the performance measure Y. The two options are additive overdisp(add)
and multiplicative overdisp(mult).

ypercent multiplies y-axis values by 100 to present the values on a whole number scale. ypercent is available
for subcommands prop and smr

ratedenom(#) multiplies y-axis values by the specified amount (e.g. ratedenom(1000)) and divides x-axis values by that specified amount, indicating
that the values represent a rate. ratedenom(#) is available for subcommand rate.

ymin(#) presents observations and CI values only above the specified minimum value. For example, measures that cannot be negative (e.g. proportions and rates)
may be truncated at zero rather than showing lower CI values which may be negative.

ymax(#) presents observations and CI values only below the specified maximum value.

figure[(twoway_options)] allows the user to modify the default graph settings using all available twoway options.

Examples

1) funnelinst for a proportion

Set-up

. use cabg.dta, clear

Data are 30-day mortality rates following coronary artery bypass grafts in New York state, 1997–1999, for 33 hospitals. The numerator here is
the unadjusted number of deaths and the denominator is the volume of cases. We specify that the maximum value on the x-axis is 2000, the minimum value
on the y-axis is 0 and the maximum value on the y-axis is 8, and that values on the y-axis be multiplied by 100 to present percentages as whole numbers.

. funnelinst prop deaths cases, xrangehigh(2000) ymin(0) ymax(8) ypercent

same as above but compute exact statistics instead of normal approximation, and add Y and X titles.

. funnelinst prop deaths cases, xrangehigh(2000) ymin(0) ymax(8) ypercent exact figure(ytitle(% mortality) xtitle(Volume of cases))

same as above but add marker labels to identify outliers.

. funnelinst prop deaths cases, xrangehigh(2000) ymin(0) ymax(8) ypercent exact figure(ytitle(% mortality) xtitle(Volume of cases) mlabel( hospital))

Winsorize 10% of high and low observations and adjust for multiplicative overdispersion.

. funnelinst prop deaths cases, xrangehigh(2000) ymin(0) ymax(8) ypercent winsor(10) over(mult) figure(ytitle(% mortality) xtitle(Volume of cases))

2) funnelinst for a rate

Set-up

. use teenage.dta, clear

Data are the under-18 conception rate for English health authorities, 2000–2001. The numerator here is the unadjusted number of teenage pregnancies as a rate per 1000.
The denominator is the population size. We specify that the maximum value on the x-axis is 26000 and the minimum displayed value on the y-axis is 0. We indicate
that the rate should be computed and displayed as "per thousand".

. funnelinst rate No2001 pop2001, xrangehigh(26000) pval(0.025 0.001) ymin(0) ratedenom(1000) figure(ytitle(Rate per 1000) xtitle(Population ('000s)))

Same as above but adjust for additive overdispersion after winsorizing 10% of the observations.

. funnelinst rate No2001 pop2001, xrangehigh(26000) pval(0.025 0.001) ymin(0) ratedenom(1000) wins(10) over(add) figure(ytitle(Rate per 1000) xtitle(Population ('000s)))

Same as above but specifies "madjust" to account for multiple testing rather than use the specified pval(#).

. funnelinst rate No2001 pop2001, xrangehigh(26000) ymin(0) ratedenom(1000) wins(10) over(add) figure(ytitle(Rate per 1000) xtitle(Population ('000s))) madjust

3) funnelinst for an indirect standardized rate (SMR)

Set-up

. use cabg.dta, clear

Data are 30-day adjusted and unadjusted mortality rates following coronary artery bypass grafts in New York state, 1997–1999. Here we use all funnelinst
default settings.

. funnelinst smr deaths expected

Here we set a target range to be between 0 and 2, and modify the titles on the X and Y axes.

. funnelinst smr deaths expected, target(0 2) figure(ytitle(Observed number of deaths) xtitle(Expected number of deaths))

Same as above, but we specify that 'exact' statistics be used to construct the CIs

. funnelinst smr deaths expected, target(0 2) figure(ytitle(Observed number of deaths) xtitle(Expected number of deaths)) exact

4) funnelinst for a continuous measure

Set-up

. use mean.dta, clear

These are artificial data of a hypothetical performance indicator measured on a continuous scale. We start by using all funnelinst
default settings.

. funnelinst mean means se

We now adjust for additive overdispersion, winsoring 10% of the observations.

. funnelinst mean means se, wins(10) over(add)

Same as above but we set the target to range between 4.5 and 5.5.

. funnelinst mean means se, wins(10) over(add) target(4.5 6.5)

5) funnelinst for a change over time

Set-up

. use teenage.dta, clear

We revisit the data for the under-18 conception rate for English health authorities. However, we now compare the
change in conception rates from 1998 to 2001 as a ratio (2001/1998). We use the government's target of a 7.5% reduction
(ratio = 0.925) as the target. The resulting funnel plot is exactly as that presented in Figure 4 of Spiegelhalter (2005a).

. funnelinst smr O1 E1 O2 E2, xrangehigh(1200) ymin(0.6) ymax(1.6) target(0.925) ///

figure(xlabel(0(200)1200) xtitle("Average observed") ytitle("Ratio of Standardised rates (2001/1998)")) exact

Here we compute the ratio of the proportions between 1998 and 2001.

. funnelinst prop No1998 pop1998 No2001 pop2001 , xrangehigh(30000) figure(xtitle("Average observed") ///

ytitle("Ratio of proportions (2001/1998)")) chtype(ratio) ymax(2)

Here we compute the difference in the proportions between 1998 and 2001.

. funnelinst prop No1998 pop1998 No2001 pop2001 , xrangehigh(30000) figure(xtitle("Average observed") ///

ytitle("Difference in proportions (2001 vs 1998)")) chtype(diff) ymax(2)

Stored results

funnelinst stores the following in r():

Scalars

r(N) The number of units

r(target) The target (if there is only one)

r(target1) The lower specified target (if there are two specified)

r(target2) The higher specified target (if there are two specified)

r(winsor) The percent of observations winsorized

r(phi) The dispersion value

r(tau) The random-effects standard deviation (from additive overdispersion)

r(pbon) The Bonferroni adjusted p-value (if madjust is specified)

r(pfdr) The FDR adjusted p-value (if madjust is specified)

References

Jones, H. E., Ohlssen, D. I. and D. J. Spiegelhalter. 2008. Use of the false discovery rate when comparing multiple health care providers.
Journal of Clinical Epidemiology 61: 232-240.

Spiegelhalter, D. J. 2005a. Funnel plots for comparing institutional performance. Statistics in Medicine 24: 1185-1202.

Spiegelhalter, D. J. 2005b. Handling over-dispersion of performance indicators. BMJ Quality & Safety 14: 347-351.

Spiegelhalter, D. J., Sherlaw-Johnson, C., Bardsley, M., Blunt, I., Wood, C. and O. Grigg. 2012. Statistical methods for healthcare regulation: rating, screening and surveillance.
Journal of the Royal Statistical Society Series A: Statistics in Society 175: 1-47.

Storey, J. D. and R. Tibshirani. 2003. Statistical significance for genomewide studies. Proceedings of the National Academy of Sciences 100: 9440-9445.

Citation of funnelinst

funnelinst is not an official Stata command. It is a free contribution
to the research community, like a paper. Please cite it as such:

Linden A. & D. J. Spiegelhalter (2024). FUNNELINST: Stata module for generating a funnel plot to compare institutional performance.
https://ideas.repec.org/c/boc/bocode/s459283.html

Authors

Ariel Linden

Linden Consulting Group, LLC

alinden@lindenconsulting.org

David J Spiegelhalter

Statistical Laboratory

Centre for Mathematical Sciences

University of Cambridge

Acknowledgments

We thank Daniel Klein for providing the mata code to compute the inverse CDF of the binomial distribution and Poisson distribution for the "exact" option. We thank
Hayley Jones for responding to questions about implementing FDR. We also wish to thank John Moran for advocating that we write this package.

Also see

Online: meta funnelplot, funnelcompar if installed