Fit dose-response models (with drc) and plot, for many samples

do_drm.R

# these two functions perform the following:
# do_drm() --> performs drm (from the drc package) on a long dataframe, using the LL.4 log-logistic model for describing dose-response relationships
# do_drm_plot() --> plots the output of do_drm()
# to use the functions just paste this file in your R session and source it.
# try it out with
# do_drm(S.alba, Dose, DryMatter, Herbicide) %>% do_drm_plot(ed50 = TRUE, color = ~Herbicide) + scale_x_log10()

##=====================
# do_drm() usage
# do_drm(df, d, r, x, y, ...)

# the fourth (and subsequent) arguments can be used for grouping by sample, treatment etc. before modelling
# if the fourth argument is not supplied, it fits a drm to all the dose-response values

# df is the long dataframe
# d is the name of the dose column
# r is the name of the response column
# x, y, etc are grouping variables (columns containing grouping information)

# for example
# do_drm(S.alba, Dose, DryMatter, Herbicide)

# if you want a table with the coefficients, pipe the output of this function to
#  unnest(coefs) %>% spread(names, x)

# if you want a nice plot, pipe it to
# do_drm_plot()
##=====================
do_drm <- function(df, d, r, ...) {
# using eclipsis (...) for group_by(...), so that arbitrary number of grouping variables can be used
# like this dplyr works in the function, no need for quosure and !!
  
require(tidyverse)
require(drc)
require(modelr)
require(broom)
  
# rename columns to dose and response in order for the drm to work (formula interface problems)
d <- deparse(substitute(d))
r <- deparse(substitute(r))
  
df$dose <- df[[d]]
df$response <- df[[r]]
  
  drm.func <- function(x) {
    drm(response ~ dose, 
      fct = LL.4(names = c("Slope", "Lower Limit", "Upper Limit", "ED50")),
      data = x)
}

predict.fun <- function(x) {
  add_predictions(data.frame(dose = seq(0,max(df$dose))), x)
}

coefs.fun <- function(x) {coef(x) %>% tidy}

df %>% group_by(...) %>% nest() %>% 
  mutate(drmod = map(data, drm.func), 
         pred = map(drmod, predict.fun),
         coefs = map(drmod, coefs.fun))

}

# do_drm_plot
# this function plots the results from the
# do_drm() function, which is a df with list-columns
#=========================
# do_drm_plot() usage
# do_drm_plot(df, ed50 = FALSE , aes_arguments)

# for example
# do_drm_plot(df, ed50 = TRUE, color = ~Herbicide) + facet_wrap(...) and so on
# note the ~Herbicide notation here (check the manual for aes_) 
#=========================

do_drm_plot <- function(df, ed50 = FALSE, ...) {
  
  require(tidyverse)
  
  p <- df %>% unnest(data) %>% 
    ggplot() + 
    geom_point(aes_(~dose,  ~response, ...), alpha = 0.6) +  
    geom_line(aes_(~dose, ~pred, ...), data = df %>% unnest(pred)) +
    theme_bw()
  ifelse(ed50 == FALSE, 
         p, 
         p <- p + geom_vline(aes_(xintercept = ~`ED50:(Intercept)`, ...), 
                             linetype = 5,
                             data = df %>% unnest(coefs) %>% spread(names, x))
  ) 
  return(p)
  
}