Skip to content

Instantly share code, notes, and snippets.

@angelovangel
Last active Nov 23, 2017
Embed
What would you like to do?
Fit dose-response models (with drc) and plot, for many samples
# 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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment