Using models from the drc package (for dose-response curves) with ggplot

Using models from the drc package with ggplot

# 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)

hi
This ggplot approach is very helpful. It worked with my data very well until the last stage where I encountered errors regarding the Coefs. I have it defined as yours but somehow R does not recognize it.
Again, when I run the 'df2 %>% unnest(coefs) %>% spread(names, x)' code the upper limit of my data did not show. There was a warning with respect to that.
Please help. Thank you

time insecticide mortality  sex age

33 15 dieldrin 3 male 8
34 15 dieldrin 4 male 8
35 15 dieldrin 6 male 8
36 15 dieldrin 5 male 8
37 30 dieldrin 2 male 8
39 30 dieldrin 2 male 8
41 45 dieldrin 4 male 8
42 45 dieldrin 2 male 8
44 45 dieldrin 2 male 8
45 60 dieldrin 3 male 8
46 60 dieldrin 2 male 8
161 60 ddt 1 male 8
162 60 ddt 2 male 8
163 60 ddt 2 male 8
164 60 ddt 1 male 8
165 120 ddt 3 male 8
167 120 ddt 2 male 8
169 180 ddt 3 male 8
170 180 ddt 1 male 8
171 180 ddt 1 male 8
173 240 ddt 3 male 8
174 240 ddt 8 male 8
175 240 ddt 6 male 8
176 240 ddt 5 male 8

Thank you! Your code worked really well with my data! I've got one more query: How could I add the confidence intervals to the plot?

You could use my package tidydrc.
Alternatively, you can use broom which now I think supports drc models...

First, thank you so much for uploading your code. It helped me a lot. Second, I have the same question as @dandradem had. I tried your suggested packages, but I could not figure it out. Is there already a code or do you maybe have some suggestions? Thank you in advance !

Hi Emma-S-lab,
couple of options with increasing level of difficulty:

• you could try a Shiny app I have written which is running here
• use the function tidydrc_plot function from the tidydrc package with confint = TRUE, e.g.

tidydrc_model(Puromycin, conc, rate, model = MM.3(), state) %>% 
tidydrc_plot(color = ~state, confint = T)

• add the confidence intervals yourself using geom_ribbon from ggplot2, you can take a look here for an example

Hope this helps.