Types of responses

We define a dose (metameter) as any pre-specified amount of biological, chemical, or radiation stress eliciting a certain, well-defined response. Other kinds of exposure or stress could also be imagined, e.g., time elapsed in germination experiments. The dose is a non-negative quantity and it is often but not always assumed to be measured without error as is often the case in designed experiments [8].

Specifically, we define the response to a given dose as the quantification of a biologically relevant effect and as such it is subject to random variation. The most common type is a continuous response such as biomass, enzyme activity, or optical density. A binary or aggregated binary (binomial) response is also frequently used to describe results such as dead/alive, immobile/mobile, or present/absent [9]. The response may also be discrete as in number of events observed in a specific time interval such as number of juveniles, offspring, or roots [10]; this setting also encompassses Wadley’s problem considered in Example 6 in S1 File [11]. Dose-response curves may also be used to summarize experiments where the response is the event time, that is the time elapsed before some specific event is being observed [12, 13].

Modelling the mean

The full specification of a statistical dose-response (regression) model involves specifying how the mean is described by a parametric function of dose as well as specifying assumptions about the distribution of the response.

We will focus on ways to model the mean trends through mostly s-shaped or related biphasic functions because these functions have in common that they reflect a basic understanding about the causal relationship between the dose and the response, e.g., when a dose increases the response monotonically decreases or increases one way or another towards minimum or maximum response limits, respectively. Consequently, these functions have turned out to be extremely versatile for describing various biological mechanisms involving model parameters that allow the interpretation of observed effects within a biologically plausible framework.

So we define dose-response models to be a collection of statistical models having a certain mean structure in common; this is not a strict mathematical definition, but rather a definition driven by applications. Consequently, dose-response models encompass a range of statistical models from nonlinear regression, generalized (non)linear regression, and parametric survival analysis.

Let y denote an observed response value, possibly aggregated in some way, corresponding to a dose value x ≥ 0. The values of y are often positive but may take arbitrary positive or negative values. Furthermore, we will assume that observation of y is subject to sampling variation, necessiating the specification of a statistical model describing the random variation. Specifically, we will focus on characterizing the mean of y (denoted E(y) below) in terms of a model function f that depends on the dose x: (1) So, for a given dose x, the corresponding observed response values will be distributed around f(x, β). The function f is completely known as it reflects the assumed relationship between x and y, except for the values of the model parameters β = (β₁, …, β_p), which will be estimated from the data to obtain the best fitting function. The remaining distributional assumptions on y will depend on the type of response. For instance, for a continuous response the normal distribution is commonly assumed whereas for a binary or quantal response the binomial distribution is commonly assumed.

Built-in dose-response models.

A large number of more or less well-known model functions are built-in in drc (see Table 1). These models are parameterized using a unified structure with a coefficient b denoting the steepness of the dose-response curve, c, d the lower and upper asymptotes or limits of the response, and, for some models, e the effective dose ED50.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. List of model functions and corresponding names of some of the most important built-in models available in drc.

https://doi.org/10.1371/journal.pone.0146021.t001

By far the log-logistic models are the most used dose-response models [9]. The four-parameter log-logistic model corresponds to the model function: (2) Two slightly different parameterizations are available: one where ED50 is a model parameter, that is e in Eq (2), and another where the logarithm of ED50 denoted by , say, is a model parameter as in the following Eq (3): (3) This second version of the four-parameter log-logistic model may be preferred for dose-response analysis involving very small datasets (<15–20) where a normally distributed parameter estimate may more reasonably be assumed on a logarithm-transformed dose scale than the original dose scale.

Model functions derived from the log-logistic model include the generalized five-parameter log-logistic model, which results in an asymmetric dose-response curve [14], the five-parameter baroflex model (with two slope parameters b) [15], and the non-monotonous hormesis models [16, 17]. Log-logistic type fractional polynomial dose-response models were introduced for a binomial response [18]. These models have been implemented in drc in novel, extended four-parameter versions that are also equally suitable for describing a dose-response curve for a continuous response. Biphastic functions obtained as the sum of two four-parameter log-logistic models may also be fitted using drc [19]. Recently, other types of biphasic dose-response models were proposed in the context of biosensors [20].

Log-normal models, which result in dose-responses curves very similar to curves obtained from the corresponding log-logistic models, and two types of asymmetric Weibull models are also available in drc. In passing we note that there is a close link between the two-parameter log-logistic, log-normal, and Weibull type I models available in drc (where only the two parameters b and e are not fixed) and the log-logistic, log-normal, and Weibull models available in the package survival [21]. Generalized four- and five-parameter versions of the gamma and (quadratic) multistage models, respectively, are also implemented [22]. Another built-in model is the so-called “no effect concentration” (NEC) threshold model [23].

Most of these functions are scale invariant in the sense that the magnitude of doses is accommodated by the model itself through the parameter b. For instance, for log-logistic functions the slope parameter b acts as a scaling factor, centering doses around 1. In contrast the Brain-Cousens and Cedergreen-Ritz-Streibig models are sensitive to the magnitudes of the doses, which may need to be manually up- or downscaled appropriately prior to model fitting. Moreover, some of these functions may be used for fitting both decreasing and increasing dose-response curves. One exception is the hormesis models, which can only describe decreasing trends, i.e., they are so-called inverse j-shaped hormesis models. However, modified functions for u-shaped hormesis models for increasing dose-response curves are also available [24]. It is also possible to define your own model function [25].

Nonlinear least squares

Assume that x₁, …, x_n are the dose values and y₁, …, y_n the corresponding observed response values. The nonlinear least squares estimates are obtained by minimizing the following sum of squares: (4) where w_i’s are user-specified weights (often left unspecified, i.e., equal to 1). In drm() weights are specified through the argument weights and they should be on the same scale as the response, e.g., expressed as standard deviations and not empirical variances. Weights may be used for addressing variance heterogeneity in the response. However, the transform-both-sides approach should be preferred over using often very imprecisely determined weights (see the next subsection).

Eq (4) has to be solved numerically in an iterative manner. One approach is to use iteratively weighted least squares as is done in nls(), which is part of standard installation in R [30]. Another approach is to use a general-purpose minimizer directly, as in drm() where optim() is used in combination with some pre-scaling of parameters. In our experience drm() is more robust than nls(); lack of robustness of nls() has also been pointed out previously [31].

The estimated variance-covariance of the parameter estimates () is obtained as the scaled inverse of the observed information matrix, which consists of second-order partial derivatives of f w.r.t. the model parameters β₁, …, β_p: (5) The scaling factor is the residual standard error, which is estimated in the same way as in linear regression based on residuals. The observed information matrix (“hessian”) in Eq (5) is approximated numerically in optim() upon convergence. However, it is not even a requirement that the inverse is available: if it is not then NAs are returned for standard errors; see Example 1 in S1 File.

Transform-both-sides approach

The transform-both-sides approach implies that both the response and the nonlinear model function are shifted by a constant C and then transformed by a suitable function g_λ [32]. The resulting model for the mean of the transformed response looks like this: (6) with the transformed response on the left-hand side and the transformed model function on the right-hand side. Usually the function g_λ is taken to be the Box-Cox transformation g_λ(y) = (y^λ − 1)/λ for some suitable choice of . The value λ = 1 implies no transformation whilst λ = 0 corresponds to the logarithm transformation. All other values of λ correspond to power transformations. It is noteworthy that the Box-Cox transformation may alleviate variance heterogeneity and some skewness in the distribution of the response and thus recover a normal distribution, but it may not remedy other problems with the distributional assumptions such as counts observed with ties [10].

The package drc provides a method boxcox, which has been implemented in much the same way as the corresponding method for linear models available in the package MASS [33]. There is, however, a choice between a profiling approach as used for linear models in MASS or a more robust analysis of variance (ANOVA) approach where the optimal λ is estimated from a more general ANOVA model, i.e., a linear model, and not from the specified dose-response model; the latter requires replicate observations for at least some doses.

Manual specification of the transformation is also possible through the arguments bcVal and bcAdd, which correspond to λ and C in Eq (6), respectively. Using a specific transformation (i.e., a specific value of λ and possibly C) based on previous experience should be preferred over using a data-driven choice.

Maximum likelihood estimation

In case of maximum likelihood estimation the following expression has to be maximized in terms of β values in order to obtain parameter estimates: (7) where l denotes the logarithm-transformed likelihood function for a single measurement or observation. The likelihood in turn depends on the distributional assumptions made for the dose-response model. The estimated variance-covariance of the parameter estimates is found as the inverse of the observed information matrix [34]: (8) Maximum likelihood estimation is used for fitting dose-response models to bionomial, count, and event-time data, including models for binomial data where natural immunity and/or natural mortality is present [11].

Robust estimation

One way to handle observations that seem to be extreme or deviating in one way or another as compared to the bulk of the observations is to apply a robust estimation procedure, which will weigh down the influence of such observations.

Robust nonlinear regression is available through the function nlrob() in the R package robustbase. This function relies on nls() through an iterated weighted least squares approach. However, no self starter functionality is currently available. The commercial software GraphPad also has some limited functionality for robust nonlinear regression, but the resulting model fits are provided without standard errors of parameter estimates [35].

In drc robust dose-response analysis with a continuous response is available using the same models as considered in the case of robust linear regression [33], except for the model based on Hampel’s ψ, which is currently not implemented. However, it may be difficult to fit such robust nonlinear models unless fairly accurate starting values for the model parameters are provided.

The estimated “information”-type variance-covariance has the following form [36, 37]: where the function ρ controls the influence of observations on the estimation procedure (e.g., ρ(y) = y² corresponds to ordinary nonlinear least squares estimation). Once parameter estimates have been obtained all methods and extractors available in drc may be used as if the model fit had been obtained using maximum likelihood estimation.

Constrained estimation

By default unconstrained estimation is carried out by drm() although for most models there are actually constraints on at least some of the model parameters, e.g., ED50 needs to be positive. Likewise, there may also be constraints needed to ensure uniqueness of the model fits, e.g., the slope parameter b should be negative for an increasing dose-response relationship to ensure the interpretation that c is the lower asymptote and d the upper asymptote. In practice such constraints are indirectly enforced through the appropriate choice of starting values.

However, in principle any of the aforementioned estimation procedures may be combined with the use of constraints where the range of one or more parameters is restricted, e.g., to certain intervals by setting lower and upper bounds (so-called box constraints) that are different from −∞ and ∞, respectively. This is possible with drm() using the arguments lowerl and upperl, which, if specified, will invoke constrained optimisation through the method “L-BFGS-B” in optim().

Self starter functions

The estimation procedures introduced in this section are iterative and to ensure convergence towards optimal parameter estimates they need to be initiated through the provision of starting values for the model parameters. The choice of starting values may crucially affect whether convergence is eventually achieved. Thus, availability of good starting values facilitate parameter estimation in nonlinear models. Starting values may be obtained either by using parameter estimates previously reported for similar experiments or, in a data-driven way, by using the dose-response data themselves to elicit relevant information. For instance, within weed science, routine calculation of starting values for log-logistic models from simple linear regression models of transformed responses dates back at least to [38], utilizing extending linearization techniques previously used for probit analysis [39].

The flexibility of R has made it possible to implement so-called self starter functions that return data-driven starting values for the model parameters. In R self starter functions are available for the function nls(). In drc the above-mentioned data-driven linearization technique has been extended to other types of dose-response models and it is now available for all built-in models. It is the default setting, which may, however, be overruled by supplying starting values manually through the argument start. The function getInitial() may be used to obtain the starting values that were actually used for obtaining a particular model fit.

More sophisticated self starter functions have been proposed [40], but so far they have not implemented in drc except for some log-logistics models.

Choice of model

Ideally, either a single dose-response model is fitted based on established knowledge or a number of plausible candidate models are identified and a model-averaging approach is adopted. However, in reality a decision-tree approach is typically used to find the best-fitting model among several candidates. For this approach the function mselect() may be useful as it returns values for information criteria or, if models are nested, p-values from lack-of-fit/goodness-of-fit tests as implemented in the function modelFit().

Inverse regression

Estimating doses corresponding to specific response levels is often of particular interest in dose-response analysis. Such estimates are obtained by solving an inverse regression problem, which may be approached either through after-fitting or re-parameterization. After-fitting means solving the inverse regression problem after having fitted the model. For some models the solution may be derived explicitly (e.g., log-logistic models [2, 44]) whereas for other models only numerical solutions are available (e.g., hormesis models [17]). After-fitting relies on the so-called delta method, which enables the calculation of approximate standard errors for derived parameter estimates [45]. Reparameterization means fitting the model using a reparameterized model function where the solution is an actual model parameter.

A practical implication of the after-fitting approach is that it suffices to fit a dose-response model once, in a parameterization that has proven to be the most stable for estimation. In our opinion, after-fitting provides a major improvement over the often quite instable re-parameterization approaches still being used with other statistical software programmes [46].

Prediction

There are a number of functions available in drc for obtaining predicted values from a dose-response model fit. Notably, there are the general-purpose methods fitted() and predict(). The latter may also provide confidence or prediction intervals (controlled through the argument type). The function backfit() may be used for obtaining the estimated doses that correspond to the observed average response at each dose.

The predict() method is useful for obtaining data for construction of plots using the traditional graphical functionality in R. In particular, the plot() method may be used to show the original or summarized data together with the fitted dose-response curve(s) superimposed. Such plots are often shown in articles as they are useful for assessing model fit and they could be used together with residual and normal probability plots. By default plot() uses a logarithmic dose axis, which may be switched off with the argument log = “”. The argument broken allows dose 0 also to be shown in plot together with the logarithmic axis. Normalization of the response relative to the control (dose 0) may be achieved using normal = TRUE and possibly the argument normref for setting the reference, which is by default 1, as suggested by [52]. This normalization is based on the model fit and it does not involve modification of the original data.

For more sophisticated plots, however, we suggest to use the package ggplot2 as shown in Example 2 in S1 File.