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) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hi Emma-S-lab,
couple of options with increasing level of difficulty:
tidydrc_plot
function from the tidydrc package withconfint = TRUE
, e.g.geom_ribbon
fromggplot2
, you can take a look here for an exampleHope this helps.