Package 'drc'

Description

Data from an experiment where the chemicals acifluorfen and diquat tested on Lemna minor. The dataset has 7 mixtures used in 8 dilutions with three replicates and 12 common controls, in total 180 observations.

Usage

data(acidiq)

Format

A data frame with 180 observations on the following 3 variables.

dose

a numeric vector of dose values

pct

a numeric vector denoting the grouping according to the mixtures percentages

rgr

a numeric vector of response values (relative growth rates)

Details

The dataset is analysed in Soerensen et al (2007). Hewlett's symmetric model seems appropriate for this dataset.

Source

The dataset is kindly provided by Nina Cedergreen, Department of Agricultural Sciences, Royal Veterinary and Agricultural University, Denmark.

References

Soerensen, H. and Cedergreen, N. and Skovgaard, I. M. and Streibig, J. C. (2007) An isobole-based statistical model and test for synergism/antagonism in binary mixture toxicity experiments, Environmental and Ecological Statistics, 14, 383–397.

Examples

## Fitting the model with freely varying ED50 values
## Ooops: Box-Cox transformation is needed
acidiq.free <- drm(rgr ~ dose, pct, data = acidiq, fct = LL.4(),
pmodels = list(~factor(pct), ~1, ~1, ~factor(pct) - 1))

## Lack-of-fit test
modelFit(acidiq.free)
summary(acidiq.free)

## Plotting isobole structure
isobole(acidiq.free, xlim = c(0, 400), ylim = c(0, 450))

## Fitting the concentration addition model
acidiq.ca <- mixture(acidiq.free, model = "CA")

## Comparing to model with freely varying e parameter
anova(acidiq.ca, acidiq.free)  # rejected

## Plotting isobole based on concentration addition -- poor fit
isobole(acidiq.free, acidiq.ca, xlim = c(0, 420), ylim = c(0, 450))  # poor fit

## Fitting the Hewlett model
acidiq.hew <- mixture(acidiq.free, model = "Hewlett")

## Comparing to model with freely varying e parameter
anova(acidiq.free, acidiq.hew)  # accepted
summary(acidiq.hew)

## Plotting isobole based on the Hewlett model
isobole(acidiq.free, acidiq.hew, xlim = c(0, 400), ylim = c(0, 450))  # good fit

Description

Dataset from an experiment exploring the effect of increasing concentrations of a herbicide on the volume of the treated algae.

Usage

data(algae)

Format

A data frame with 14 observations on the following 2 variables.

conc

a numeric vector of concentrations.

vol

a numeric vector of response values, that is relative change in volume.

Details

This datasets requires a cubic root transformation in order to stabilise the variance.

Source

Meister, R. and van den Brink, P. (2000) The Analysis of Laboratory Toxicity Experiments, Chapter 4 in Statistics in Ecotoxicology, Editor: T. Sparks, New York: John Wiley \& Sons, (pp. 114–116).

Examples

algae.m1 <- drm(vol~conc, data=algae, fct=LL.3())
summary(algae.m1)

algae.m2 <- boxcox(algae.m1)
summary(algae.m2)

Description

'anova' produces an analysis of variance table for one or two non-linear model fits.

Usage

## S3 method for class 'drc'
anova(object, ..., details = TRUE, test = NULL)

Arguments

Details

Specifying only a single object gives a test for lack-of-fit, comparing the non-linear regression model to a more general one-way or two-way ANOVA model.

If two objects are specified a test for reduction from the larger to the smaller model is given. (This only makes statistical sense if the models are nested, that is: one model is a submodel of the other model.)

Value

An object of class 'anova'.

Author(s)

Christian Ritz

References

Bates, D. M. and Watts, D. G. (1988) Nonlinear Regression Analysis and Its Applications, New York: Wiley \& Sons (pp. 103–104)

See Also

For comparison of nested or non-nested model the function mselectcan also be used.

The function anova.lm for linear models.

Examples

## Comparing a Gompertz three- and four-parameter models using an F test
ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, fct = W1.4())
ryegrass.m2 <- drm(rootl ~ conc, data = ryegrass, fct = W1.3())
anova(ryegrass.m2, ryegrass.m1)  # reduction to 'W1.3' not possible (highly significant)   

anova(ryegrass.m2, ryegrass.m1, details = FALSE)  # without details

Description

Providing the mean function and the corresponding self starter function for the asymptotic regression model.

Usage

AR.2(fixed = c(NA, NA), names = c("d", "e"), ...)

  AR.3(fixed = c(NA, NA, NA), names = c("c", "d", "e"), ...)

Arguments

Details

The asymptotic regression model is a three-parameter model with mean function:

$f(x) = c + (d-c)(1-\exp(-x/e))$

The parameter $c$ is the lower limit (at $x=0$), the parameter $d$ is the upper limit and the parameter $e>0$ is determining the steepness of the increase as $x$.

Value

A list of class drcMean, containing the mean function, the self starter function, the parameter names and other components such as derivatives and a function for calculating ED values.

Note

The functions are for use with the function drm.

Author(s)

Christian Ritz

See Also

A very similar, but monotonously decreasing model is the exponential decay model: EXD.2 and EXD.3.

Examples

## First model
met.as.m1<-drm(gain ~ dose, product, data = methionine, fct = AR.3(), 
pmodels = list(~1, ~factor(product), ~factor(product)))
plot(met.as.m1, log = "", ylim = c(1450, 1800))
summary(met.as.m1)

## Calculating bioefficacy: approach 1
coef(met.as.m1)[5] / coef(met.as.m1)[4] * 100

## Calculating bioefficacy: approach 2
EDcomp(met.as.m1, c(50,50))

## Simplified models
met.as.m2<-drm(gain ~ dose, product, data = methionine, fct = AR.3(), 
pmodels = list(~1, ~1, ~factor(product)))
anova(met.as.m2, met.as.m1)  # simplification not possible

met.as.m3 <- drm(gain ~ dose, product, data = methionine, fct = AR.3(), 
pmodels = list(~1, ~factor(product), ~1))
anova(met.as.m3, met.as.m1)  # simplification not possible

Description

MCPA, 2,4-D, mecorprop and dichorlprop were applied either as technical grades materials (h = 1, 2, 3, 4) or as commercial formulations (herb = 5, 6, 7, 8). Each experimental unit consisted of five 1-week old seedlings grown together in a pot of nutrient solution during 14 days.

Usage

data(auxins)

Format

A data frame with 150 observations on the following 5 variables.

r

a numeric vector

h

a numeric vector

w

a numeric vector

y

a numeric vector

dose

a numeric vector

Details

Data are parts of a larger joint action experiment with various herbicides.

The eight herbicide preparations are naturally grouped into four pairs: (1, 5), (2, 6), (3, 7), and (4, 8), and in each pair of herbicides should have the same active ingredients but different formulation constituents, which were assumed to be biologically inert. The data consist of the 150 observations y of dry weights, each observation being the weight of five plants grown in the same pot. All the eight herbicide preparations have essentially the same mode of action in the plant; they all act like the plant auxins, which are plant regulators that affect cell enlongation an other essential metabolic pathways. One of the objects of the experiment was to test if the response functions were identical except for a multiplicative factor in the dose. This is a necessary, but not a sufficient, condition for a similar mode of action for the herbicides.

Source

Streibig, J. C. (1987). Joint action of root-absorbed mixtures of auxin herbicides in Sinapis alba L. and barley (Hordeum vulgare L.) Weed Research, 27, 337–347.

References

Rudemo, M., Ruppert, D., and Streibig, J. C. (1989). Random-Effect Models in Nonlinear Regression with Applications to Bioassay. Biometrics, 45, 349–362.

Examples

## Fitting model with varying lower limits
auxins.m1 <- boxcox(drm(y ~ dose, h, 
pmodels = data.frame(h, h, 1, h), fct = LL.4(), data = auxins), method = "anova")

## Fitting model with common lower limit
auxins.m2 <- boxcox(drm(y ~ dose, h, 
pmodels = data.frame(h, 1, 1, h), fct = LL.4(), data = auxins), method = "anova")

## Comparing the two models
anova(auxins.m2, auxins.m1)

Description

By inverse regression backfitted dose values are calculated for the mean response per dose.

Usage

backfit(drcObject)

Arguments

Value

Two columns with the original dose values and the corresponding backfitted values using the fitted dose-response model. For extreme dose values (e.g., high dose ) the backfitted values may not be well-defined (see the example below).

Author(s)

Christian Ritz after a suggestion from Keld Sorensen.

References

??

See Also

A related function is ED.drc.

Examples

ryegrass.LL.4 <- drm(rootl~conc, data=ryegrass, fct=LL.4())

backfit(ryegrass.LL.4)

Description

'baro5' allows specification of the baroreflex 5-parameter dose response function, under various constraints on the parameters.

Usage

baro5(fixed = c(NA, NA, NA, NA, NA), names = c("b1", "b2", "c", "d", "e"),
  method = c("1", "2", "3", "4"), ssfct = NULL)

Arguments

Details

The five-parameter function given by the expression

$y = c + \frac{d-c}{1+f\exp(b1(\log(x)-\log(e))) + (1-f)\exp(b2(\log(x)-\log(e)))}$

$f = 1/( 1 + \exp((2b1b2/|b1+b2|)(\log(x)-\log(e))))$

If the difference between the parameters b1 and b2 is different from 0 then the function is asymmetric.

Value

The value returned is a list containing the nonlinear model function, the self starter function and the parameter names.

Note

See the example for the dataset heartrate.

Author(s)

Christian Ritz

References

Ricketts, J. H. and Head, G. A. (1999) A five-parameter logistic equation for investigating asymmetry of curvature in baroreflex studies. Am. J. Physiol. (Regulatory Integrative Comp. Physiol. 46), 277, 441–454.

Description

'BC.4' and 'BC.5' provide the Brain-Cousens modified log-logistic models for describing u-shaped hormesis.

Usage

BC.5(fixed = c(NA, NA, NA, NA, NA), names = c("b", "c", "d", "e", "f"), ...)

  BC.4(fixed = c(NA, NA, NA, NA), names = c("b", "d", "e", "f"), ...)

Arguments

Details

The model function for the Brain-Cousens model (Brain and Cousens, 1989) is

$f(x, b,c,d,e,f) = c + \frac{d-c+fx}{1+\exp(b(\log(x)-\log(e)))}$

,

and it is a five-parameter model, obtained by extending the four-parameter log-logistic model (LL.4 to take into account inverse u-shaped hormesis effects.

The parameters have the following interpretations

• $b$: Not direct interpretation

• $c$: Lower horizontal asymptote

• $d$: Upper horizontal asymptote

• $e$: Not direct interpretation

• $f$: Size of the hormesis effect: the larger the value the larger is the hormesis effect. $f=0$ corresponds to no hormesis effect and the resulting model is the four-parameter log-logistic model. This parameter should be positive in order for the model to make sense.

Fixing the lower limit at 0 yields the four-parameter model

$f(x) = 0 + \frac{d-0+fx}{1+\exp(b(\log(x)-\log(e)))}$

used by van Ewijk and Hoekstra (1993).

Value

See braincousens.

Note

This function is for use with the function drm.

Author(s)

Christian Ritz

References

Brain, P. and Cousens, R. (1989) An equation to describe dose responses where there is stimulation of growth at low doses, Weed Research, 29, 93–96.

van Ewijk, P. H. and Hoekstra, J. A. (1993) Calculation of the EC50 and its Confidence Interval When Subtoxic Stimulus Is Present, Ecotoxicology and Environmental Safety, 25, 25–32.

See Also

More details are found for the general model function braincousens.

Examples

## Fitting the data in van Ewijk and Hoekstra (1993)
lettuce.bcm1 <- drm(weight ~ conc, data = lettuce, fct=BC.5())
modelFit(lettuce.bcm1)
plot(lettuce.bcm1)

lettuce.bcm2 <- drm(weight ~conc, data = lettuce, fct=BC.4())
summary(lettuce.bcm2)
ED(lettuce.bcm2, c(50))  
# compare the parameter estimate and 
# its estimated standard error 
# to the values in the paper by 
# van Ewijk and Hoekstra (1993)


## Brain-Cousens model with the constraint b>3
ryegrass.bcm1 <- drm(rootl ~conc, data = ryegrass, fct = BC.5(), 
lower = c(3, -Inf, -Inf, -Inf, -Inf), control = drmc(constr=TRUE))

summary(ryegrass.bcm1)

## Brain-Cousens model with the constraint f>0 
## (no effect as the estimate of f is positive anyway)
ryegrass.bcm2 <- drm(rootl ~conc, data = ryegrass, fct = BC.5(), 
lower = c(-Inf, -Inf, -Inf, -Inf, 0), control = drmc(constr=TRUE))

summary(ryegrass.bcm2)

Description

Finds the optimal Box-Cox transformation for non-linear regression.

Usage

## S3 method for class 'drc'
boxcox(object, lambda = seq(-2, 2, by = 0.25), plotit = TRUE, bcAdd = 0, 
method = c("ml", "anova"),  level = 0.95, eps = 1/50, 
xlab = expression(lambda), ylab = "log-Likelihood", ...)

Arguments

Details

The optimal lambda value is determined using a profile likelihood approach: For each lambda value the dose-response regression model is fitted and the lambda value (and corresponding model fit) resulting in the largest value of the log likelihood function is chosen.

Value

An object of class "drc" (returned invisibly). If plotit = TRUE a plot of loglik vs lambda is shown indicating a confidence interval (by default 95 the optimal lambda value.

Author(s)

Christian Ritz

References

Carroll, R. J. and Ruppert, D. (1988) Transformation and Weighting in Regression, New York: Chapman and Hall (Chapter 4).

See Also

For linear regression the analogue is boxcox.

Examples

## Fitting log-logistic model without transformation
ryegrass.m1 <- drm(ryegrass, fct = LL.4())
summary(ryegrass.m1)

## Fitting the same model with the optimal Box-Cox transformation
ryegrass.m2 <- boxcox(ryegrass.m1)
summary(ryegrass.m2)

Description

'braincousens' provides a very general way of specifying Brain-Cousens' modified log- logistic model for describing hormesis, under various constraints on the parameters.

Usage

braincousens(fixed = c(NA, NA, NA, NA, NA), 
  names = c("b", "c", "d", "e", "f"),
  method = c("1", "2", "3", "4"), ssfct = NULL,
  fctName, fctText)

Arguments

Details

The Brain-Cousens model is given by the expression

$f(x) = c + \frac{d-c+fx}{1+\exp(b(\log(x)-\log(e)))}$

which is a five-parameter model.

It is a modification of the four-parameter logistic curve to take hormesis into account proposed by Brain and Cousens (1989).

Value

The value returned is a list containing the non-linear function, the self starter function, the parameter names and additional model specific objects.

Note

This function is for use with the function drm.

The convenience functions of braincousens are BC.4 and BC.5. These functions should be used rather than braincousens directly.

Author(s)

Christian Ritz

References

Brain, P. and Cousens, R. (1989) An equation to describe dose responses where there is stimulation of growth at low doses, Weed Research, 29, 93–96.

Description

Bread and meat for the sandwich estimator of the variance-covariance.

Usage

bread.drc(x, ...)
  
  estfun.drc(x, ...)

Arguments

Details

The details are provided by Zeileis (2006).

Value

The unscaled hessian is returned by bread.drc, whereas estfun.drc returns the estimating function evaluated at the data and the parameter estimates.

By default no clustering is assumed, corresponding to robust standard errors under independence. If a cluster variable is provided the log likelihood contributions provided by estfun are summed up for each cluster.

Author(s)

Christian Ritz

References

Zeileis, A. (2006) Object-oriented Computation of Sandwich Estimators, J. Statist. Software, 16, Issue 9.

Examples

## The lines below requires that the packages
## 'lmtest' and 'sandwich' are installed
# library(lmtest)
# library(sandwich)

# ryegrass.m1<-drm(rootl ~ conc, data = ryegrass, fct = LL.4())

# Standard summary output
# coeftest(ryegrass.m1)  

# Output with robust standard errors
# coeftest(ryegrass.m1, vcov = sandwich)

Description

'cedergreen' provides a very general way of specifying then Cedergreen-Ritz-Streibig modified log-logistic model for describing hormesis, under various constraints on the parameters.

CRS.6 is the extension of link{cedergreen} with freely varying alpha parameter.

For u-shaped hormesis data 'ucedergreen' provides a very general way of specifying the Cedergreen-Ritz-Streibig modified log-logistic model, under various constraints on the parameters.

Usage

cedergreen(fixed = c(NA, NA, NA, NA, NA), 
  names = c("b", "c", "d", "e", "f"), 
  method = c("1", "2", "3", "4"), ssfct = NULL,
  alpha, fctName, fctText)
  
  CRS.6(fixed = c(NA, NA, NA, NA, NA, NA), 
  names = c("b","c","d","e","f","g"),
  method = c("1", "2", "3", "4"), ssfct = NULL)
  
  ucedergreen(fixed = c(NA, NA, NA, NA, NA), 
  names = c("b", "c", "d", "e", "f"), 
  method = c("1", "2", "3", "4"), ssfct = NULL,
  alpha)

Arguments

Details

The model is given by the expression

$f(x) = c + \frac{d-c+f exp(-1/(x^{\alpha}))}{1+exp(b(log(x)-log(e)))}$

which is a five-parameter model (alpha is fixed or freely varying). Not all features (eg EC/ED calculation) are available for the model with freely varying alpha.

It is a modification of the four-parameter logistic curve to take hormesis into account.

The u-shaped model is given by the expression

$f(x) = cd - \frac{d-c+f \exp(-1/x^{\alpha})}{1+\exp(b(\log(x)-\log(e)))}$

Value

The value returned is a list containing the non-linear function, the self starter function and the parameter names.

Note

The functions are for use with the functions drm.

Author(s)

Christian Ritz

References

Cedergreen, N. and Ritz, C. and Streibig, J. C. (2005) Improved empirical models describing hormesis, Environmental Toxicology and Chemistry 24, 3166–3172.

See Also

For fixed alpha, several special cases are handled by the following convenience functions CRS.4a, CRS.4b, CRS.4c, CRS.5a, CRS.5b, CRS.5c, UCRS.4a, UCRS.4b, UCRS.4c, UCRS.5a, UCRS.5b, UCRS.5c where a, b and c correspond to the pre-specified alpha values 1, 0.5 and 0.25, respectively.

Examples

## Estimating CRS model with alpha unknown
lettuce.crsm1 <- drm(weight~conc, data = lettuce, fct = CRS.6())
summary(lettuce.crsm1)
plot(lettuce.crsm1)  # oops: not increasing until hormesis peak

Description

Germination data from tests of chickweed seeds from chlorsulfuron resistant and sensitive biotypes

Usage

data(chickweed)

Format

A data frame with 35 observations on the following 3 variables.

start

a numeric vector of left endpoints of the monitoring intervals

end

a numeric vector of right endpoints of the monitoring intervals

count

a numeric vector of the number of seeds germinated in the interval between start and end

time

a numeric vector of the non-zero left endpoints of the monitoring intervals (often used for recording in practice)

Details

The germination tests of chickweed seeds from chlorsulfuron resistant and sensitive biotypes in central Zealand were done in petri dishes (diameter: 9.0cm) in a dark growth cabinet at a temperature of 5 degrees Celsius. The seeds were incubated for 24 hours in a 0.3% solution of potassium nitrate in order to imbibe seeds prior to the test. A total of 200 seeds were placed on filter plate. After initialization of the tests, the number of germinated seeds was recorded and removed at 34 consecutive inspection times. Definition of a germinated seed was the breakthrough of the seed testa by the radicle.

Chickweed is known to have dormant seeds and therefore we would not expect 100% germination. It means that the upper limit of the proportion germinated has to be incorporated as a parameter into a model, which adequately reflects the experimental design as well as any expectations about the resulting outcome.

Source

Data are kindly provided by Lisa Borggaard (formerly at the Faculty of Life Sciences, University of Copenhagen).

References

Ritz, C., Pipper, C. B. and Streibig, J. C. (2013) Analysis of germination data from agricultural experiments, Europ. J. Agronomy, 45, 1–6.

Examples

## Incorrect analysis using a logistic regression model
## (treating event times as binomial data)
## The argument "type" specifies that binomial data are supplied 
chickweed.m0a <- drm(count/200 ~ time, weights = rep(200, 34), 
data = chickweed0, fct = LL.3(), type = "binomial")
summary(chickweed.m0a) # showing a summmary of the model fit (including parameter estimates)

## Incorrect analysis based on nonlinear regression
## LL.3() refers to the three-parameter log-logistic model
## As the argument "type" is not specified it is assumed that the data type
##  is continuous and nonlinear regression based on least squares estimation is carried out
chickweed.m0b <- drm(count/200 ~ time, data = chickweed0, fct = LL.3())
summary(chickweed.m0b)  # showing a summmary of the model fit (including parameter estimates)

## How to re-arrange the data for fitting the event-time model
## (only for illustration of the steps needed for converting a dataset, 
##  but in this case not needed as both datasets are already provided in "drc")
#chickweed <- data.frame(start = c(0, chickweed0$time), end = c(chickweed0$time, Inf)) 
#chickweed$count <- c(0, diff(chickweed0$count), 200 - tail(chickweed0$count, 1))
#head(chickweed)  # showing top 6 lines of the dataset
#tail(chickweed)  # showing bottom 6 lines

## Fitting the event-time model (by specifying the argument type explicitly)
chickweed.m1 <- drm(count~start+end, data = chickweed, fct = LL.3(), type = "event")
summary(chickweed.m1)  # showing a summmary of the model fit (including parameter estimates)

## Summary output with robust standard errors
## library(lmtest)
## library(sandwich)
## coeftest(chickweed.m1, vcov = sandwich)

## Calculating t10, t50, t90 for the distribution of viable seeds
ED(chickweed.m1, c(10, 50, 90))

## Plotting data and fitted regression curve
plot(chickweed.m1, xlab = "Time (hours)", ylab = "Proportion germinated", 
xlim=c(0, 340), ylim=c(0, 0.25), log="", lwd=2, cex=1.2)  
## Adding the fitted curve obtained using nonlinear regression
plot(chickweed.m0b, add = TRUE, lty = 2, xlim=c(0, 340), 
ylim=c(0, 0.25), log="", lwd=2, cex=1.2)
# Note: the event-time model has slightly better fit at the upper limit

## Enhancing the plot (to look like in the reference paper)
abline(h = 0.20011, lty = 3, lwd = 2)
text(-15, 0.21, "Upper limit: d", pos = 4, cex = 1.5)

segments(0,0.1,196,0.1, lty = 3, lwd = 2)
segments(196,0.1, 196, -0.1, lty = 3, lwd = 2)
text(200, -0.004, expression(paste("50% germination: ", t[50])), pos = 4, cex = 1.5)

abline(a = 0.20011/2-0.20011*20.77/4, b = 0.20011*20.77/4/196, lty = 3, lwd = 2)
#text(200, 0.1, expression(paste("Slope: ", b*(-d/(4*t[50])))), pos = 4, cex = 1.5)
text(200, 0.1, expression("Slope: b" %.% "constant"), pos = 4, cex = 1.5)
points(196, 0.1, cex = 2, pch = 0)


## Adding confidence intervals

## Predictions from the event-time model
#coefVec <- coef(chickweed.m1)
#names(coefVec) <- c("b","d","e")
#
#predFct <- function(tival)
#{
#    as.numeric(deltaMethod(coefVec, paste("d/(1+exp(b*(log(",tival,")-log(e))))"), 
#    vcov(chickweed.m1)))
#}
#predFctv <- Vectorize(predFct, "tival")
#
#etpred <- t(predFctv(0:340))
#lines(0:340, etpred[,1]-1.96*etpred[,2], lty=1, lwd=2, col="darkgray")
#lines(0:340, etpred[,1]+1.96*etpred[,2], lty=1, lwd=2, col="darkgray")
#
### Predictions from the nonlinear regression model
#nrpred <- predict(chickweed.m0b, data.frame(time=0:340), interval="confidence")
#lines(0:340, nrpred[,2], lty=2, lwd=2, col="darkgray")
#lines(0:340, nrpred[,3], lty=2, lwd=2, col="darkgray")

Description

For single mixture data combination indices for effective doses as well as effects may be calculated and visualized.

Usage

CIcomp(mixProp, modelList, EDvec)

CIcompX(mixProp, modelList, EDvec, EDonly = FALSE)

plotFACI(effList, indAxis = c("ED", "EF"), caRef = TRUE, 
showPoints = FALSE, add = FALSE, ylim, ...)

Arguments

Details

CIcomp calculates the classical combination index for effective doses whereas CIcompX calculates the combination index also for effects as proposed by Martin-Betancor et al. (2015); for details and examples using "drc" see the supplementary material of this paper. The function plotFACI may be used to visualize the calculated combination index as a function of the fraction affected.

Value

CIcomp returns a matrix which one row per ED value. Columns contain estimated combination indices, their standard errors and 95% confidence intervals, p-value for testing CI=1, estimated ED values for the mixture data and assuming concentration addition (CA) with corresponding standard errors.

CIcompX returns similar output both for effective doses and effects (as a list of matrices).

Author(s)

Christian Ritz and Ismael Rodea-Palomares

References

Martin-Betancor, K. and Ritz, C. and Fernandez-Pinas, F. and Leganes, F. and Rodea-Palomares, I. (2015) Defining an additivity framework for mixture research in inducible whole-cell biosensors, Scientific Reports 17200.

See Also

See mixture for simultaneous modelling of several mixture ratios, but only at the ED50 level.

See also the help page for metals.

Examples

## Fitting marginal models for the 2 pure substances
acidiq.0 <- drm(rgr ~ dose, data = subset(acidiq, pct == 999 | pct == 0), fct = LL.4())
acidiq.100 <- drm(rgr ~ dose, data = subset(acidiq, pct == 999 | pct == 100), fct = LL.4())

## Fitting model for single mixture with ratio 17:83
acidiq.17 <- drm(rgr ~ dose, data = subset(acidiq, pct == 17 | pct == 0), fct = LL.4())
    
## Calculation of combination indices based on ED10, ED20, ED50
CIcomp(0.17, list(acidiq.17, acidiq.0, acidiq.100), c(10, 20, 50))  
## CI>1 significantly for ED10 and ED20, but not so for ED50

Description

Extract parameter estimates.

Usage

## S3 method for class 'drc'
coef(object, ...)

Arguments

Value

A vector of parameter coefficients which are extracted from the model object 'object'.

Note

This function may work even in cases where 'summary' does not, because the parameter estimates are retrieved directly from the model fit object without any additional computations of summary statistics and standard errors.

Author(s)

Christian Ritz

See Also

A more comprehensive summary is obtained using summary.drc.

Examples

## Fitting a four-parameter log-logistic model
ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())
coef(ryegrass.m1)

Description

Comparison of a pair of effective dose values from independent experiments where only the estimates and their standard errors are reported.

Usage

comped(est, se, log = TRUE, interval = TRUE, operator = c("-", "/"), 
  level = 0.95, df = NULL)

Arguments

Details

The choice "/" for the argument operator and FALSE for log will result in estimation of a socalled relative potency (sometimes also called a selectivity index).

The combination TRUE for log and "/" for operator only influences the confidence interval, that is no ratio is calculated based on logarithm-transformed effective dose values.

By default confidence interval relies on percentiles in the normal distribution.

In case the entire dataset is available the functions drm and (subsequently) EDcomp should be used instead.

Value

A matrix with the estimated difference or ratio and the associated standard error and the resulting confidence interval (unless not requested).

Note

The development of the function comped is a side effect of the project on statistical analysis of toxicity data funded by the Danish EPA ("Statistisk analyse og biologisk tolkning af toksicitetsdata", MST j.nr. 669-00079).

Author(s)

Christian Ritz

References

Wheeler, M. W. and Park, R. M. and Bailer, A. J. (2006) Comparing median lethal concentration values using confidence interval overlap or ratio tests, Environmental Toxicology and Chemistry, 25, 1441–1441.

See Also

The function ED.drc calculates arbitrary effective dose values based on a model fit. The function EDcomp calculates relative potencies based on arbitrary effective dose values.

Examples

## Fitting the model
S.alba.m1 <- boxcox(drm(DryMatter~Dose, Herbicide, data=S.alba, fct = LL.4(),
pmodels=data.frame(Herbicide,1,1,Herbicide)), method = "anova") 

## Displaying estimated ED values
ED(S.alba.m1, c(10, 90))

## Making comparisons of ED50 in two ways and for both differences and ratios
compParm(S.alba.m1, "e", "/")

comped(c(28.396147, 65.573335), c(1.874598, 5.618945), log=FALSE, operator = "/")
# similar result

compParm(S.alba.m1, "e", "-")

comped(c(28.396147, 65.573335), c(1.874598, 5.618945), log=FALSE, operator = "-")
# similar result

## Making comparisons of ED10 and ED90
comped(c(21.173, 44.718), c(11.87, 8.42), log=FALSE, operator = "/")

comped(c(21.173, 44.718), c(11.87, 8.42), log=FALSE, operator = "/", interval = FALSE)

comped(c(21.173, 44.718), c(11.87, 8.42), log=FALSE, operator = "-")

Description

Compare parameters from different assays, either by means of ratios or differences.

Usage

compParm(object, strVal, operator = "/", vcov. = vcov, od = FALSE, 
  pool = TRUE, display = TRUE)

Arguments

Details

The function compares actual parameter estimates, and therefore the results depend on the parameterisation used. Probably it is most useful in combination with the argument collapse in drm for specifying parameter constraints in models, either through data frames or lists with formulas without intercept (-1).

Value

A matrix with columns containing the estimates, estimated standard errors, values of t-statistics and p-values for the null hypothesis that the ratio equals 1 or that the difference equals 0 (depending on the operator argument).

Author(s)

Christian Ritz

Examples

# Fitting a model with names assigned to the parameters!
spinach.m1 <- drm(SLOPE~DOSE, CURVE, data = spinach,
fct = LL.4(names = c("b", "lower", "upper", "ed50")))

## Calculating ratios of parameter estimates for the parameter named "ed50"
compParm(spinach.m1, "ed50")      

## Calculating differences between parameter estimates for the parameter named "ed50"
compParm(spinach.m1, "ed50", "-")

Description

Computes confidence intervals for one or more parameters in a model of class 'drc'.

Usage

## S3 method for class 'drc'
confint(object, parm, level = 0.95, pool = TRUE, ...)

Arguments

Details

For binomial and Poisson data the confidence intervals are based on the normal distribution, whereas t distributions are used of for continuous/quantitative data.

Value

A matrix (or vector) with columns giving lower and upper confidence limits for each parameter. These will be labelled as (1-level)/2 and 1 - (1-level)/2 in

Author(s)

Christian Ritz

Examples

## Fitting a four-parameter log-logistic model
ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())

## Confidence intervals for all parameters
confint(ryegrass.m1)

## Confidence interval for a single parameter
confint(ryegrass.m1, "e")

Description

'CRS.4a', 'CRS.4b' and 'CRS.4c' provide the Cedergreen-Ritz-Streibig modified log-logistic model for describing hormesis with the lower limit equal to 0.

'UCRS.4a', 'UCRS.4b' and 'UCRS.4c' provide the Cedergreen-Ritz-Streibig modified log-logistic model for describing u-shaped hormesis with the lower limit equal to 0.

Usage

CRS.4a(names = c("b", "d", "e", "f"), ...)

  UCRS.4a(names = c("b", "d", "e", "f"), ...)

Arguments

Details

The model is given by the expression

$f(x) = 0 + \frac{d-0+f \exp(-1/x)}{1+\exp(b(\log(x)-\log(e)))}$

which is a five-parameter model.

It is a modification of the four-parameter logistic curve to take hormesis into account.

The u-shaped model is given by the expression

$f(x) = 0 + d - \frac{d-0+f \exp(-1/x^{\alpha})}{1+\exp(b(\log(x)-\log(e)))}$

The a,b,c models are obtained by setting alpha equal to 1, 0.5 and 0.25, respectively.

Value

See cedergreen.

Note

This function is for use with the function drm.

Author(s)

Christian Ritz

References

See the reference under cedergreen.

See Also

Similar functions are CRS.5a and UCRS.5a, but with an extra parameter for the lower limit.

Examples

## Fitting modified logistic models
lettuce.crsm1 <- drm(lettuce[,c(2,1)], fct=CRS.4a())
summary(lettuce.crsm1)
ED(lettuce.crsm1, c(50))

## Need to explicitly specify that the upper limit
## is the reference in order to get ED10 and ED90 right
ED(lettuce.crsm1, c(10, 50, 90), reference = "upper")

lettuce.crsm2 <- drm(lettuce[,c(2,1)], fct=CRS.4b())
summary(lettuce.crsm2)
ED(lettuce.crsm2, c(50))

lettuce.crsm3 <- drm(lettuce[,c(2,1)], fct=CRS.4c())
summary(lettuce.crsm3)
ED(lettuce.crsm3, c(50))

Description

'CRS.5a', 'CRS.5b' and 'CRS.5c' provide the Cedergreen-Ritz-Streibig modified log-logistic model for describing (inverse u-shaped or j-shaped) hormesis.

'UCRS.5a', 'UCRS.5b' and 'UCRS.5c' provide the Cedergreen-Ritz-Streibig modified log-logistic model for describing u-shaped hormesis.

Usage

CRS.5a(names = c("b", "c", "d", "e", "f"), ...)

  UCRS.5a(names = c("b", "c", "d", "e", "f"), ...)

Arguments

Details

The model function for inverse u-shaped hormetic patterns is

$f(x) = c + \frac{d-c+f \exp(-1/x^{\alpha})}{1+\exp(b(\log(x)-\log(e)))}$

,

which is a five-parameter model. It is a modification of the four-parameter log-logistic curve to take hormesis into account.

The parameters have the following interpretations

• $b$: Not direct interpretation

• $c$: Lower horizontal asymptote

• $d$: Upper horizontal asymptote

• $e$: Not direct interpretation

• $f$: Size of the hormesis effect: the larger the value the larger is the hormesis effect. $f=0$ corresponds to no hormesis effect and the resulting model is the four-parameter log-logistic model. This parameter should be positive in order for the model to make sense.

The model function for u-shaped hormetic patterns is

$f(x) = c + d - \frac{d-c+f \exp(-1/x^{\alpha})}{1+\exp(b(\log(x)-\log(e)))}$

This model also simplifies to the four-parameter log-logistic model in case $f=0$ (in a slightly different parameterization as compared to the one used in LL.4).

The models denoted a,b,c are obtained by fixing the alpha parameter at 1, 0.5 and 0.25, respectively.

Value

See cedergreen.

Note

This function is for use with the function drm.

Author(s)

Christian Ritz

References

See the reference under cedergreen.

See Also

Similar functions are CRS.4a and UCRS.4a, but with the lower limit (the parameter $c$) fixed at 0 (one parameter less to be estimated).

Examples

## Modified logistic model
lettuce.m1 <- drm(lettuce[,c(2,1)], fct=CRS.5a())
summary(lettuce.m1)
ED(lettuce.m1, c(50))

lettuce.m2 <- drm(lettuce[,c(2,1)], fct=CRS.5b())
summary(lettuce.m2)
ED(lettuce.m2, c(50))

lettuce.m3 <- drm(lettuce[,c(2,1)], fct=CRS.5c())
summary(lettuce.m3)
ED(lettuce.m3, c(50))

Description

The number of immobile daphnids –in contrast to mobile daphnids– out of a total of 20 daphnids was counted for several concentrations of a toxic substance.

Usage

data(daphnids)

Format

A data frame with 16 observations on the following 4 variables.

dose

a numeric vector

no

a numeric vector

total

a numeric vector

time

a factor with levels 24h 48h

Details

The same daphnids were counted at 24h and later again at 48h.

Source

Nina Cedergreen, Faculty of Life Sciences, University of Copenhagen, Denmark.

Examples

## Fitting a model with different parameters
## for different curves
daphnids.m1 <- drm(no/total~dose, time, weights = total, 
data = daphnids, fct = LL.2(), type = "binomial")

## Goodness-of-fit test
modelFit(daphnids.m1)

## Summary of the data
summary(daphnids.m1)

## Fitting a model with a common intercept parameter
daphnids.m2 <- drm(no/total~dose, time, weights = total, 
data = daphnids, fct = LL.2(), type = "binomial", 
pmodels = list(~1, ~time))

Description

The two decontaminants 1-hexadecylpyridium chloride and oxalic acid were used. Additionally there was a control group (coded as concentration 0 and only included under oxalic acid).

Usage

data("decontaminants")

Format

A data frame with 128 observations on the following 3 variables.

conc

a numeric vector of percentage weight per volume

count

a numeric vector of numbers of M. bovis colonies at stationarity

group

a factor with levels hpc and oxalic of the decontaminants used

Details

These data examplify Wadley's problem: counts where the maximum number is not known. The data were analyzed by Trajstman (1989) using a three-parameter logistic model and then re-analyzed by Morgan and Smith (1992) using a three-parameter Weibull type II model. In both cases the authors adjusted for overdispersion (in different ways).

It seems that Morgan and Smith (1992) fitted separate models for the two decontaminants and using the control group for both model fits. In the example below a joint model is fitted where the control group is used once to determine a shared upper limit at concentration 0.

Source

Trajstman, A. C. (1989) Indices for Comparing Decontaminants when Data Come from Dose-Response Survival and Contamination Experiments, Applied Statistics, 38, 481–494.

References

Morgan, B. J. T. and Smith, D. M. (1992) A Note on Wadley's Problem with Overdispersion, Applied Statistics, 41, 349–354.

Examples

## Wadley's problem using a three-parameter log-logistic model
decon.LL.3.1 <- drm(count~conc, group, data = decontaminants, fct = LL.3(), 
type = "Poisson", pmodels = list(~group, ~1, ~group))

summary(decon.LL.3.1)

plot(decon.LL.3.1)


## Same model fit in another parameterization (no intercepts)
decon.LL.3.2 <- drm(count~conc, group, data = decontaminants, fct=LL.3(), 
type = "Poisson", pmodels = list(~group-1, ~1, ~group-1))

summary(decon.LL.3.2)

Description

Quantal assay data from an experiment where the insectide deguelin was applied to Macrosiphoniella sanborni.

Usage

data(deguelin)

Format

A data frame with 6 observations on the following 4 variables.

dose

a numeric vector of doses applied

log10dose

a numeric vector of logarithm-transformed doses

r

a numeric vector contained number of dead insects

n

a numeric vector contained the total number of insects

Details

The log-logistic model provides an inadequate fit.

The dataset is used in Nottingham and Birch (2000) to illustrate a semiparametric approach to dose-response modelling.

Source

Morgan, B. J. T. (1992) Analysis of Quantal Response Data, London: Chapman \& Hall/CRC (Table 3.9, p. 117).

References

Notttingham, Q. J. and Birch, J. B. (2000) A semiparametric approach to analysing dose-response data, Statist. Med., 19, 389–404.

Examples

## Log-logistic fit
deguelin.m1 <- drm(r/n~dose, weights=n, data=deguelin, fct=LL.2(), type="binomial")
modelFit(deguelin.m1)
summary(deguelin.m1)

## Loess fit
deguelin.m2 <- loess(r/n~dose, data=deguelin, degree=1)

## Plot of data with fits superimposed
plot(deguelin.m1, ylim=c(0.2,1))
lines(1:60, predict(deguelin.m2, newdata=data.frame(dose=1:60)), col = 2, lty = 2)

lines(1:60, 0.95*predict(deguelin.m2, 
newdata=data.frame(dose=1:60))+0.05*predict(deguelin.m1, newdata=data.frame(dose=1:60), se = FALSE),
col = 3, lty=3)

Description

A general model fitting function for analysis of concentration/dose/time-effect/response data.

Usage

drm(formula, curveid, pmodels, weights, data = NULL, subset, fct, 
  type = c("continuous", "binomial", "Poisson", "quantal", "event"), 
  bcVal = NULL, bcAdd = 0, 
  start, na.action = na.omit, robust = "mean", logDose = NULL, 
  control = drmc(), lowerl = NULL, upperl = NULL, separate = FALSE, pshifts = NULL)

Arguments

Details

This function relies on the general optimiser function optim for the minimisation of negative log likelihood function. For a continuous response this reduces to least squares estimation, which is carried out by minimising the following sums of squares

$\sum_{i=1}^N [w_i (y_i-f_i)]^2$

where $y_i$, $f_i$, and $w_i$ correspond to the observed value, expected value, and the weight respectively, for the ith observation (from 1 to $N$).

The control arguments are specified using the function drmc.

Setting lowerl and/or upperl automatically invokes constrained optimisation.

The columns of a data frame argument to pmodels are automatically converted into factors. This does not happen if a list is specified.

Value

An object of class 'drc'.

Note

For robust estimation MAD (median abslolute deviance) is used to estimate the residual variance.

Author(s)

Christian Ritz and Jens C. Streibig

See Also

Examples using drm found in the help pages of ryegrass (continuous data), secalonic (continuous data), and selenium (binomial data), as well as for a number of other datasets and functions in drc.

Description

Set control arguments in the control argument in the function 'drm'.

Usage

drmc(constr = FALSE, errorm = TRUE, maxIt = 500, method="BFGS", 
  noMessage = FALSE, relTol = 1e-07, rmNA=FALSE, useD = FALSE, 
  trace = FALSE, otrace = FALSE, warnVal = -1, dscaleThres = 1e-15, rscaleThres = 1e-15)

Arguments

Value

A list with 8 components, one for each of the above arguments.

Note

The use of a non-zero constant bcAdd may in some cases make it more difficult to obtain convergence of the estimation procedure.

Author(s)

Christian Ritz

Examples

### Displaying the default settings
drmc()  

### Using 'method' argument
model1 <- drm(ryegrass, fct = LL.4())

model2 <- drm(ryegrass, fct = LL.4(),
control = drmc(method = "Nelder-Mead"))

Description

The dataset was obtained from a toxicity test using earthworms, and it contains the number of earthworms remaining in a container that was contaminated with a toxic substance (not disclosed) at various doses; so the number of earthworms not migrating to the neighbouring uncontaminated container.

Usage

data(earthworms)

Format

A data frame with 35 observations on the following 3 variables.

dose

a numeric vector of dose values

number

a numeric vector containing counts of remaining earthworms in the container

total

a numeric vector containing total number of earthworms put in the containers

Details

At dose 0 around half of the earthworms is expected be in each of the two containers. Thus it is not appropriate to fit an ordinary logistic regression with log(dose) as explanatory variable to these data as it implies an upper limit of 1 at dose 0 and in fact this model does not utilise the observations at dose 0 (see the example section below).

Source

The dataset is kindly provided by Nina Cedergreen, Faculty of Life Sciences, University of Copenhagen, Denmark.

Examples

## Fitting a logistic regression model
earthworms.m1 <- drm(number/total~dose, weights = total, data = earthworms,
fct = LL.2(), type = "binomial")
modelFit(earthworms.m1)  # a crude goodness-of-fit test

## Fitting an extended logistic regression model 
##  where the upper limit is estimated
earthworms.m2 <- drm(number/total~dose, weights = total, data = earthworms,
fct = LL.3(), type = "binomial")
modelFit(earthworms.m2)  # goodness-of-fit test
# improvement not visible in test!!!

## Comparing model1 and model2 
## (Can the first model be reduced to the second model?)
anova(earthworms.m1, earthworms.m2)

Description

ED estimates effective doses (ECp/EDp/ICp) for given reponse levels.

Usage

## S3 method for class 'drc'
ED(object, respLev, interval = c("none", "delta", "fls", "tfls"), 
  clevel = NULL, level = ifelse(!(interval == "none"), 0.95, NULL),
  reference = c("control", "upper"), type = c("relative", "absolute"), lref, uref,
  bound = TRUE, od = FALSE, vcov. = vcov, display = TRUE, pool = TRUE, logBase = NULL, 
  multcomp = FALSE, ...)

Arguments

Details

For hormesis models (braincousens and cedergreen), the additional arguments lower and upper may be supplied. These arguments specify the lower and upper limits of the bisection method used to find the ED values. The lower and upper limits need to be smaller/larger than the EDx level to be calculated. The default limits are 0.001 and 1000 for braincousens and 0.0001 and 10000 for cedergreen and ucedergreen, but this may need to be modified (for cedergreen the upper limit may need to be increased and for ucedergreen the lower limit may need to be increased). Note that the lower limit should not be set to 0 (use instead something like 1e-3, 1e-6, ...).

Value

An invisible matrix containing the shown matrix with two or more columns, containing the estimates and the corresponding estimated standard errors and possibly lower and upper confidence limits. Or, alternatively, a list with elements that may be plugged directly into parm in the package multcomp (in case the argument multcomp is TRUE).

Author(s)

Christian Ritz

See Also

backfit, isobole, and maED use ED for specific calculations involving estimated ED values.

The related function EDcomp may be used for estimating differences and ratios of ED values, whereas compParm may be used to compare other model parameters.

Examples

## Fitting 4-parameter log-logistic model
ryegrass.m1 <- drm(ryegrass, fct = LL.4())

## Calculating EC/ED values
ED(ryegrass.m1, c(10, 50, 90)) 
## first column: the estimates of ED10, ED50 and ED90
## second column: the corresponding estimated standard errors 

### How to use the argument 'ci'

## Also displaying 95% confidence intervals
ED(ryegrass.m1, c(10, 50, 90), interval = "delta")

## Comparing delta method and back-transformed 
##  confidence intervals for ED values

## Fitting 4-parameter log-logistic 
##  in different parameterisation (using LL2.4)
ryegrass.m2 <- drm(ryegrass, fct = LL2.4())  

ED(ryegrass.m1, c(10, 50, 90), interval = "fls")
ED(ryegrass.m2, c(10, 50, 90), interval = "delta")


### How to use the argument 'bound'

## Fitting the Brain-Cousens model
lettuce.m1 <- drm(weight ~ conc, 
data = lettuce, fct = BC.4())

### Calculating ED[-10]

# This does not work
#ED(lettuce.m1, -10)  

## Now it does work
ED(lettuce.m1, -10, bound = FALSE)  # works
ED(lettuce.m1, -20, bound = FALSE)  # works

## The following does not work for another reason: ED[-30] does not exist 
#ED(lettuce.m1, -30, bound = FALSE)

Description

Relative potencies (also called selectivity indices) for arbitrary doses are compared between fitted dose-response curves.

Usage

EDcomp(object, percVec, percMat = NULL, compMatch = NULL, od = FALSE, vcov. = vcov, 
  reverse = FALSE, 
  interval = c("none", "delta", "fieller", "fls"), 
  level = ifelse(!(interval == "none"), 0.95, NULL),
  reference = c("control", "upper"),  
  type = c("relative", "absolute"),
  display = TRUE, pool = TRUE, logBase = NULL, 
  multcomp = FALSE, ...)
  
  relpot(object, plotit = TRUE, compMatch = NULL, percVec = NULL, interval = "none", 
  type = c("relative", "absolute"), 
  scale = c("original", "percent", "unconstrained"), ...)

Arguments

Details

The function relpot is a convenience function, which is useful for assessing how the relative potency changes as a function of the response level (e.g., for plotting as outlined by Ritz et al (2006)).

Fieller's theorem is incorporated using the formulas provided by Kotz and Johnson (1983) and Finney (1978).

For objects of class 'braincousens' or 'mlogistic' the additional argument may be the 'upper' argument or the 'interval' argument. The 'upper' argument specifies the upper limit of the bisection method. The upper limits needs to be larger than the EDx level to be calculated. The default limit is 1000. The 'interval' argument should specify a rough interval in which the dose yielding the maximum hormetical response lies. The default interval is 'c(0.001, 1000)'. Notice that the lower limit should not be set to 0 (use something like 1e-3, 1e-6, ...).

Value

An invisible matrix containing the shown matrix with two or more columns, containing the estimates and the corresponding estimated standard errors and possibly lower and upper confidence limits. Or, alternatively, a list with elements that may be plugged directly into parm in the package multcomp (in case the argument multcomp is TRUE).

Note

This function only works for the following built-in functions available in the package drc: braincousens, cedergreen, ucedergreen, llogistic, and weibull1.

Author(s)

Christian Ritz

References

Finney, D. J. (1978) Statistical method in Biological Assay, London: Charles Griffin House, 3rd edition (pp. 80–82).

Kotz, S. and Johnson, N. L. (1983) Encyclopedia of Statistical Sciences Volume 3, New York: Wiley \& Sons (pp. 86–87).

Ritz, C. and Cedergreen, N. and Jensen, J. E. and Streibig, J. C. (2006) Relative potency in nonsimilar dose-response curves, Weed Science, 54, 407–412.

See Also

A related function is ED.drc (used for calculating effective doses).

Examples

spinach.LL.4 <- drm(SLOPE~DOSE, CURVE, data = spinach, fct = LL.4())

EDcomp(spinach.LL.4, c(50,50))
EDcomp(spinach.LL.4, c(10,50))
EDcomp(spinach.LL.4, c(10,50), reverse = TRUE)

## Using the package multcomp
#sires <- SI(spinach.LL.4, c(25, 50, 75))
#library(multcomp)
#summary(glht(parm(sires[[2]][[1]], sires[[2]][[2]]), rhs = 1))

## Comparing specific ratios: 25/25, 50/50, 75/75
#sires2 <- SI(spinach.LL.4, c(25, 50, 75), matrix(c(1, 1, 2, 2, 3, 3), 3, 2, byrow = TRUE))
#library(multcomp)
#summary(glht(parm(sires2[[2]][[1]], sires2[[2]][[2]]), rhs = 1))


## Relative potency of two herbicides
m2 <- drm(DryMatter~Dose, Herbicide, 
data = S.alba, fct = LL.3())

EDcomp(m2, c(50, 50))  
EDcomp(m2, c(50, 50), interval = "delta")
EDcomp(m2, c(50, 50), interval = "fieller")

## Comparison based on an absolute
##  response level

m3 <- drm(SLOPE~DOSE, CURVE,
data = spinach, fct = LL.4())

EDcomp(m3, c(0.5,0.5), compMatch = c(2,4), type = "absolute", interval = "fieller")

EDcomp(m3, c(55,80), compMatch = c(2,4))  
# same comparison using a relative response level


## Relative potency transformed from log scale
m4 <- drm(drymatter~log(dose), treatment, data=G.aparine[-c(1:40), ], 
pmodels = data.frame(treatment,treatment,1,treatment), fct = LL2.4())

EDcomp(m4, c(50,50), interval = "fls", logBase = exp(1))

Description

Relative growth rate in biomass of mixed sewage microorganisms (per hour) as a function of increasing concentrations of the antibiotic erythromycin (mg/l).

Usage

data(etmotc)

Format

A data frame with 57 observations on the following 4 variables.

cell

a numeric vector

dose1

a numeric vector

pct1

a numeric vector

rgr1

a numeric vector

Details

Data stem from an experiment investigating the effect of pharmaceuticals, that are used in human and veterinary medicine and that are being released into the aquatic environment through waste water or through manure used for fertilising agricultural land. The experiment constitutes a typical dose-response situation. The dose is concentration of the antibiotic erythromycin (mg/l), which is an antibiotic that can be used by persons or animals showing allergy to penicillin, and the measured response is the relative growth rate in biomass of mixed sewage microorganisms (per hour), measured as turbidity two hours after exposure by means of a spectrophotometer. The experiment was designed in such a way that eight replicates were assigned to the control (dose 0), but no replicates were assigned to the 7 non-zero doses. Further details are found in Christensen et al (2006).

Source

Christensen, A. M. and Ingerslev, F. and Baun, A. 2006 Ecotoxicity of mixtures of antibiotics used in aquacultures, Environmental Toxicology and Chemistry, 25, 2208–2215.

Examples

etmotc.m1<-drm(rgr1~dose1, data=etmotc[1:15,], fct=LL.4())
plot(etmotc.m1)
modelFit(etmotc.m1)
summary(etmotc.m1)

etmotc.m2<-drm(rgr1~dose1, data=etmotc[1:15,], fct=W2.4())
plot(etmotc.m2, add = TRUE)
modelFit(etmotc.m2)
summary(etmotc.m2)

etmotc.m3<-drm(rgr1~dose1, data=etmotc[1:15,], fct=W2.3())
plot(etmotc.m3, add = TRUE)
modelFit(etmotc.m3)
summary(etmotc.m3)

Description

Exponential decay model with or without a nonzero lower limit.

Usage

EXD.2(fixed = c(NA, NA), names = c("d", "e"), ...)

  EXD.3(fixed = c(NA, NA, NA), names = c("c", "d", "e"), ...)

Arguments

Details

The exponential decay model is a three-parameter model with mean function:

$f(x) = c + (d-c)(\exp(-x/e))$

The parameter init is the upper limit (attained at $x=0$), the parameter plateau is the lower limit reached for x going to infinity and the parameter $e>0$ is determining the steepness of the decay. The curve is monotonously decreasing in $x$.

Value

A list of class drcMean, containing the mean function, the self starter function, the parameter names and other components such as derivatives and a function for calculating ED values.

Author(s)

Christian Ritz

References

Organisation for Economic Co-operation and Development (OECD) (2006) Current approaches in the statistical analysis of ecotoxicity data: A guidance to application - annexes, Paris: OECD (p. 80).

See Also

Similar models giving exponential increasing curves are AR.2 and AR.3.

Examples

## Fitting an exponential decay model
ryegrass.m1<-drm(rootl~conc, data=ryegrass, fct=EXD.3())

plot(ryegrass.m1)

summary(ryegrass.m1)

Description

For each of six concentration of an insecticid the number of insects affected (out of the number of insects) was recorded.

Usage

data(finney71)

Format

A data frame with 6 observations on the following 3 variables.

dose

a numeric vector

total

a numeric vector

affected

a numeric vector

Source

Finney, D. J. (1971) Probit Analysis, Cambridge: Cambridge University Press.

Examples

## Model with ED50 as a parameter
finney71.m1 <- drm(affected/total ~ dose, weights = total,
data = finney71, fct = LL.2(), type = "binomial")

summary(finney71.m1)
plot(finney71.m1, broken = TRUE, bp = 0.1, lwd = 2)

ED(finney71.m1, c(10, 20, 50), interval = "delta", reference = "control")

## Model fitted with 'glm'
#fitl.glm <- glm(cbind(affected, total-affected) ~ log(dose),
#family=binomial(link = logit), data=finney71[finney71$dose != 0, ])
#summary(fitl.glm)  # p-value almost agree for the b parameter
#
#xp <- dose.p(fitl.glm, p=c(0.50, 0.90, 0.95))  # from MASS
#xp.ci <- xp + attr(xp, "SE") %*% matrix(qnorm(1 - 0.05/2)*c(-1,1), nrow=1)
#zp.est <- exp(cbind(xp.ci[,1],xp,xp.ci[,2]))
#dimnames(zp.est)[[2]] <- c("zp.lcl","zp","zp.ucl")
#zp.est  # not far from above results with 'ED'

## Model with log(ED50) as a parameter
finney71.m2 <- drm(affected/total ~ dose, weights = total,
data = finney71, fct = LL2.2(), type = "binomial")

## Confidence intervals based on back-transformation
##  complete agreement with results based on 'glm'
ED(finney71.m2, c(10, 20, 50), interval = "fls", reference = "control")

Description

Extracts fitted values from an object of class 'drc'.

Usage

## S3 method for class 'drc'
fitted(object, ...)

Arguments

Value

Fitted values extracted from 'object'.

Author(s)

Christian Ritz

Examples

ryegrass.m1 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())
plot(fitted(ryegrass.m1), residuals(ryegrass.m1))  # a residual plot

Description

Model function for specifying dose-response models that are a combination of a logistic model and an appropriate class of fractional polynomials.

Usage

fplogistic(p1, p2, fixed = c(NA, NA, NA, NA), names = c("b", "c", "d", "e"), 
method = c("1", "2", "3", "4"), ssfct = NULL, fctName, fctText)

FPL.4(p1, p2, fixed = c(NA, NA, NA, NA), names = c("b", "c", "d", "e"), ...)

Arguments

Details

The fractional polynomial dose-response models introduced by Namata et al. (2008) are implemented using the logistic model as base.

Value

The value returned is a list containing the nonlinear function, the self starter function and the parameter names.

Author(s)

Christian Ritz

References

Namata, Harriet and Aerts, Marc and Faes, Christel and Teunis, Peter (2008) Model Averaging in Microbial Risk Assessment Using Fractional Polynomials, Risk Analysis 28, 891–905.

See Also

Examples are found maED.

Description

Small plants of Galium aparine, growing in pots in a green house, were sprayed with the technical grade phenmidipham herbicide either alone or in mixture with an ester of oleic acid. The plants were allowed to grow in the green house for 14 days after herbicide treatment. Then the dry matter was measured per pot.

Usage

data(G.aparine)

Format

A data frame with 240 observations on the following 3 variables.

dose

a numeric vector of dose value (g/ha)

drymatter

a numeric vector of dry matter weights (mg/pot)

treatment

a numeric vector giving the grouping: 0: control, 1,2: herbicide formulations

Source

Cabanne, F., Gaudry, J. C. and Streibig, J. C. (1999) Influence of alkyl oleates on efficacy of phenmedipham applied as an acetone:water solution on Galium aparine, Weed Research, 39, 57–67.

Examples

## Fitting a model with a common control (so a single upper limit: "1")
G.aparine.m1 <- drm(drymatter ~ dose, treatment, data = G.aparine, 
pmodels = data.frame(treatment, treatment, 1, treatment), fct = LL.4())

## Visual inspection of fit
plot(G.aparine.m1, broken = TRUE)

## Lack of fit test
modelFit(G.aparine.m1)

## Summary output
summary(G.aparine.m1)

## Predicted values with se and confidence intervals
#predict(G.aparine.m1, interval = "confidence")
# long output

## Calculating the relative potency
EDcomp(G.aparine.m1, c(50,50))

## Showing the relative potency as a
## function of the response level
relpot(G.aparine.m1)
relpot(G.aparine.m1, interval = "delta")
# appears constant!

## Response level in percent
relpot(G.aparine.m1, scale = "percent")

## Fitting a reduced model (with a common slope parameter)
G.aparine.m2 <- drm(drymatter ~ dose, treatment, data = G.aparine, 
pmodels = data.frame(1, treatment, 1, treatment), fct = LL.4())

anova(G.aparine.m2, G.aparine.m1)

## Showing the relative potency 
relpot(G.aparine.m2)

## Fitting the same model in a different parameterisation
G.aparine.m3 <- drm(drymatter ~ dose, treatment, data = G.aparine, 
pmodels = data.frame(treatment, treatment, 1, treatment), fct = LL2.4())

EDcomp(G.aparine.m3, c(50, 50), logBase = exp(1))

Description

The gamma dose-response model is a four-parameter model derived from the cumulative distribution function of the gamma distribution.

Usage

gammadr(fixed = c(NA, NA, NA, NA), 
  names = c("b", "c", "d", "e"), fctName, fctText)

Arguments

Details

Following Wheeler and Bailer (2009) the model function is defined as follows:

$f(x) = c + (d-c) * pgamma(b*x, e, 1)$

This model is only suitable for increasing dose-response data.

Value

The value returned is a list containing the nonlinear function, the self starter function and the parameter names.

Author(s)

Christian Ritz

References

Wheeler, M. W., Bailer, A. J. (2009) Comparing model averaging with other model selection strategies for benchmark dose estimation, Environmental and Ecological Statistics, 16, 37–51.

Description

Model functions for fitting symmetric or skewed bell-shaped/biphasic dose-response patterns.

Usage

gaussian(fixed = c(NA, NA, NA, NA, NA), names = c("b", "c", "d", "e", "f"), 
    method = c("1", "2", "3", "4"), ssfct = NULL, fctName, fctText, loge = FALSE)
    
    lgaussian(fixed = c(NA, NA, NA, NA, NA), names = c("b", "c", 
    "d", "e", "f"), method = c("1", "2", "3", "4"), ssfct = NULL, 
    fctName, fctText, loge = FALSE)

Arguments

Details

Details yet to be provided.

Value

The value returned is a list containing the nonlinear function, the self starter function and the parameter names.

Note

The functions are for use with the function drm.

Author(s)

Christian Ritz

Description

Germination data were obtained from experiments involving the three species mungbean, rice, and wheat, which were opposed to different temperatures between 10 and 40 degrees Celsius. Experiments lasted at most 18 days.

Usage

data(germination)

Format

A data frame with 192 observations on the following 5 variables.

temp

a numeric vector of temperatures that seeds were exposed to

species

a factor with levels mungbean rice wheat

start

a numeric vector of left endpoints of the monitoring intervals

end

a numeric vector of right endpoints of the monitoring intervals

germinated

a numeric vector giving the numbers of seeds germinated

Details

For each of the three species mungbean, rice, and wheat, a total of 20 seeds were uniformly distributed on filter paper in a petri dish (diameter: 9.0cm) and then placed in dark climate cabinets with different temperatures (10, 16, 22, 28, 34, 40 degrees Celsius). Not all of the temperatures were applied to all species. The germinated seeds were counted and removed from the petri dish on a daily basis up to 18 days (or until all seeds had germinated). I

n this experiment we also assume that the upper limit of the proportion germinated is a parameter that has to be estimated from the data. Moreover, we assume that different combinations of species and temperature may lead to different germination curves with respect to slope, time required for 50% germination, and upper limit.

References

Ritz, C., Pipper, C. B. and Streibig, J. C. (2013) Analysis of germination data from agricultural experiments, Europ. J. Agronomy, 45, 1–6.

See Also

Analysis of a single germination curve is shown for chickweed.

Examples

## Fitting two-parameter log-logistic curves to each combination of species and temperature 
##  (upper limit fixed at 1)
## Note: Rows 24 and 62 are omitted from the dataset (all mungbean seeds germinated 
##  and thus no right-censoring in this case)

## germLL.2 <- drm(germinated ~ start + end, species:factor(temp), 
## data = germination[c(1:23, 25:61, 63:192), ], fct = LL.2(), type = "event")
## plot(germLL.2, ylim=c(0, 1.5), legendPos=c(2.5,1.5))  # plotting the fitted curves and the data
## summary(germLL.2)  # showing the parameter estimates

## Fitting two-parameter log-logistic curves to each combination of species and temperature
## Note: the argument "start" may be used for providing sensible initial 
##  parameter values for estimation procedure (is needed occasionally)
##  (initial values were obtained from the model fit germLL.2)
## Note also: the argument "upper" ensures that the upper limit cannot exceed 1
## (however, no restrictions are imposed on the two remaining parameters 
## (as indicated by an infinite value)

## germLL.3 <- drm(germinated~start+end, species:factor(temp), 
## data = germination[c(1:23, 25:61, 63:192), ], fct = LL.3(), type = "event",
## start = c(coef(germLL.2)[1:13], rep(0.7,13), coef(germLL.2)[14:26]), 
## upper = c(rep(Inf, 13), rep(1, 13), rep(Inf, 13)))

## Plotting the fitted curves and the data
## plot(germLL.3, ylim = c(0, 1.5), legendPos = c(2.5,1.5))

## Showing the parameter estimates
## summary(germLL.3)

## Showing the parameter estimates with robust standard errors
## library(lmtest)
## coeftest(germLL.3, vcov = sandwich) 

## Calculating t50 with associated standard errors
## ED(germLL.3, 50)

## Calculating t10, t20, t50 with 95% confidence intervals
## ED(germLL.3, c(10, 20, 50), interval = "delta")

## Comparing t50 between combinations by means of approximate t-tests
## compParm(germLL.3, "e", "-")

## Making plots of fitted regression curves for each species

## Plot for mungbean
#plot(germLL.3, log="", ylim=c(0, 1), xlim=c(0, 20), 
#level=c("mungbean:10", "mungbean:16"), 
#lty=2:3, lwd = 1.5,
#xlab="Time (days)", 
#ylab="Proportion germinated",  
#main="Mungbean",
#legendPos=c(3, 1.05), legendText=c(expression(10*degree), expression(16*degree)))

## Plot for rice
#plot(germLL.3, log="", ylim=c(0, 1), xlim=c(0, 20), 
#level=c("rice:16", "rice:22", "rice:28", "rice:34", "rice:40"), 
#lty=2:6, lwd = 1.5,
#xlab="Time (days)", 
#ylab="Proportion germinated",
#main="Rice",   
#pch=2:6,
#legendPos=c(3, 1.05), legendText=c(expression(16*degree), expression(22*degree), 
#expression(28*degree), expression(34*degree), expression(40*degree)))

## Plot for wheat
#plot(germLL.3, log="", ylim=c(0, 1), xlim=c(0, 20), 
#level=c("wheat:10", "wheat:16", "wheat:22", "wheat:28", "wheat:34", "wheat:40"), 
#lty=c("dashed","dotted","dotdash","longdash","twodash","232A"), lwd = 1.5,
#xlab="Time (days)", 
#ylab="Proportion germinated", 
#main="Wheat",
#legendPos=c(3, 1.05), 
#legendText=c(expression(10*degree), expression(16*degree), expression(22*degree), 
#expression(28*degree), expression(34*degree), expression(40*degree)))

Description

Function for showing the starting values of the model parameters used when fitting a dose-response model

Usage

getInitial(object)

Arguments

Value

A vector of starting values for the model parameters used to initialize the estimation procedure.

Note

This function is masking the standard function in the stats package.

Author(s)

Christian Ritz

Description

Display information about available, built-in dose-response models.

Usage

getMeanFunctions(noParm = NA, fname = NULL, flist = NULL, display =TRUE)

Arguments

Details

The arguments noParm and fname can be combined.

Value

An invisible list of functions or a list of strings with brief function descriptions is returned.

Author(s)

Christian Ritz

Examples

## Listing all functions
getMeanFunctions()

## Listing all functions with 4 parameters
getMeanFunctions(4)

## Listing all (log-)logistic functions
getMeanFunctions(fname = "L")

## Listing all three-parameter (log-)logistic or Weibull functions
getMeanFunctions(3, fname = c("LL", "W"))

## Listing all four-parameter (log-)logistic or Weibull functions
getMeanFunctions(4, fname = c("LL", "W"))

Description

The dataset has 7 mixtures, 8 dilutions, two replicates and 5 common control controls. Four observations are missing, giving a total of 113 observations.

Usage

data(glymet)

Format

A data frame with 113 observations on the following 3 variables.

dose

a numeric vector of dose values

pct

a numeric vector denoting the grouping according to the mixtures percentages

rgr

a numeric vector of response values (relative growth rates)

Details

The dataset is analysed in Soerensen et al (2007). The concentration addition model can be entertained for this dataset.

Source

The dataset is kindly provided by Nina Cedergreen, Department of Agricultural Sciences, Royal Veterinary and Agricultural University, Denmark.

References

Soerensen, H. and Cedergreen, N. and Skovgaard, I. M. and Streibig, J. C. (2007) An isobole-based statistical model and test for synergism/antagonism in binary mixture toxicity experiments, Environmental and Ecological Statistics, 14, 383–397.

Examples

## Fitting the model with freely varying ED50 values
glymet.free <- drm(rgr~dose, pct, data = glymet, 
fct = LL.3(), pmodels = list(~factor(pct) , ~1, ~factor(pct))) 

## Lack-of-fit test
modelFit(glymet.free)  # acceptable
summary(glymet.free)

## Plotting isobole structure
isobole(glymet.free, exchange=0.01)

## Fitting the concentration addition model
glymet.ca <- mixture(glymet.free, model = "CA")

## Comparing to model with freely varying e parameter
anova(glymet.ca, glymet.free)  # borderline accepted

## Plotting isobole based on concentration addition
isobole(glymet.free, glymet.ca, exchange = 0.01)  # acceptable fit

## Fitting the Hewlett model
glymet.hew <- mixture(glymet.free, model = "Hewlett")

### Comparing to model with freely varying e parameter
anova(glymet.ca, glymet.hew)  
# borderline accepted
# the Hewlett model offers no improvement over concentration addition

## Plotting isobole based on the Hewlett model
isobole(glymet.free, glymet.hew, exchange = 0.01)  
# no improvement over concentration addition

Description

This function provides a very general way of specifying the mean function of the decreasing or incresing Gompertz dose-response or growth curve models.

Usage

gompertz(fixed = c(NA, NA, NA, NA), names = c("b", "c", "d", "e"), 
  method = c("1", "2", "3", "4"), ssfct = NULL,
  fctName, fctText)

Arguments

Details

The Gompertz model is given by the mean function

$f(x) = c + (d-c)(\exp(-\exp(b(x-e))))$

and it is a dose-response/growth curve on the entire real axis, that is it is not limited to non-negative values even though this is the range for most dose-response and growth data. One consequence is that the curve needs not reach the lower asymptote at dose 0.

If

$b<0$

the mean function is increasing and it is decreasing for

$b>0$

. The decreasing Gompertz model is not a well-defined dose-response model and other dose-response models such as the Weibull models should be used instead.

Various re-parameterisations of the model are used in practice.

Value

The value returned is a list containing the non-linear function, the self starter function and the parameter names.

Note

The function is for use with the function drm, but typically the convenience functions G.2, G.3, G.3u, and G.4 should be used.

Author(s)

Christian Ritz

References

Seber, G. A. F. and Wild, C. J. (1989) Nonlinear Regression, New York: Wiley \& Sons (p. 331).

See Also

The Weibull model weibull2 is closely related to the Gompertz model.

Description

'gompertzd' provides a way of specifying the derivative of the Gompertz function as a dose-response model.

Usage

gompertzd(fixed = c(NA, NA), names = c("a", "b"))

Arguments

Details

The derivative of the Gompertz function is defined as

$f(x) = a \exp(bx-a/b(exp(bx)-1))$

For $a>0$ and $b$ not 0, the function is decreasing, equaling $a$ at $x=0$ and approaching 0 at plus infinity.

Value

The value returned is a list containing the model function, the self starter function and the parameter names.

Note

This function is for use with the function drm.

Author(s)

Christian Ritz

Description

For three days, moths of the tobacco budworm (Heliothis virescens) were exposed to doses of the pyrethroid trans-cypermethrin.

Usage

data(H.virescens)

Format

A data frame with 12 observations on the following 4 variables.

dose

a numeric vector of dose values ($\mu g$)

numdead

a numeric vector of dead or knocked-down moths

total

a numeric vector of total number of moths

sex

a factor with levels F M denoting a grouping according to sex

Details

In Venables and Ripley (2002), these data are analysed using a logistic regression with base-2 logarithm of dose as explanatory variable.

Source

Venables, W. N. and Ripley, B. D (2002) Modern Applied Statistics with S, New York: Springer (fourth edition).

Examples

## Fitting dose-response model (log-logistic with common slope)
Hv.m1 <- drm(numdead/total~dose, sex, weights = total, data = H.virescens, fct = LL.2(), 
pmodels = list(~ 1, ~ sex - 1), type = "binomial")
summary(Hv.m1)

## Fitting the same model as in Venables and Riply (2002)
Hv.m2 <- glm(cbind(numdead, total-numdead) ~ sex + I(log2(dose)) - 1, data = H.virescens, 
family = binomial)

## Comparing the fits
logLik(Hv.m1)
logLik(Hv.m2)

## Estimated ED values (matching those given in MASS)
ED(Hv.m1, c(25, 50, 75))

Description

Hat values (leverage values) and Cook's distance are provided for nonlinear dose-response model fits using the same formulas as in linear regression but based on the corresponding but approximate quantities available for nonlinear models.

Usage

## S3 method for class 'drc'
cooks.distance(model, ...)

  ## S3 method for class 'drc'
hatvalues(model, ...)

Arguments

Details

Hat values and Cook's distance are calculated using the formula given by Cook et al. (1986) and McCullagh and Nelder (1989).

The output values can be assessed in the same way as in linear regression.

Value

A vector of leverage values (hat values) or values of Cook's distance (one value per observation).

Author(s)

Christian Ritz

References

Cook, R. D. and Tsai, C.-L. and Wei, B. C. (1986) Bias in Nonlinear Regression, Biometrika 73, 615–623.

McCullagh, P. and Nelder, J. A. (1989) emphGeneralized Linear Models, Second edition, Chapman \& Hall/CRC.

Examples

ryegrass.LL.4 <- drm(rootl ~ conc, data = ryegrass, fct = LL.4())

hatvalues(ryegrass.LL.4)

cooks.distance(ryegrass.LL.4)

Description

The dataset contains measurements of mean arterial pressure (mmHG) and heart rate (b/min) for a baroreflex curve.

Usage

data(heartrate)

Format

A data frame with 18 observations on the following 2 variables.

pressure

a numeric vector containing measurements of arterial pressure.

rate

a numeric vector containing measurements of heart rate.

Details

The dataset is an example of an asymmetric dose-response curve, that is not easily handled using the log-logistic or Weibull models (LL.4, LL.5, W1.4 and W2.4), whereas the baro5 model provides a nice fit.

Source

Ricketts, J. H. and Head, G. A. (1999) A five-parameter logistic equation for investigating asymmetry of curvature in baroreflex studies, Am. J. Physiol. (Regulatory Integrative Comp. Physiol. 46), 277, 441–454.

Examples

## Fitting the baro5 model
heartrate.m1 <- drm(rate~pressure, data=heartrate, fct=baro5())
plot(heartrate.m1)

coef(heartrate.m1)

#Output:
#b1:(Intercept) b2:(Intercept)  c:(Intercept)  d:(Intercept)  e:(Intercept)
#      11.07984       46.67492      150.33588      351.29613       75.59392

## Inserting the estimated baro5 model function in deriv()
baro5Derivative <- deriv(~ 150.33588 + ((351.29613 - 150.33588)/
(1 + (1/(1 + exp((2 * 11.07984 * 46.67492/(11.07984 + 46.67492)) * 
(log(x) - log(75.59392 ))))) * (exp(11.07984 * (log(x) - log(75.59392)))) + 
(1 - (1/(1 + exp((2 * 11.07984 * 46.67492/(11.07984 + 46.67492)) * 
(log(x) - log(75.59392 )))))) * (exp(46.67492 * (log(x) - log(75.59392 )))))), "x", function(x){})

## Plotting the derivative
#pressureVector <- 50:100
pressureVector <- seq(50, 100, length.out=300)
derivativeVector <- attr(baro5Derivative(pressureVector), "gradient")
plot(pressureVector, derivativeVector, type = "l")

## Finding the minimum
pressureVector[which.min(derivativeVector)]

Description

'isobole' displays isobole based on EC/ED50 estimates from a log-logistic model. Additionally isoboles determined by the concentration addition model, Hewlett's model and Voelund's model can be added to the plot.

Usage

isobole(object1, object2, exchange = 1, cifactor = 2, ename = "e", 
xaxis = "100", xlab, ylab, xlim, ylim, ...)

Arguments

Details

The model fits to be supplied as first and optionally second argument are obtained using mixture and drm.

Value

No value is returned. Only used for the side effect: the isobologram shown.

Author(s)

Christian Ritz

References

Ritz, C. and Streibig, J. C. (2014) From additivity to synergism - A modelling perspective Synergy, 1, 22–29.

See Also

The examples in acidiq, glymet and mecter.

Description

In an experiment barley was grown in a hydroponic solution with a herbicide.

Usage

data(leaflength)

Format

A data frame with 42 observations on the following 2 variables.

Dose

a numeric vector

DW

a numeric vector

Details

The dataset exhibits a large hormetical effect.

Source

Nina Cedergreen, Royal Veterinary and Agricultural University, Denmark.

Examples

## Fitting a hormesis model
leaflength.crs4c1 <- drm(DW ~ Dose, data = leaflength, fct = CRS.4c())
plot(fitted(leaflength.crs4c1), residuals(leaflength.crs4c1))

leaflength.crs4c2 <- boxcox(drm(DW ~ Dose, data = leaflength, fct = CRS.4c()), 
method = "anova", plotit = FALSE)
summary(leaflength.crs4c2)

## Plottinf fitted curve and original data
plot(leaflength.crs4c2, broken = TRUE, conLevel = 0.001, type = "all", legend = FALSE, 
ylab = "Produced leaf length (cm)", xlab = "Metsulfuron-methyl (mg/l)",
main = "Hormesis: leaf length of barley")