Last active
June 23, 2023 13:43
-
-
Save angelovangel/c69d78a27c4360df3057b40f4a705837 to your computer and use it in GitHub Desktop.
Using models from the drc package (for dose-response curves) with ggplot
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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 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 withconfint = 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
fromggplot2
, you can take a look here for an example
Hope this helps.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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 !