Dose Response








Dose-Response curves

Finney (Bioassay and the Practice of Statistical Inference, International Statistical Review / Revue Internationale de Statistique, Vol. 47, No. 1 (Apr., 1979), pp. 1-12) wrote: “As a branch of applied statistics, biological assay appears to have a somewhat specialized appeal. Although a few statisticians have worked in it intensively, to the majority it appears as a topic that can be neglected, either because of its difficulties or because it is trivial……. I am convinced that many features of bioassay are outstandingly good for concentrating the mind on important parts of biometric practice and statistical inference.”
In biology there are four general curves that describe many biological phenomena. Three of them, the exponential, the Michaelis-Menten (rectangular hyperbola in weed science) and the sigmoid curves are having either one or two asymptotes; that is, they have upper and/or lower limits when the independent variable approaches very small or very large values.
In most text books, non-linear regression is not very well explained although it is easy for practitioners to implement. The difference between linear and non-linear regression can often be boiled down to the very fact: linear regression can be solved analytically by least square analysis, whereas the nonlinear regression cannot.
Before fitting a nonlinear regression model, one has to do some guesstimates about the initial regression parameters before fitting. In some cases a nonlinear relationship can be converted to a linear one. The most common one is the exponential curve (Figure 1) that turns into a linear relationship when taking the logarithm of the y provided the lower limit is zero. For the rectangular hyperbola it is usually with no lower limit and then it can also be linearize by manipulating as in ezymology and biochemistry.The rectangular hyperbola is in other disciplines called the Michaelis Menten model.
This chapter will be devoted to the use of the add-on package `drc´ and will illustrate some of the facilities built into the package. For further information consult the latest publication on the package: Dose-Response Analysis Using R in PlosOne
List of model functions and corresponding names of some of the most important built-in models available in drc.

Figure 1. Common regressions models in the biological sciences.
In the
For the sigmoid curves, which can take several forms we have to use nonlinear regression. One of the most common curves is the symmetric log-logistic model

y = c+\frac{d-c}{1+exp(b(log(x)-log(ED50)))} .


y is the response, c denotes the lower limit of the response when the dose x approaches infinity or zero depending on the sign of b. b which denotes the slope around the point of inflection, which is the ED_{50}, i.e. the dose required to reduce the response half-way between the upper and lower limit. If ´b´is negative the upper limit is at infinite x and if bis positive the upper limit is close to zero and the lower limit is when x approaches infinity ; d is the upper limit when the dose x approaches 0.
Guessing the parameters of a sigmoid curve can, for the less experienced person, be tiresome and frustrating. This has now been overcome by the development of a package named drc, which was developed at the University of Copenhagen. The purpose of the drc package was to ease the analysis of dose-response curves from greenhouse and field experiments in Weed and Pesticide Science courses and to assess the selectivity of herbicides and efficacy of pesticides. On the basis of a quick and dirty survey on the net the drcpackage commonly is used in ecotoxicology, toxicology, pharmacology, medicine and chemistry and various other disciplines.
Gradually, drc has been extended with numerous nonlinear curves commonly used in the life sciences and elsewhere. The engine of drc is the drm(y~x, fct=...) function with self-starter that automatically calculates initial parameters. The fct=argument identifies the particular dose-response curve to be used. An overview of the various functions available can be accessed, after loading the drc package, by typing ?getMeanFunctions() or consult the PlosOne paper.

One Dose-Response Curve

Let us look at how to fit a dose-response curve using a training data set in the drc package.

library(drc)

Before curve fitting, it is important to get an idea of how the data look.

head(ryegrass,3)
##      rootl conc
## 1 7.580000    0
## 2 8.000000    0
## 3 8.328571    0

and create a plot with the x variable both on original scale and on log scale.

op <- par(mfrow = c(1, 2)) #make two plots arranged in two columns 
plot(rootl ~ conc, data = ryegrass, main="Original Dose Scale")
plot(rootl ~ log(conc+.1), data = ryegrass, main="Logarithmic Dose Scale")


Figure 2. There are two ways of plotting dose-response data,note with the logarithmic dose scale the untreated control, 0, is not defined. That is why we add 0.1 to the conc before taking the logarithm.
Below we fit a four-parameter log-logistic model with user-defined parameter names. The default names of the parameters (b, c, d, and e) included in the drm() function might not make sense to many weed scientists, but the names=c() argument can be used to facilitate sharing output with less seasoned drc users. The four parameter log-logistic curve has an upper limit, d, lower limit, c, the ED_{50} is denoted e, and finally the slope, b, around the ED_{50}.

ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, 
                   fct = LL.4(names = c("Slope", "Lower Limit", "Upper Limit", "ED50")))
summary(ryegrass.m1)
## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##                         Estimate Std. Error  t-value p-value
## Slope:(Intercept)        2.98222    0.46506  6.41251  0.0000
## Lower Limit:(Intercept)  0.48141    0.21219  2.26876  0.0345
## Upper Limit:(Intercept)  7.79296    0.18857 41.32722  0.0000
## ED50:(Intercept)         3.05795    0.18573 16.46440  0.0000
## 
## Residual standard error:
## 
##  0.5196256 (20 degrees of freedom)

The summary of the curve fitting shows the estimates of each of the four parameters and their standard errors. The Slope is positive and the upper limit represents the untreated control and the lower limit denotes the Length of Roots at infinite doses. The p-values tell us if the parameters are different from zero Figure 3 . In this instance all four parameters are significantly different from zero and as seen on the graph the log-logistic curve seems to fit well to data.

op <- par(mfrow = c(1, 2))
plot(ryegrass.m1, broken=TRUE, 
     xlab="Concentration of Ferulic Acid", ylab="Length of Roots")
plot(ryegrass.m1, broken=TRUE, 
     xlab="Concentration of Ferulic Acid", ylab="Length of Roots",type="all")


Figure 3. Plot of regression and averages of observations within each concentration (left), and all observation (right). The default in the drc is on logarithmic scale for x.
The slope of the dose-response curve at ED_{50} has the opposite sign as compared to the sign of the parameter b. This is the consequence of the parameterization used in drc for the log-logistic model, a choice that is in part rooted in what was commonly used in the past. The actual slope of the tangent of the curve at ED_{50} is:

\frac{-b}{(d-c)/(4*e)}


So apart from the sign there is a scaling factor that converts the parameter b into the slope, therefore, the use of the word “relative”. Note that the scaling factor may be viewed as a kind of normalization factor involving the range of the response. As a consequence estimated b values often lie in the range 0.5 to 20 in absolute value, regardless of the assay or experiment generating the data.
At the end of the ryegrass help file ?ryegrass you can see examples, which give you some ideas of what can be done. The ryegrass help provides examples of various sigmoid curves to show the differences among them (see later). Note that we use the argument broken=TRUE to break the x-axis, because on a logarithmic scale, which is the default in drc, there is no such thing as a zero dose.
There are various ways to check if a regression model gives a satisfactory picture of the variation in data. Before doing further analysis, the best way to judge it is to look at the graph to see if there are some serious discrepancies between the regression and mean values. Being satisfied with the results of the fit, one way of checking is to test for lack of fit against the most general model, an ANOVA. It is done automatically in drc with the function modelFit().

#Test for lack of fit
modelFit(ryegrass.m1)
## Lack-of-fit test
## 
##           ModelDf    RSS Df F value p value
## ANOVA          17 5.1799                   
## DRC model      20 5.4002  3  0.2411  0.8665

Obviously, the test for lack of fit is non-significant and therefore we can tentatively entertain the idea that the regression analysis was just as good to describe the variation in data as was an ANOVA, which is the most general model we have; it only looks at means not relationships. The test for lack of fit is not enough to ensure the fit is reasonable, because of few degrees of freedom, in this instance only 3 degrees of freedom. Another issue is that the more doses we use the less likely it is that the test for lack of fit would become non-significant; the mere fact being that the log-logistic curve, like any other sigmoid curve, is just an approximation of the real relationship.
The distribution of the residuals and the hopefully normal distribution of residuals are informative to make sure that the prerequisites for doing a regression in the first place are not violated. The assumptions we must consider include:

  1. Correct regression model
  2. Variance Homogeneity
  3. Normally distributed measurement errors
  4. Mutually independent measurement error \varepsilon

Number 1 can always be discussed and we will not go too much into detail. This is generally determined by researchers familiar with the phenomena being described. Number 2 is important if you want to infer on the variability of parameters. Normally, the parameter estimates do not change much whether there are homogeneous variance or not; it is the standard errors that change and that is why we want to come so close to homogeneity of variance as possible. This is important when comparing curves, or evaluating the certainty we have about various parameters (like the effective dose or upper limit). Number 3 can be graphically shown and is only in some instances critical. Number 4 is important because we assume the observations at various doses are independent of each other, a question that can be debated. When we do dose-response experiments we often use the same stock solution. It means that any systematic error can be carried all the way though the dilution process and be part of the dose. Below we illustrate items 2 and item 3.

#Graphical analysis of residuals
op <- par(mfrow = c(1, 2)) #put two graphs together
    plot(residuals(ryegrass.m1) ~ fitted(ryegrass.m1), main="Residuals vs Fitted")
abline(h=0)
qqnorm(residuals(ryegrass.m1))
 qqline(residuals(ryegrass.m1))


Figure 4. Analysis of residuals. To the left the distribution of residuals against the fitted value (Variance Homogeneity) . To the right illustration of if the residuals are normally distribution (Normally distributed measurement errors).
The residual plot indicates that the model was able to catch the variation in data. The distribution of residuals should look like a shotgun shot on a wall. Apparently, there are some systematic violations in the residual plot. At small fitted values the residuals are close to zero than at high fitted values where the variation is much more pronounced. This is called a funnel type pattern. We also see it in the graphs with the fit and all the observations (Figure 3). The normal Q-Q plot looks all right, no serious deviation from the assumption of normal distribution of residuals. Later we will look at this problem, which is a statistical problem not a biological one.

Note on use of R^2

One of the most popular ways of assessing goodness of fit in regression, be it linear or nonlinear, is to use R^2. The higher the better, but the R^2 does not tell you how well the curve fits to the data. In fact, if there are large differences between the maximum and the minimum observation, the R^2 will be larger than if the differences are small. It is also important to note that the most common interpretation of R^2 for linear regression does not hold true for nonlinear regression; in the nonlinear case, the R^2 value is not the amount of variability in Y explained by X. At its core, the R^2 is a comparison of 2 models](http://www.jstor.org/stable/2684259); in linear regression the two models being compared are the regression model and a null model indicating no relationship between X and Y. For nonlinear regression, the choice of a ‘null’ model is not as simple, and therefore the R^2 type measure will be different depending on the chosen reduced model. The appropriate ‘null’ model may differ (and perhaps should differ) depending on the nonlinear model being used. Therefore, the R^2 is not part of the output from drc. To give an example we fit a linear model to the ryegrass data.

op <- par(mfrow = c(1, 2)) #put two graphs together
plot(rootl ~ conc, data = ryegrass)
linear.m1 <- lm(rootl ~ conc, data = ryegrass)
summary(linear.m1)
## 
## Call:
## lm(formula = rootl ~ conc, data = ryegrass)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4399 -1.7055  0.9623  1.7532  2.3575 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.24176    0.52973  11.783 5.64e-11 ***
## conc        -0.25929    0.04326  -5.994 4.94e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.07 on 22 degrees of freedom
## Multiple R-squared:  0.6202, Adjusted R-squared:  0.603 
## F-statistic: 35.93 on 1 and 22 DF,  p-value: 4.939e-06
abline(linear.m1)
plot(linear.m1, which=1)


Figure 5 Linear regression y=-0.25x+6.24 with highly significant parameters and residual plot. The R^2 is not neglect-able, 0.6, but the model is obviously not correct.
The regression parameters are highly significantly different from zero so it looks good as does R^2 of 0.62 . This is not a small value, but obviously the goodness of fit is not good, let alone the systematic distribution of residuals. When trying to determine whether the nonlinear model is better or worse than the linear model, Akaike Information Criterion (AIC) can be used. Smaller values of AIC indicate better fit to the data.

AIC(ryegrass.m1, linear.m1)
##             df       AIC
## ryegrass.m1  5  42.31029
## linear.m1    3 106.95109

The AIC value for the nonlinear model ryegrass.m1 is much lower than the linear fit linear.m1 that definitely violated item 1. This confirms that the nonlinear model is a better fit to the data compared with the linear model.

Which Dose-Response Model?

There exists numerous sigmoid curves of which the log-logistic is the most common one. It is symmetric and the ED_{50} is a parameter of the model. Several other models exist, particularly, the non symmetric ones where ED_{50} is not a “natural” parameter. For example the Weibull-1 model.

y = c+(d-c)exp(-exp(b(log(x)-log(e)))) .


The parameter b is the steepness of the curve and the parameter e is the dose where the inflection point of the dose- response curve is located, which is not the ED_{50}.The Weibull-1 descends slowly from the upper limit d and approach the lower limit c rapidly.
Another model, Weibull-2 is somewhat different from Weibull-1, because it has a different asymmetry, with a rapid descent from the upper limit d and a slow approach toward the lower limit Figure 6.

y = c+(d-c)(1-exp(-exp(b(log(x)-log(e))))) .


Ritz 2009 has in detailed described those asymmetric dose-response curves and also illustrated it as shown in Figure 6

library(  drc)
## Comparing log-logistic and Weibull models
## (Figure 2 in Ritz (2009))
ryegrass.m0 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())
ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, fct = W1.4())
ryegrass.m2 <- drm(rootl ~ conc, data = ryegrass, fct = W2.4())

plot(ryegrass.m0, broken=TRUE, xlab="Dose (mM)", ylab="Root length (cm)", lwd=2, 
     cex=1.2, cex.axis=1.2, cex.lab=1.2)
plot(ryegrass.m1, add=TRUE, broken=TRUE, lty=2, lwd=2)
plot(ryegrass.m2, add=TRUE, broken=TRUE, lty=3, lwd=2)

arrows(3, 7.5, 1.4, 7.5, 0.15, lwd=2)
text(3,7.5, "Weibull-2", pos=4, cex=1.2)

arrows(2.5, 0.9, 5.7, 0.9, 0.15, lwd=2)
text(3,0.9, "Weibull-1", pos=2, cex=1.2)


Figure 6. Comparison among the two asymmetric Weibull models and the log-logistic model. The differences in this instance are at either the upper limit or the lower limit of the curves.
The ED_{50} for the Weibull-1 and Weibull-2 model has to be derived from the fit by using the function ED().

 ED(ryegrass.m1,50,interval="delta")
## 
## Estimated effective doses
## 
##        Estimate Std. Error   Lower   Upper
## e:1:50  3.08896    0.17331 2.72744 3.45048

A relevant questions would be how do the ED_{x} change when fitting the different sigmoid curves?
To illustrate this we use the ggplot2 package and a function to combine graphs. First we have to calculate the ED_{x} and combine the results in a data.frame as done below.

edLL<-data.frame(ED(ryegrass.m0,c(10,50,90),interval="delta", display=FALSE),ll="Log-logistic")            
edW1<-data.frame(ED(ryegrass.m1,c(10,50,90),interval="delta", display=FALSE),ll="Weibull 1")
edW2<-data.frame(ED(ryegrass.m2,c(10,50,90),interval="delta", display=FALSE),ll="Weibull 2")
CombED<-rbind(edLL,edW1,edW2)

then we activate the function multiplot() in the file multiplot.R (if you wish to see the code you just paste the URL below in your browser). it is a function that combines various independent ggplots into one plot.

source("http://bioassay.dk/up to 18-6-15/Rstats4ag.Files/multiplot.R")
p1 <- ggplot(data=CombED[c(1,4,7),], aes(x=ll, y=Estimate))+
 geom_bar(stat="identity", fill="lightgreen", colour="black")+
  geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0.1) +
  ylab("ED10")+
  xlab("")
p2 <- ggplot(data=CombED[c(2,5,8),], aes(x=ll, y=Estimate))+
  geom_bar(stat="identity", fill="lightgreen", colour="black")+
  geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0.1) +
  ylab("ED50")+
  xlab("")

p3 <- ggplot(data=CombED[c(3,6,9),], aes(x=ll, y=Estimate))+
geom_bar(stat="identity", fill="lightgreen", colour="black")+
  geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0.1) +
  ylab("ED90")+
  xlab("Sigmoid four parameter models")

multiplot(p1,p2,p3)
## Loading required package: grid


Figure 7. Comparison of ED_{10}, ED_{50} and ED_{90} for various sigmoid dose-response curves with associated 95% confidence intervals.
Obviously, the ED-levels in Figure 7 in this instance did not change much among the three sigmoid curves. Sometimes they do and residual plots show that one of the asymmetrical curves is more appropriate than the common symmetric log-logistic, although it can not be substantiated by an ordinary test for lack of fit.drcalso has a function comped()where you can test the differences between the ED_{x} from individually fit as shown below

comped(CombED[c(1,4),1],CombED[c(1,4),2],log=F,operator = "-")
## 
## Estimated difference of effective doses
## 
##       Estimate Std. Error     Lower Upper
## [1,]  0.057726   0.314927 -0.559519 0.675
comped(CombED[c(1,7),1],CombED[c(1,7),2],log=F,operator = "-")
## 
## Estimated difference of effective doses
## 
##      Estimate Std. Error    Lower  Upper
## [1,] -0.16457    0.23335 -0.62193 0.2928
comped(CombED[c(3,6),1],CombED[c(3,6),2],log=F,operator = "-")
## 
## Estimated difference of effective doses
## 
##      Estimate Std. Error    Lower  Upper
## [1,]  1.28762    0.99032 -0.65337 3.2286
comped(CombED[c(6,9),1],CombED[c(6,9),2],log=F,operator = "-")
## 
## Estimated difference of effective doses
## 
##      Estimate Std. Error   Lower  Upper
## [1,]  -2.7048     1.4939 -5.6327 0.2231
comped(CombED[c(3,9),1],CombED[c(3,9),2],log=F,operator = "-")
## 
## Estimated difference of effective doses
## 
##      Estimate Std. Error   Lower Upper
## [1,]  -1.4172     1.6368 -4.6253 1.791
comped(CombED[c(6,9),1],CombED[c(6,9),2],log=F,operator = "-")
## 
## Estimated difference of effective doses
## 
##      Estimate Std. Error   Lower  Upper
## [1,]  -2.7048     1.4939 -5.6327 0.2231

None of the differences are significantly different form zero, lower confidence limits are negative in all instance and upper confidence limits of the differences are positive.

More Dose-Response Curves

When it comes to assessment of herbicide selectivity, one dose-response curve does not tell you anything. Selectivity of a herbicide is always dose-dependent and a bioassay with a pre-fixed harvest day only gives a snapshot of the effect of a herbicide. When comparing dose-response curves of say two herbicides on the same plant species or one herbicide on two plant species, it is imperative that the curves you compare were run at the same point in time or at least close to and at the same stage of plant development. The training data set S.alba in the drc package consists of two dose-response curves, one for bentazon and one for glyphosate.

head(S.alba,n=2)
##   Dose  Herbicide DryMatter
## 1    0 Glyphosate       4.7
## 2    0 Glyphosate       4.6
tail(S.alba,n=2)
##    Dose Herbicide DryMatter
## 67  640 Bentazone       0.6
## 68  640 Bentazone       0.8
## Fitting a log-logistic model with
S.alba.m1 <- drm(DryMatter ~ Dose, Herbicide, data=S.alba, fct = LL.4())

Note that the third argument in the drm() is the classification variable Herbicide.

modelFit(S.alba.m1)
## Lack-of-fit test
## 
##           ModelDf    RSS Df F value p value
## ANOVA          53 8.0800                   
## DRC model      60 8.3479  7  0.2511  0.9696
summary(S.alba.m1)
## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##               Estimate Std. Error   t-value p-value
## b:Glyphosate  2.715409   0.748279  3.628873   6e-04
## b:Bentazone   5.134810   1.130949  4.540266   0e+00
## c:Glyphosate  0.891238   0.194703  4.577429   0e+00
## c:Bentazone   0.681845   0.095111  7.168925   0e+00
## d:Glyphosate  3.875759   0.107463 36.066087   0e+00
## d:Bentazone   3.805791   0.110341 34.491101   0e+00
## e:Glyphosate 62.087606   6.611444  9.390929   0e+00
## e:Bentazone  29.268444   2.237090 13.083268   0e+00
## 
## Residual standard error:
## 
##  0.3730047 (60 degrees of freedom)

The test for lack of fit shows non-significance, and furthermore, the summary shows that the upper limits ds are very close to each other and the same applies to the lower limit cs. Consequently, we could entertain the idea that they are similar for both curves. We use the argument in the drm() function pmodels=data.frame(..)to make the upper limit and the lower limit similar for both curves. The argument follows the alphabetical order of the names of the parameters, b, c, d, e. Finally the e parameters can vary freely as well.

##  common lower and upper limits
S.alba.m2 <- drm(DryMatter ~ Dose, Herbicide, data=S.alba, fct = LL.4(),pmodels=data.frame(Herbicide, 1, 1, Herbicide)) 
summary(S.alba.m2)
## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##                Estimate Std. Error   t-value p-value
## b:Bentazone    5.046141   1.040135  4.851430       0
## b:Glyphosate   2.390218   0.495959  4.819387       0
## c:(Intercept)  0.716559   0.089245  8.029117       0
## d:(Intercept)  3.854861   0.076255 50.551925       0
## e:Bentazone   28.632355   2.038098 14.048566       0
## e:Glyphosate  66.890545   5.968819 11.206663       0
## 
## Residual standard error:
## 
##  0.3705151 (62 degrees of freedom)
plot(S.alba.m1, broken=TRUE)
plot(S.alba.m2, col="red", add=TRUE)


Figure 8. The two fits, one without assuming common upper and lower limits for the two curves and one where the curves have individual upper and lower limits.

anova(S.alba.m1, S.alba.m2) # Test for lack of fit not signifiancant
## 
## 1st model
##  fct:     LL.4()
##  pmodels: Herbicide, 1, 1, Herbicide
## 2nd model
##  fct:     LL.4()
##  pmodels: Herbicide (for all parameters)
## ANOVA table
## 
##           ModelDf    RSS Df F value p value
## 2nd model      62 8.5114                   
## 1st model      60 8.3479  2  0.5876  0.5588

A test for lack of fit was non-significant and, therefore, we entertain the idea that the curves have similar upper an lower limits. Now we have the ideal framework for comparing the relative potency between the two herbicides with similar upper and lower limit,i.e., we have the same reference for both curves.
Comparing the fit without freely estimated upper and lower limit and with common upper and lower limit did not change the ED_{50} much, but did reduce the standard errors somewhat.
The next step could be to assume that the two curves also have the same slope, i.e. the two curves are similar and only differ in their relative displacement along the x-axis.

##  common lower and upper limits and slope
S.alba.m3 <- drm(DryMatter~Dose, Herbicide, data=S.alba, fct = LL.4(),
                 pmodels=data.frame(1, 1, 1, Herbicide)) 
anova(S.alba.m2, S.alba.m3)
## 
## 1st model
##  fct:     LL.4()
##  pmodels: 1, 1, 1, Herbicide
## 2nd model
##  fct:     LL.4()
##  pmodels: Herbicide, 1, 1, Herbicide
## ANOVA table
## 
##           ModelDf    RSS Df F value p value
## 2nd model      63 9.4946                   
## 1st model      62 8.5114  1  7.1614  0.0095

Apparently, now the test for lack of fit is significant, so the assumption of identical curves (similar curves) except for the ED_{50} cannot be entertained.

plot(S.alba.m3, broken=TRUE)


Figure 9. The fit of the two curves are similar except for their relative displacement on the x-axis. However, here the test for lack of fit is highly significant so we have to resort to accepting only the upper and lower limits are common for the two curves.
Comparing Figures 8 and 9 shows that the fits are forced so they do not capture the responses properly at the lower limits.
As seen below the relative potency with common upper and lower limit but different slopes yield different relative potencies: ED_{10}, ED_{50}, and ED_{90}.

EDcomp(S.alba.m2, c(10, 10), interval="delta", reverse=TRUE)
## 
## Estimated ratios of effect doses
## 
##                            Estimate   Lower   Upper
## Glyphosate/Bentazone:10/10  1.44007 0.81836 2.06177
EDcomp(S.alba.m2, c(50, 50), interval="delta", reverse=TRUE)
## 
## Estimated ratios of effect doses
## 
##                            Estimate  Lower  Upper
## Glyphosate/Bentazone:50/50   2.3362 1.8571 2.8153
EDcomp(S.alba.m2, c(90, 90), interval="delta", reverse=TRUE)
## 
## Estimated ratios of effect doses
## 
##                            Estimate  Lower  Upper
## Glyphosate/Bentazone:90/90   3.7899 2.0793 5.5006

Above that the relative potency varies depending of the ED-level. This is of course due to the fact that the two curves have different slopes, b. From a biological point of view the herbicides have different mode of action and experience shows that PS II herbicides, such as bentazon, would have a steeper slope than does the EPSP inhibitor, glyphosate.
In drc there is a function relpot() that can illustrate the relationship between changes of the relative potency and the predicted responses.

relpot(S.alba.m2, interval = "delta")


Figure 10.The relative potency varies depending of response level, with the 95% confidence intervals. The Horizontal line shows the relative potency at the ED_{50} for glyphosate and bentazon.
The most precise ED-level is at the ED_{50}; towards the extremes close to the upper and lower limit the variation becomes huge. This is a general trend and explains why we often wish to quantify toxicity of a compound by using ED_{50}. In weed science ED_{90} is more appropriate as cont

When Upper and Lower Limits are not similar

In real life, e.g., when testing herbicide tolerant/resistant biotypes of weeds and different weed species, we cannot expect that the upper and lower limits are similar among curves and that the curves have the same slopes. On the contrary, if there is a penalty for acquiring resistance (fitness cost) then the untreated control for each biotype would be expected to have different biomass. Depending on the time from spraying and ambient growth condition at harvest; the lower limit may also differ.
Whatever the purpose of assessing potencies among plants the problem of different upper and lower limits is often “solved” in the literature by scaling the raw data to percentage of untreated control. This “solution” often results in smaller but incorrect standard error of ED_{50} being inflated by the scaling.
We do not want to go for low standard errors, but go for the correct standard errors. Therefore, we recommend to do all the curve-fitting on original raw data!
Below we have two accessions of Echonichlora oryzoides where both the upper limit and the lower limits are different.

op <- par(mfrow = c(1, 1))
TwoCurves<-read.csv("http://rstats4ag.org/data/twoCurves.csv")
head(TwoCurves,3)
##   xn Dose Accesion Dry.weight.per.pot  X
## 1 33    0  SAM-E-8              1.968 NA
## 2 34    0  SAM-E-8              2.030 NA
## 3 35    0  SAM-E-8              2.129 NA
Bispyribac.sodium.m1 <- drm(Dry.weight.per.pot ~ Dose, Accesion,
                            data = TwoCurves, fct=LL.4())
modelFit(Bispyribac.sodium.m1)
## Lack-of-fit test
## 
##           ModelDf    RSS Df F value p value
## ANOVA          48 1.7678                   
## DRC model      56 2.2799  8  1.7384  0.1137
summary(Bispyribac.sodium.m1)
## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##             Estimate Std. Error   t-value p-value
## b:SAM-E-8   2.047171   0.327585  6.249291  0.0000
## b:TEK-E-10  0.876716   0.315014  2.783102  0.0073
## c:SAM-E-8   0.245913   0.073938  3.325933  0.0016
## c:TEK-E-10 -0.053077   0.091626 -0.579276  0.5647
## d:SAM-E-8   2.242014   0.081575 27.484144  0.0000
## d:TEK-E-10  1.827553   0.101102 18.076407  0.0000
## e:SAM-E-8  21.702693   2.276642  9.532763  0.0000
## e:TEK-E-10  2.666463   1.026905  2.596602  0.0120
## 
## Residual standard error:
## 
##  0.2017751 (56 degrees of freedom)
plot(Bispyribac.sodium.m1, broken=TRUE, bp=.1, ylab="Dry Matter(g/pot)",
     xlab="Bispyribac sodium (g/ha)")

# drawing lines on the plot
ED.SAM.E.8 <- 0.24591 + (2.24201 - 0.24591) * 0.5
ED.TEK.E.10 <- -0.05308 + (1.82755 +- 0.05308) * 0.5
arrows(0.1, ED.SAM.E.8, 21.7, ED.SAM.E.8, code=0)
arrows(21.7, ED.SAM.E.8, 21.7, 0, code=0)
arrows(0.1, ED.TEK.E.10, 2.67, ED.TEK.E.10, code=0, lty=2)
arrows(2.67, ED.TEK.E.10, 2.67, 0, code=0, lty=2)


Figure 11. Fit of regression for two accessions with different upper and lower limit.
The fit looks reasonable and the test for lack of fit is non-significant and the ED_{50} s are shown. Everything seems to be all right, but we notice that distribution of observations for TEK-E-10 could be better to describe the mid range of the responses. The summary above reveals that the lower limit, c, for TEK-E-10 is less than zero, not much and not different from zero though, but still illogical. We re-fit the model using the lowerl= argument which will allow us to restrict model parameters to biologically relevant estimates. In this case, we will restrict the lower limit to values of zero or greater. All other parameters will still be allowed to vary without restriction.

Bispyribac.sodium.m2 <- drm(Dry.weight.per.pot ~ Dose, Accesion,
                            data = TwoCurves, fct=LL.4(),
                            lowerl=c(NA,0,NA,NA))
summary(Bispyribac.sodium.m2)
## 
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
## 
## Parameter estimates:
## 
##             Estimate Std. Error   t-value p-value
## b:SAM-E-8   2.042454   0.329052  6.207086  0.0000
## b:TEK-E-10  1.015226   0.313919  3.234041  0.0020
## c:SAM-E-8   0.243411   0.074596  3.263056  0.0019
## c:TEK-E-10  0.000000   0.071617  0.000000  1.0000
## d:SAM-E-8   2.242077   0.082057 27.323391  0.0000
## d:TEK-E-10  1.826855   0.101703 17.962711  0.0000
## e:SAM-E-8  21.740911   2.295727  9.470165  0.0000
## e:TEK-E-10  2.730186   0.961385  2.839846  0.0063
## 
## Residual standard error:
## 
##  0.2028892 (56 degrees of freedom)

We see that the standard error dropped a bit and now we use the fit and calculate the selectivity factor at ED_{50}:

EDcomp(Bispyribac.sodium.m2, c(50, 50),interval="delta")
## 
## Estimated ratios of effect doses
## 
##                        Estimate   Lower   Upper
## SAM-E-8/TEK-E-10:50/50   7.9632  2.0988 13.8275

Apparently, the resistance factor is about 8 and significantly different from 1. Note this is one of the few times we do not test the null hypothesis but the hypothesis, that the potency is the same for the two accessions; that is, the relative potency is 1.0. Another important issue is that the dry matter for the ED_{50} differs because of the different upper and lower limits. This will not be solved by scaling the observation to the individual accessions’ untreated control.
When using absolute dry matter of say 1.24, the response at ED_{50} for SAM-E-8 and TEK-E-10 is seen below.

ED(Bispyribac.sodium.m2, 1.24, type="absolute", interval="delta")
## 
## Estimated effective doses
## 
##                  Estimate Std. Error     Lower     Upper
## e:SAM-E-8:1.24  21.799437   2.302921 17.186132 26.412741
## e:TEK-E-10:1.24  1.306694   0.698047 -0.091663  2.705051

Looking at the distribution of observation for TEK-E-10, there are no observations to determine the shape of the curve above the ED_{50}.
If this is the right way to compare absolute response level in Figure 11 we could use ED_{1.24} .

plot(Bispyribac.sodium.m2,broken=TRUE,bp=.1, ylab="g Dry Matter/pot",
     xlab="Dose (g Bispyribac sodium/ha)")
arrows(  0.1, 1.24, 21.79, 1.24, code=0, lty=1, col="red")
arrows(21.79, 1.24, 21.79,    0, code=0, lty=1, col="red")
arrows(  1.31, 1.24,   1.31,    0, code=0, lty=2, col="red")