angelovangel / tidydrc

tidy modelling of dose-response relationships with the drc package

tidydrc

tidydrc

Description

Tidy modelling of dose-response relationships with the drc package in R.
This is a wrapper for drc by Christian Ritz, which is probably the best package for modelling dose-response.
tidydrc contains two functions which make it easier to generate and plot these models. Install it with:

devtools::install_github("angelovangel/tidydrc")

Examples

The tidydrc_model() function returns a dataframe with list-columns (the data, predictions and coefficients). It is thus easy to implement in tidy workflows. For example, to fit Michaelis-Menten kinetics models for treated and untreated samples in the Puromycin dataset (built-in) and get the Km values with std. error:

mm <- tidydrc_model(Puromycin, conc, rate, model = MM.3(), state)
names(mm$drmod) <- as.character(mm$state)
map(mm$drmod, ED, 50) %>% map_df(as_tibble, .id = "sample")

The list-column dataframe can be directly piped to tidydrc_plot()

tidydrc_model(S.alba, Dose, DryMatter, model = LL.4(), Herbicide) %>%
tidydrc_plot(ed50 = TRUE, color = ~Herbicide) + 
scale_x_log10()

A more involved example...

References

Ritz, C., Baty, F., Streibig, J. C., Gerhard, D. (2015) Dose-Response Analysis Using R PLOS ONE, 10(12), e0146021

tidydrc's Issues

ED50 extrapolation is not where the curve intersects y=50

Hey Angel - Thanks for the great work on this package!

This issue relates to the ED50 calculation being off from where the curve intersects at y=50. Here, I use normalized data (constrained within 0-100) and attempt to extract an ED50 closer to the ground truth.

library(tidydrc)
#Make a data frame of dose-response results (response signifies pre-normalized %Cell Death values)
#As a ground truth, we have dose of 59.375 equals response of 50.000
mydr <- data.frame("dose"=c(1900,950,475,237.5,118.75,59.375,29.68,14.84,7.42,3.71), 
                   "response"=c(98.38,96.267,93.300,91.1,53.56,50.000,44.2,24.1,15.5,9.01))

#Apply the LL.4 model to fit a log-logistic curve
mydr.LL4 <- tidydrc_model(mydr, dose, response, model=LL.4())

#Extract the EC50 value
mydr_EC50 <- as.data.frame(as.numeric(mydr.LL4[[4]][[1]][4,2]))

#Plot the curve
#ed50=TRUE plots a dashed black line according to the extrapolated EC50 
#ground truth lines are in red. y-axis=50, and x-axis=59.375

tidydrc_plot(mydr.LL4,ed50=TRUE) + 
  scale_x_log10(labels=log10) + 
  xlab("Log10[concentration]")+
  ylab("% Cell Death") +
  geom_hline(aes(yintercept=50), linetype=1, color="red") +
  geom_vline(aes(xintercept=59.375), linetype=1, color="red")+
  geom_point(aes(x=dose, y=response), size=2.4, color="blue") +
  geom_text(aes(x=10,y=100,label=paste0("original_EC50: ",round(mydr_EC50,3))))

We can see the dashed line is off, calling the ec50 as 67.29, when we know it is 59.375
Could this be due to extrapolating 50% of the curves maximal effect, and not considering the constriants of 0-100 imposed by the pre-normalized data?

#Look at the predicted values of the fit to see if we can get closer to the 'true' ec50 value
mydr_predvals <- mydr.LL4$pred[[1]]
plot(mydr_predvals$pred)

#From the predicted values, get the dose where the response is closest to 50
new_EC50 <- mydr_predvals[which.min(abs(50-mydr_predvals$pred)),1] #note which.min alone returns the *position*, so doing it like this returns the actual value

#Re-make the plot but add new_EC50 (green) and see if it gets any closer
tidydrc_plot(mydr.LL4,ed50=TRUE) + 
  scale_x_log10(labels=log10) + 
  xlab("Log10[concentration]")+
  ylab("% Cell Death") +
  geom_hline(aes(yintercept=50), linetype=1, color="red") +
  geom_vline(aes(xintercept=59.375), linetype=1, color="red")+
  geom_point(aes(x=dose, y=response), size=2.4, color="blue") +
  geom_vline(aes(xintercept=new_EC50), linetype=5, color="darkgreen") +
  geom_text(aes(x=10,y=100,label=paste0("original_EC50: ",round(mydr_EC50,3)))) +
  geom_text(aes(x=10.9,y=94, label=paste0("new_EC50: ", round(new_EC50,3))))

The new EC50 gets much closer.
Thoeretically, some datasets may have a response point closer to 50 than exists in the predicted values

#Merge the original response with the predicted response and take whichever value of the combined results is closest to 50:
merged <- data.frame("dose"=abind(mydr_predvals$dose,mydr$dose), "response"=abind(mydr_predvals$pred, mydr$response))
new2_EC50 <- merged[which.min(abs(50-merged$response)),1]

#Re-make the plot and add the new2_EC50 value
tidydrc_plot(mydr.LL4,ed50=TRUE) + 
  scale_x_log10(labels=log10) + 
  xlab("Log10[concentration]")+
  ylab("% Cell Death") +
  geom_hline(aes(yintercept=50), linetype=1, color="red") +
  geom_vline(aes(xintercept=59.375), linetype=1, color="red")+
  geom_point(aes(x=dose, y=response), size=2.4, color="blue") +
  geom_vline(aes(xintercept=new_EC50), linetype=5, color="darkgreen") +
  geom_vline(aes(xintercept=new2_EC50), linetype=5, color="cyan", size=1.2)+
  geom_text(aes(x=10,y=100,label=paste0("original_EC50: ",round(mydr_EC50,3)))) +
  geom_text(aes(x=10.9,y=94, label=paste0("new_EC50: ", round(new_EC50,3)))) +
  geom_text(aes(x=10.6,y=88, label=paste0("new2_EC50: ", round(new2_EC50,3))))

Now we've achieved the ground truth. I realize that the likelihood of a single point in the observed data being closer to 50 than the best predicted value is pretty low, but merging both datasets and looking for the closest dose where response=50 seems to fix this. Either way, the two values aren't super different.

Best,
Charlie

Warning: Recycling array of length 1 in array-vector arithmetic is deprecated

Hi, first of all thank you for this package!

I created the issue because when I run the command below, I get the following warning. I guess that could be easily solved with replacing as.vector() with c().

tidydrc_model(response = efficacy, dose = concentration, data=compounddata, model=LL.4(fixed=c(NA,0,1,NA), names = c("Slope", "Lower Limit", "Upper Limit", "ED50")), robust="mean")

Warning message:
There were 180 warnings in `dplyr::mutate()`.
The first warning was:
ℹ In argument: `pred = purrr::map(drmod, predict.fun)`.
ℹ In group 1: `robust = "mean"`.
Caused by warning in `(tquan * sqrt(varVal + sumObjRV)) * c(-1, 1)`:
! Recycling array of length 1 in array-vector arithmetic is deprecated.
  Use c() or as.vector() instead.