### ABCBenchmark function ###

# A function so that the ABC (Achievable Benchmark of Care) can be calculated within a dplyr pipe workflow. 
# The function is based on calculating the Adjusted Performance Fraction (APF):

# Kiefe CI, Allison JJ, Williams OD, Person SD, Weaver MT, Weissman NW. 
# Improving quality improvement using achievable benchmarks for physician feedback: randomized controlled trial. 
# JAMA. 2001 Jun 13;285(22):2871-9.

### AND ###

# Paddock SM. 
# Statistical benchmarks for health care provider performance assessment: a comparison of standard approaches to a heirachical Bayesian histogram-based method. 
# HEALTH SERVICES RESEARCH. 2014 Jun; 49(3):1056-73.

ABCBenchmark <- function(.data, complete, total, grouping.variables, ...){
  require(dplyr)
  require(rlang)
  min10 <- function(x) min(x[x>=10], na.rm = T) # Get the minimum, given values are above 10  
  
  # Allow the function to be called in a dplyr workflow (magrittr pipes)
  if (dplyr::is_grouped_df(.data)) {
    return(dplyr::do(.data, ABCBenchmark(., ...))) 
  } 
  
  # Enquo the input variabels so they are read as columns
  complete = enquo(complete)
  total = enquo(total)

  if(hasArg('grouping.variables')){
  .data <- .data %>%
    left_join(.data %>%
    mutate(APF = (1+!!complete)/(2+!!total)) %>% #Calculate adjusted performance fraction
    arrange(desc(-APF)) %>%
    group_by_at(vars(!! grouping.variables)) %>%
    mutate(total.cumsum = cumsum(!!total),
           perc.contrib = 100*total.cumsum/sum(!!total, na.rm = T),
           min.contrib = min10(perc.contrib)) %>%
    filter(perc.contrib < 10 | perc.contrib == min.contrib) %>% #Use the top performing sites that account for at least 10 % of the patients
    mutate(ABCBenchmark = 100*sum(!!complete, na.rm = T)/sum(!!total, na.rm = T)) %>%
      select_at(vars(!! grouping.variables, ABCBenchmark)),
    by = grouping.variables)
  }
  
  else{
    .data$ABCBenchmark <- .data %>%
                  mutate(APF = (1+!!complete)/(2+!!total)) %>% #Calculate adjusted performance fraction
                  arrange(desc(-APF)) %>%
                  mutate(total.cumsum = cumsum(!!total),
                         perc.contrib = 100*total.cumsum/sum(!!total, na.rm = T),
                         min.contrib = min10(perc.contrib)) %>%
                  filter(perc.contrib < 10 | perc.contrib == min.contrib) %>% #Use the top performing sites that account for at least 10 % of the patients
                  mutate(ABCBenchmark = 100*sum(!!complete, na.rm = T)/sum(!!total, na.rm = T)) %>%
                  .$ABCBenchmark %>% unique()
  }
  
  .data
}

### EXAMPLE ###

# Example of how the function works

## Make the dataframe
healthdata <- data.frame(
ClinicianID = sample(1:100, size = 100, replace = F),
pop.complete = sample(1:10, size = 100, replace = T), #Mock numerator data
pop.total = sample(10:25, size = 100, replace = T), #Mock denominator data
site = sample(LETTERS[1:10], size = 100, replace = T),
region = sample(c("North", "West", "South", "East"), size = 100, replace = T))

# Using the function
healthdata_with_aspirational_target <- ABCBenchmark(healthdata, pop.complete, pop.total, grouping.variables = c('site', 'region'))
# Or using pipes
healthdata_with_aspirational_target <- healthdata %>% ABCBenchmark(pop.complete, pop.total, grouping.variables = c('site', 'region'))
# Without grouping variables
healthdata_with_aspirational_target <- healthdata %>% ABCBenchmark(pop.complete, pop.total)