Hi Danilo,

There are standard models for dose response cytotoxicity curves, which I suggest you use. For comparing toxicicty between the different algae, I would suggest, you calculate IC10, IC50 and IC90 values, which can than be compared using one-way ANOVA.

To do this, you can use this equation:

ICF= ((100-F)/F)1/H * IC50

Where F equals the fraction of maximal response, e.g. for IC90, F=90 and H is the Hill's slope of the cytotoxicity curve.

Why re-invent the wheel?

Cheers,

Kirsten

Hi Danilo,

The drc package will allow you to compare many models besides logistic and allow estimates to be extracted and compared between curves. It has an advantage over other modeling methods particularly for binomial data because it doesn't correct background mortality at zero and low doses so you can test for differences in asymptotes of curves, in your case, between species. Once you have a model/s you're happy with you can compare slopes, IC10,50,90's, and asymptotes.

I've attached a copy of the manual here, which goes into details of modeling multiple curves. I think in your case, you may be using the wrong model. From memory LL.5 deals with hormetic curves where low concentrations can cause an increase in the response. You might consider using LL.3 or Wiebull three parameter models with the lower asymptote fixed at zero survival.

The most recent Ritz et al paper 2015 also has some nice code for fancy dose response plots using ggplot

Thanks Kirsten for the suggestion.

Hi Alan,

thanks for the info. I am comparing the models using the AICc number and using visual observation of the curves in the plot. You are right, i tried LL.3 function and fits better the curves.

I will have a look to Ritz et al paper to plot it using ggplot.

Thanks again

Danilo

G'day all,

Danilo, I hope that my opinions and advice are useful to you. I don't do very much dose response work, precisely, but nls and drc have proven very useful in other physiological applications that are characterised by slope-threshold functions.

Anyway, my advice to using drc is that you should become familiar with the whole concept of non-linear regression in general. Christian Ritz and Jens Streibig have written a fantastic book on this (see below), which explains the longhand way of doing many of the operations that you're asking about. Once you understand how the approach works, you'll have a better feeling for exactly what drc is doing, because it is essentially just automating many of the processes.

Ritz, C. and Streibig, J.C. (2008) Nonlinear Regression with R, Springer.

The simple answer to your question is that nls (and drc) operates in a different manner to linear regression. Instead of being an ANOVA-oriented approach, nls analysis operates by model-comparison, as you have found with your AIC. What this means is that a global model is fitted first, and then the data are grouped by factors (like your species) and subsequent models are tested against the global model to see if the goodness of fit has increased. In the drc, this is done by using the "comparm" command. You can fit unique models to each of your species, but in doing so you will have no statistical basis to understand whether the models really are different or not because you cannot compare different nls models by, say, ANCOVA. Also, by using "comparm" it is possible to interrogate which parameters in your fits are actually different, offering mechanistic insight into what these differences mean. For example, should your different species have different upper and lower asymptotes to your fits? I wouldn't have thought that the lower asymptote should be different, because dead is dead in all species, but as for the rest, I don't know. The scaling factor (slope) of each function might reasonably be different, and also the "e" factor, representing the IC50. As Alan said above, once you have these species-specific models resolved, you can interrogate their derivatives for other critical thresholds, such as IC10 and IC90.

The key point here is that this isn't reinventing the wheel. Calculating these thresholds alone is theoretically rooted in the assumption that dose response functions are non-linear, but the calculation of these thresholds, as I understand it, relies upon a normal and linear pattern of variance. Neither of those things is typical in untransformed nonlinear data. Therefore, it is actually more accurate to fit the drc function and interpolate from it. Furthermore, thresholds alone provide limited insight into biological responses to physiological stressors, because substantial differences can exist between species, with substantial ecological implications, well within very similar tolerance thresholds. Without fitting the drc function, and by just calculating threshold values by hand a great deal of insight is lost. I think you're going about this very much in the correct manner, but maybe you'll be more across what commands are available and what they are doing after a bit more reading. I'd suggest the following:

Ritz, C. et al. (2015) Dose-Response Analysis Using R. PLoS ONE 10, e0146021. doi:10.1371/journal.pone.0146021.

Ritz, C. and Streibig, J.C. (2012) Dose response curves and other nonlinear curves in Weed Science and Ecotoxicology with the add-on package drc in R.

Ritz, C. and Streibig, J.C., Package 'drc': Analysis of dose response curves in biology, environmental sciences, medicine, pharmacology, toxicology. R package version 2.3-96., 2013.

Good luck with it.

Sean

Hi Danilo,

 I enjoy working with the drc as it provides a very power method for analysing nonlinear data, especially those that are defined by increasing/ decreasing doses. I also agree with others that various thresholds are of great importance (such as the ED10, ED50 and ED90), but before you quantify these points you need to justify statistically the appropriate curve fit and in conjunction with other parameters, especially the slope! The drc is extremely user-friendly and a powerful tool for this type of analysis.

In your drc approach make sure you are using the right function, as the drc can provide a number of sigmoidal fits, ie. Weibull, logistic, loglogistic. The help package gives a great overview of the available functions. Then the next step is to describe your data in the most parsimonious way -  this means that you are required to do a model reduction. In your example you have shown a five parameter loglogistic function, which will have to be compared with a four, a three and perhaps also a two parameter model. When reducing your model, various parameters become 0. In the instance of your four parameter model, you have the upper limit (=d), lower limit (=c), slope (=b), and the inception point (=e, this parameter also is known as the ED50 or median point in the curve decay). The three parameter model assumes that your lower limit c is equal to 0, and the two parameter model you will only quantify the slope and the inception point (this is probably too simple).

 The first step is to compare the various parameter functions using the mselect() function (e.g. mselect(m1, list(LL.5(), LL.4(), LL.3(), LL.2()), linreg=TRUE, icfct=AIC). This will help you reduce the model to the appropriate parameter function, and in the example I have given you, based on the AIC, loglikelihood and lack of fit. Also note that the linreg =TRUE means that the curves are also compared against linear functions. The function will rank your curves based on the best fit. Caution here, as this is user defined, so you need to check that you are using a biologically appropriate function. You can also use mselect() to compare between different families if you need to.

 After you have done so run two models, the first with treatment as an effect, and the second that looks at the whole complete model without treatment, followed by an anova.drc() of the models. This will justify that significance of treatment (if it is a significant effect).

m1

Thanks Sean and Wolfgang for the well explained answer.

what you are suggesting is not too far from what i am doing but you just gave me some important tips that complete the analysis.

I am going to read to the paper you are suggesting.

Thanks again.

Danilo

Thanks for posting your challenges. I had similar issue but this is how I solved it.

When you install the dose response curve package in R, you might consider using install.packages("drc", dependencies =TRUE)

Setting dependencies to TRUE helped me to overcome such challenge and it worked out fine. Hope it helps.

Dear Sean Tomlinson,

An excellent answer, which did help me a lot on a similar, but different topic.

Thanks!

Kind regards, Harrie Verhoeven 

G'day Harrie,

Thanks for your response. I'm glad that it helped. Also, I tip my hat to Wolfgang, who works much more closely with drc than I do. I'm sure that his insights into using the package were probably helpful as well. Good luck with your research

Cheers

Sean

Hi all, this thread has helped me significantly with my own attempts to use drc for fitting my data, and inspired me to get a researchgate account, so thank you!

I have some related questions i've run into with my dose-response regression analysis: https://www.researchgate.net/post/Dose_response_curves_in_drc_package_in_R_Model_fitting_interpretation_and_comparison

And with the same data attempting to use drc to run a regression on time to event (mortality).

https://www.researchgate.net/post/dose_response_regression_analysis_for_time_to_mortality_data_in_drc_package_in_R

I would be grateful for any insight you may be able to provide on any aspect of these questions!