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.
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 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 which denotes the slope around the point of inflection, which is the , 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.
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
Before curve fitting, it is important to get an idea of how the data look.
## 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 is denoted
e, and finally the slope,
b, around the .
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
The slope of the dose-response curve at 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 is:
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
#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:
- Correct regression model
- Variance Homogeneity
- Normally distributed measurement errors
- Mutually independent measurement error
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
One of the most popular ways of assessing goodness of fit in regression, be it linear or nonlinear, is to use . The higher the better, but the 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 will be larger than if the differences are small. It is also important to note that the most common interpretation of for linear regression does not hold true for nonlinear regression; in the nonlinear case, the value is not the amount of variability in Y explained by X. At its core, the 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 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 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 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.
## 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 is a parameter of the model. Several other models exist, particularly, the non symmetric ones where is not a “natural” parameter. For example the Weibull-1 model.
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 .The Weibull-1 descends slowly from the upper limit
d and approach the lower limit
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.
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 for the Weibull-1 and Weibull-2 model has to be derived from the fit by using the function
## ## 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 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 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 , and 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 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.
## Dose Herbicide DryMatter ## 1 0 Glyphosate 4.7 ## 2 0 Glyphosate 4.6
## 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
## 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
## ## 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
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 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 cannot be entertained.
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: , , and .
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.
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 for glyphosate and bentazon.
The most precise ED-level is at the ; 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 . In weed science 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 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
## ## 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 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 :
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 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 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 .
If this is the right way to compare absolute response level in Figure 11 we could use .
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")