library(tidyverse)
theme_set(theme_minimal())

points <- tibble(
  age = c(38, 45, 52, 61, 80, 74), 
  prop = c(0.146, 0.241, 0.571, 0.745, 0.843, 0.738))

colors <- list(
  data = "#41414550",
  # data = "grey80",
  fit = "#414145")

ggplot(points) + 
  aes(x = age, y = prop) + 
  geom_point(size = 3.5, color = colors$data) +
  scale_x_continuous(
    name = "Age in months", 
    limits = c(0, 96), 
    # Because age is in months, I want breaks to land on multiples
    # of 12. The `Q` in `extended_breaks()` are "nice" numbers to use
    # for axis breaks.
    breaks = scales::extended_breaks(Q = c(24, 12))) + 
  scale_y_continuous(
    name = "Intelligibility",
    limits = c(0, NA),
    labels = scales::percent_format(accuracy = 1))

 sar_data_drc <- sar_raw_file %>% mutate(Mean_value = rowMeans(dplyr::select(sar_raw_file, -c("Concentration (nM)", "DEG_ID")), na.rm = TRUE))
        
drc_model <- drm(Mean_value~ `Concentration (nM)`, DEG_ID, fct=LL.4(names=c("Slope", "Lower", "Upper", "ED50")), data = sar_data_drc)

        results <- drc_model$coefficients %>%
          as.data.frame() %>%
          rownames_to_column() %>% setNames(c("coefficients", "value")) %>%
          mutate(coeff = str_split_fixed(coefficients, pattern = ":",n = 2)[,1]) %>%
          mutate(DEG_ID = str_split_fixed(coefficients, pattern = ":",n = 2)[,2]) %>%
          dplyr::select(c("value", "coeff", "DEG_ID")) %>%
          reshape2::dcast(DEG_ID~coeff) %>%
          rowwise() %>%
          mutate(IC50 = ((Upper -Lower )/(0.5-Lower ) -1)^(1/Slope) *ED50) %>%
          mutate(IC30 = ((Upper -Lower )/(0.7-Lower ) -1)^(1/Slope) *ED50) %>%
          mutate(IC10 = ((Upper -Lower )/(0.9-Lower ) -1)^(1/Slope) *ED50) %>% 
          mutate(AUC = computeAUC(c(10000, 10000/3, 10000/(3^2), 10000/(3^3), 
                                                10000/(3^4), 10000/(3^5), 10000/(3^6), 
                                                10000/(3^7), 10000/(3^8), 10000/(3^9)), 
                                              area.type = "Fitted", 
                                              conc_as_log = F, 
                                              Hill_fit = c(-Slope, Upper, ED50))) %>% 
          cbind(data.frame(drc::ED(drc_model, 50)[,1], drc::ED(drc_model, 90)[,1]) %>% as.data.frame() %>% setNames(c("EC50", "EC90"))) %>% 
          dplyr::rename("Emax" = "Lower") %>% 
          mutate_at(c("Slope", "Emax", "Upper", "ED50", "IC50", "IC30", "IC10", "AUC","EC50", "EC90"), ~round(., digits = 3))