# How to draw funnel plot using ggplot2 in R?

As title, I need to draw something like this:

Can ggplot, or other packages if ggplot is not capable, be used to draw something like this?

-
I've got a few ideas about how to do and implement this, but would appreciate having some data to play with. Any ideas on that? –  Chase Dec 7 '10 at 1:53
Yes, ggplot can easily draw a plot that is made up of points and lines ;) geom_smooth will get you 95% of the way - if you want more advice you'll need to provide more details. –  hadley Dec 7 '10 at 2:05
This is not a funnel plot. Instead, the lines evidently are constructed from estimates of standard errors based on the number of admissions. They seem intended to enclose a specified proportion of data, which would make them tolerance limits. They are likely of the form y = baseline + constant / Sqrt(# admissions * f(baseline)). You could modify the code in the existing responses to graph the lines, but you likely would need to supply your own formula to compute them: the examples I have seen plot confidence intervals for the fitted line itself. That's why they look so different. –  whuber Dec 7 '10 at 17:32
@whuber (+1) That's a very good point, indeed. I hope this might provide a good starting point anyway (even if my R code isn't that optimized). –  chl Dec 7 '10 at 20:34
Ggplot still provides stat_quantile() to put conditional quantiles on a scatterplot. You can then control the functional form of the quantile regression with the formula parameter. I'd suggest things like formula= y~ns(x,4) to get a smooth splined fit. –  Shea Parkes Mar 18 '13 at 14:40

Although there's room for improvement, here is a small attempt with simulated (heteroscedastic) data:

library(ggplot2)
set.seed(101)
x <- runif(100, min=1, max=10)
y <- rnorm(length(x), mean=5, sd=0.1*x)
df <- data.frame(x=x*70, y=y)
m <- lm(y ~ x, data=df)
fit95 <- predict(m, interval="conf", level=.95)
fit99 <- predict(m, interval="conf", level=.999)
df <- cbind.data.frame(df,
lwr95=fit95[,"lwr"],  upr95=fit95[,"upr"],
lwr99=fit99[,"lwr"],  upr99=fit99[,"upr"])

p <- ggplot(df, aes(x, y))
p + geom_point() +
geom_smooth(method="lm", colour="black", lwd=1.1, se=FALSE) +
geom_line(aes(y = upr95), color="black", linetype=2) +
geom_line(aes(y = lwr95), color="black", linetype=2) +
geom_line(aes(y = upr99), color="red", linetype=3) +
geom_line(aes(y = lwr99), color="red", linetype=3)  +
annotate("text", 100, 6.5, label="95% limit", colour="black",
size=3, hjust=0) +
annotate("text", 100, 6.4, label="99.9% limit", colour="red",
size=3, hjust=0) +
labs(x="No. admissions...", y="Percentage of patients...") +
theme_bw()


-

If you are looking for this (meta-analysis) type of funnel plot, then the following might be a starting point:

library(ggplot2)

set.seed(1)
p <- runif(100)
number <- sample(1:1000, 100, replace = TRUE)
p.se <- sqrt((p*(1-p)) / (number))
df <- data.frame(p, number, p.se)

## common effect (fixed effect model)
p.fem <- weighted.mean(p, 1/p.se^2)

## lower and upper limits for 95% and 99.9% CI, based on FEM estimator
number.seq <- seq(0.001, max(number), 0.1)
number.ll95 <- p.fem - 1.96 * sqrt((p.fem*(1-p.fem)) / (number.seq))
number.ul95 <- p.fem + 1.96 * sqrt((p.fem*(1-p.fem)) / (number.seq))
number.ll999 <- p.fem - 3.29 * sqrt((p.fem*(1-p.fem)) / (number.seq))
number.ul999 <- p.fem + 3.29 * sqrt((p.fem*(1-p.fem)) / (number.seq))
dfCI <- data.frame(number.ll95, number.ul95, number.ll999, number.ul999, number.seq, p.fem)

## draw plot
fp <- ggplot(aes(x = number, y = p), data = df) +
geom_point(shape = 1) +
geom_line(aes(x = number.seq, y = number.ll95), data = dfCI) +
geom_line(aes(x = number.seq, y = number.ul95), data = dfCI) +
geom_line(aes(x = number.seq, y = number.ll999), linetype = "dashed", data = dfCI) +
geom_line(aes(x = number.seq, y = number.ul999), linetype = "dashed", data = dfCI) +
geom_hline(aes(yintercept = p.fem), data = dfCI) +
scale_y_continuous(limits = c(0,1.1)) +
xlab("number") + ylab("p") + theme_bw()
fp


-
The presence of the linetype=2 argument inside the aes() brackets - plotting the 99% lines - gives rise to an error "continuous variable cannot be mapped to linetype" with current ggplot2 (0.9.3.1). Amending geom_line(aes(x = number.seq, y = number.ll999, linetype = 2), data = dfCI) to geom_line(aes(x = number.seq, y = number.ll999), linetype = 2, data = dfCI) works for me. Feel free to amend the original answer and lose this. –  user22153 Mar 18 '13 at 8:39
@anxiousmac Thanks, it's now fixed. –  Bernd Weiss Mar 19 '13 at 5:37

See also the cran package berryFunctions, which has a funnelPlot for proportions without using ggplot2, if anyone needs it in base graphics. http://cran.r-project.org/web/packages/berryFunctions/index.html

There is also the package extfunnel, which I haven't looked at.

-