## ----include = FALSE, echo=FALSE, fig.height= 5, fig.width=7------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.height = 5,
  fig.width = 7
)

## ----basicfunnel, warning=FALSE, error=FALSE, message=FALSE, fig.height= 5, fig.width=7----
library(ggplot2)
library(tidyr)

# Make up some data, as if it was from a regression model
# with observed and predicted (expected) events.
dt <-
  data.frame(
    observed = c(15, 40, 72, 28, 50, 66, 75),
    expected = c(13, 32, 75, 33, 54, 60, 72),
    unit = factor(c("A", "B", "c", "D", "E", "F", "G"))
  )

# Add a ratio (SR) of observed to expected, our indicator
dt$SR <- dt$observed / dt$expected

# Scatter plot in ggplot
a <-
  ggplot(dt, aes(x = expected, y = SR)) +
  geom_point()

a

# Now add a central line, as 1 is the average/expected value in this case.
a <-
  a +
  geom_hline(aes(yintercept = 1))

a

# Add a 95% Poisson limit, by using the density function to get the
# quantile value for each 'expected'.
lkup <- data.frame(id = seq(1, max(dt$expected), 1))
lkup$Upper <- (qpois(0.975, lambda = lkup$id) - 0.025) / lkup$id
lkup$lower <- (qpois(0.025, lambda = lkup$id) - 0.975) / lkup$id

lkup <- gather(lkup, key, value, -id)

a + geom_line(aes(x = id, y = value, col = key), data = lkup)


## ----install, eval=FALSE------------------------------------------------------
# devtools::install_github("https://github.com/nhs-r-community/FunnelPlotR")

## ----data, warning=FALSE, message=FALSE---------------------------------------
library(FunnelPlotR)
library(COUNT)
library(ggplot2)

data(medpar)
medpar$provnum <- factor(medpar$provnum)
medpar$los <- as.numeric(medpar$los)

mod <- glm(los ~ hmo + died + age80 + factor(type)
           , family = "poisson"
           , data = medpar)
summary(mod)

## ----prediction---------------------------------------------------------------
medpar$prds <- predict(mod, type = "response")

## ----funnel1, message=FALSE, fig.align='center', fig.retina=5, collapse=TRUE----
funnel_plot(
  medpar, numerator = los, denominator = prds, group = provnum
  , title = "Length of Stay Funnel plot for `medpar` data"
  , draw_unadjusted = TRUE, draw_adjusted = FALSE
  , label = "outlier", limit = 99
)

## ----ODcheck, message=FALSE---------------------------------------------------
sum(mod$weights * mod$residuals^2) / mod$df.residual

## ----funnel2, message=FALSE, fig.align='center', fig.retina=5, collapse=TRUE----
funnel_plot(
  medpar, numerator = los, denominator = prds, group = provnum
  , title = "Length of Stay Funnel plot for `medpar` data"
  , draw_unadjusted = FALSE, draw_adjusted = TRUE, data_type = "SR"
  , sr_method = "SHMI", label = "outlier", limit = 99
)