# the aim is to use the 'drc' package to fit models to data and then extract the data and use for plotting in ggplot
# the data could be growth curves, dose-response, Michaelis-Menten etc.
# here, the S.alba data from drc is used for dose-response

library(tidyverse)
library(broom)
library(drc)
library(modelr)
attach(S.alba) # the data used in this gist
library(egg)

df <- S.alba
# view the data in normal and log scale for dose
p1 <- df %>% ggplot() + geom_point(aes(Dose, DryMatter, color = Herbicide)) + theme_bw()
p2 <- df %>% ggplot() + geom_point(aes(log(Dose), DryMatter, color = Herbicide)) + theme_bw()
ggarrange(p1, p2)

# define drm function to use with map
drm.func <- function(x) {
  drm(DryMatter ~ 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,700)), x)
}

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

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

# plot raw data, model and ED50 line
df2 %>% unnest(data) %>%
  ggplot() +
     geom_point(aes(log(Dose), DryMatter, color = Herbicide)) +
     geom_line(aes(log(Dose), pred, color = Herbicide), data =
                                      df2 %>% unnest(pred)) +
     geom_vline(aes(xintercept = log(x), color = Herbicide),
                  linetype = 5,
                  data = df2 %>% unnest(coefs) %>% filter(names == "ED50:(Intercept)")) +
     theme_bw()

# summary of coefficients
df2 %>% unnest(coefs) %>% spread(names, x)