drda: An R package for dose-response data analysis

Analysis of dose-response data is an important step in many scientific disciplines, including but not limited to pharmacology, toxicology, and epidemiology. The R package drda is designed to facilitate the analysis of dose-response data by implementing efficient and accurate functions with a familiar interface. With drda, it is possible to fit models by the method of least squares, perform goodness of fit tests, and conduct model selection. Compared to other similar packages, drda provides, in general, more accurate estimates in the least-squares sense. This result is achieved by a smart choice of the starting point in the optimization algorithm and by implementing the Newton method with a trust region with analytical gradients and Hessian matrices. In this article, drda is presented through the description of its methodological components and examples of its user-friendly functions. Performance is finally evaluated using a real, large-scale drug sensitivity screening dataset.


Introduction
Inferring dose-response relationships is indispensable in many scientific disciplines.In cancer research, for example, estimating the magnitude of a chemical compound effect on cancer cells holds substantial promise for clinical applications.The dose-response relationship can be modeled via a nonlinear parametric function expressed as a dose-response curve.The choice of the curve is achieved by selecting the parameter values that lead to the best fit of the observations.Since conclusions about compound efficacy, for example its potency to kill the cells, are based on the estimated dose-response curve, it is of great importance to determine the curve parameters as accurately as possible.
Currently, there are multiple R (R Core Team 2022) packages that provide tools for dose-smart initialization procedure to increase the chances of converging to the global solution; (ii) applying a more advanced Newton method with a trust region and (iii) relying on analytical gradient and Hessian formulas instead of numerical approximations in problematic cases to assure proper convergence; (iv) computing confidence intervals for the whole dose-response curve using multiple comparisons test correction.Besides, similar to other packages, drda provides tools to compare the fitted curve against a linear model or other logistic models, to compute confidence intervals for the estimated parameters, and to plot multiple models in a user-friendly manner.
The most important feature of any optimization routine remains the closeness of its solution to the true optimizer of the objective function, such as the maximum likelihood estimate.One of the main disadvantages when it comes to numerical optimization is the possibility of converging to a local optimum instead of the correct answer we seek.This situation can easily happen when either the function is not well approximated by a quadratic shape in a neighborhood of the current candidate solution, or when the starting point is far from the global optimum (either the algorithm is not able to converge in a reasonable number of steps or it simply converges to a wrong solution).To cope with such scenarios, we implement here the Newton method with a trust region (Steihaug 1983), which has been shown to be a robust optimization technique for mitigating issues usually encountered in optimization problems.The method is more stable than other Newton-based methods, especially for cases when it is problematic to approximate a function with a quadratic surface (Sorensen 1982).Additionally, drda uses a multi-step initialization algorithm in order to ensure the right direction in the optimization routine.With our strategy, drda is able to find the true least squares estimate in problematic cases where the drc, nplr, and DoseFinding packages instead fail.
Once the least-squares estimate is found, drda provides the user with routines for assessing goodness-of-fit and reliability of the estimates.Assuming a Gaussian distribution with equal variance for the observed data, it is possible to compare the fitted model against, for example, a flat horizontal line or a logistic model with a different number of parameters.The drda package provides the likelihood-ratio test (LRT), the Akaike information criterion (AIC, Akaike 1974), and the Bayesian information criterion (BIC, Schwarz 1978) as a way to compare the goodness-of-fit of competing models.
The paper is organized as follows: We first describe the methodological components of drda in Section 2; show how the package is implemented in Section 3 ( including practical examples in Section 3.2); and provide a comparison of drda against packages drc, nplr, and DoseFinding using simulations and a high-throughput dose-response dataset in Section 4. We conclude the article with a summary and discussion in Section 5.

Generalized logistic and log-logistic function
Package drda implements the generalized logistic and log-logistic functions as the core models for fitting dose-response data.The generalized logistic function, also known as Richards' curve (Richards 1959), is the 5-parameter function solution to the differential equation where ψ = (α, δ, η, ϕ, ν) ⊤ and x ∈ R.
Our constraints have the benefit of giving the five parameters a clear and easy interpretation: α is the value of the function when x approaches negative infinity, δ is the (signed) height of the curve, that is α + δ is the value of the function when x approaches positive infinity, η is the steepness or growth rate of the curve, ϕ is related to the mid-value of the function, and ν regulates at which asymptote is the curve maximum growth.In dose-response data analysis, parameter δ represents the maximum effect attributable to the drug and it is often denoted in the literature as E max (Macdougall 2006).Refer to Figure 1 for a visual explanation of the five parameters.
When ν = 1 we obtain the 4-parameter logistic function.If we also set (α, δ) ⊤ = (0, 1) ⊤ or (α, δ) ⊤ = (1, −1) ⊤ we obtain the 2-parameter logistic function.When ν = 1 the parameter ϕ represents the value at which the function is equal to its midpoint, that is α + δ/2.In such a case, as a measure of drug potency, ϕ is also known as the half maximal effective logconcentration or log-EC 50 .As a measure of antagonist drug potency, ϕ is also known as the half maximal inhibitory log-concentration (log-IC 50 ).When ν → 0 we obtain the Gompertz function, i.e., lim The domain of the logistic function is the whole real line, that is x ∈ R. In dose-response data analysis, this means that the logistic function is a model for mapping log-doses to responses.When data is provided on the actual dose-scale, we can use the log-logistic function instead: where d = e x and ϕ + = e ϕ .Interpretation of parameters is equivalent to the logistic function with the only exception of ϕ + being now defined on the positive real line.Note that the log-logistic function is well-defined also for a dose of 0, where the function is equal to α.
In the literature, the log-logistic function with ν = 1 is also referred to as the E max model (Macdougall 2006).

Normal nonlinear regression
For a particular dose d k (k = 1, . . ., m) let (y ki , w ki ) ⊤ represent respectively the ith observed outcome (i = 1, . . ., n k ) and its associated positive weight.If observations have all the same importance, we simply set w ki = 1 for all k and i.We assume that each unit has expected value and variance where µ(d k ; ψ) is a nonlinear function of the dose d k and a vector of unknown parameters ψ.Parameter σ > 0 is instead the standard deviation common to all observations.In our package, µ(d k ; ψ) is simply the generalized log-logistic function (Equation 2) or the generalized logistic function (Equation 1) with the transformation By assuming the observations to be stochastically independent and normally distributed, the joint log-likelihood function is where n k is the sample size at dose k, n = k n k is the total sample size, ȳk = ( i w ki y ki )/w k.
is the weighted average corresponding to dose d k and w k. = i w ki .The maximum likelihood estimate ψ is obtained by minimizing the residual sum of squares from the means, i.e., The maximum likelihood estimate of the variance is while its unbiased estimate is where p is the total number of parameters estimated from the data.
For convenience we will use from now on the simplified notation µ k to denote the function µ(d k ; ψ).It is important to remember that µ k will always be a function of a dose d k and a particular parameter value ψ.We will also use the notation g (s) and g (st) to denote respectively the first-and second-order partial derivatives of function g(ψ), with respect first to ψ s and then ψ t .
Partial derivatives of the sum of squares g(ψ) are The gradient and Hessian of g(ψ) are therefore From the previous expressions we can easily retrieve the observed Fisher information matrix, that is the negative Hessian matrix of the log-likelihood evaluated at the maximum likelihood estimate, as where It is also worth noting that the (expected) Fisher information matrix is (5)

Statistical inference
When closed-form solutions of maximum likelihood estimates are missing, also closed-form expressions of other inferential quantities are not available.Fortunately, we can still rely on asymptotic, large sample size considerations, to obtain approximate values of quantities of interest.Obviously, the larger the sample size the better the approximation.
Using either versions (Equation 4or Equation 5) of the Fisher information matrix we can calculate approximate confidence intervals.In fact, we can think of the Fisher information matrix as an approximate precision matrix, so that we only have to invert the matrix and take diagonal elements as approximate variance estimates.In our package we use the observed Fisher information matrix (Equation 4) because it is shown to perform better with finite sample sizes (Efron and Hinkley 1978).As an example, an approximate confidence interval for generic parameter where t n−p,α is the appropriate quantile of level α of a Student's t-distribution with n − p degrees of freedom and I( ψ, σ) −1 jj is the jth element in the diagonal of the inverse observed Fisher information matrix.Using the delta method (Oehlert 1992) we can compute approximate point-wise confidence intervals for the mean function or for a new, yet to be observed, value y(d) We can also construct a (conservative and approximate) confidence band over the whole mean function µ(•; ψ) with the correction proposed by Gsteiger, Bretz, and Liu (2011): where q p,α is the appropriate quantile of level α of a χ 2 -distribution with p degrees of freedom.

Optimization by the Newton method with a trust region
A closed-form formula of the maximum likelihood estimate ψ, that is the solution of Equation 3, is in general not available for nonlinear regression models.We can, however, try to minimize numerically the sum of squares g(ψ).Suppose that our algorithm is at iteration t with current solution ψ t .We want to find a new step u such that g(ψ t + u) < g(ψ t ).We start by illustrating the standard Newton method.We approximate our function by a second-order Taylor expansion, that is The theoretical minimum is attained when the gradient with respect to u is zero, that is where 0 < γ ≤ 1 is a modifier of the step size for ensuring convergence (Armijo 1966).
When the method converges the algorithm is quadratically fast, or at least superlinear (Bonnans, Gilbert, Lemarechal, and Sagastizábal 2006): the closer g(ψ) is to a quadratic function the better is its Taylor approximation and the better are the algorithm convergence properties.
When the Hessian matrix is almost singular, it is still possible to apply quasi-Newton methods (Luenberger and Ye 2008) to (try to) avoid convergence problems.In our nonlinear regression setting we might have the extra complication of an objective function far from a quadratic shape, so that the (quasi-)Newton method might fail to converge.Although these situations can be thought to be rare, they are often encountered in real applications.For example, in Figure 2 we show a problematic surface that the quasi-Newton BFGS algorithm, as implemented by the base R function optim(), is not able to properly explore.
We will try to overcome the issues in the optimization by focusing our search only in a neighborhood of the current estimate, that is using a trust-region around the current solution ψ t .The problem to solve is now where ∆ t > 0 is the trust-region radius.Our implementation is based on the exposition of Nocedal and Wright (2006) and follows closely that of Mogensen and Riseth (2018).Briefly, at each iteration we compute the standard Newton's step and accept the new solution if it is within the trust-region.If the Newton's step is outside the admissible region we try an alternative step by a linear combination of the Newton's step and the steepest descent step, with the constraint that its length is exactly equal to the radius ∆ t (dogleg method).This new alternative step is then accepted or rejected on the basis of the actual reduction in the function value.The region radius ∆ t+1 for iteration t + 1 is adjusted according to the length and acceptance of the step just computed.For more details, we refer the reader to the extensive discussion found in Nocedal and Wright (2006).

Algorithm initialization
One of the major challenges in fitting nonlinear regression models is choosing a good starting point for initializing the optimization algorithm.Looking at the example in Figure 2, the choice of ψ 0 = (1, −1, 1, 0) ⊤ made the BFGS algorithm converge to a local optimum while a global optimum might have been found if a better starting point was instead chosen.
First of all, we present the closed-form maximum likelihood estimates α and δ when all other parameters have been fixed.For the logistic function define For the log-logistic function define h k = (d η /(d η + νϕ η + )) 1/ν .Assume h k to be known.Our mean function is now while the residual sum of squares becomes with gradient It is easy to prove that the gradient is equal to zero when Our initialization strategy is similar for both the logistic and log-logistic functions and can be summarized in the following steps.Set ν = 1 by default to focus on the 4-parameter version.Define a = min{ȳ}−ϵ and b = max{ȳ}+ϵ, where ϵ is a very small number to avoid numerical problems.By construction, is defined in the range (0, 1).When data is well-behaved and doses span properly the range of the function we can consider a and b to be approximations of α and α + δ respectively, that is z k ≈ h k .We now have the relationship Remember that ϕ = log(ϕ + ) and x k = log(d k ), thereby to avoid the logarithm of zero in the log-logistic function we use in our implementation . By fitting the simple linear regression u = β ⊤ x we obtain estimates β0 and β1 .Our initial estimates are therefore η 0 = | β1 | and ϕ 0 = − β0 / β1 .Note that for a monotonically decreasing function (δ < 0) the coefficient β1 is negative, therefore to enforce the constraint η > 0 we use the absolute value of β1 .Finally, we set α 0 and δ 0 at their maximum likelihood estimates by plugging η 0 and ϕ 0 into Equation 6.
In general, parameter ψ 0 = (α 0 , δ 0 , η 0 , ϕ 0 , 1) ⊤ is a good starting point for a numerical optimization algorithm.For data that does not behave well we might obtain a very bad starting point because of a poor approximation z k ≈ h k .To avoid problems with corner cases we also perform a grid search over the parameter space of (η, ϕ, ν) ⊤ using a space-filling maximum entropy design with 250 samples (Santner, Williams, and Notz 2018;Dupuy, Helbert, and Franco 2015).For each tested parameter (η, ϕ, ν) ⊤ in the grid we compute the corresponding values of α and δ using Equation 6.The parameter vector corresponding to the smallest residual sum of squares is chosen as our second candidate initial point.To further increase our chances of converging to the global optimum we add to the pool of initial points also two sub-optimal parameters (in the current implementation we use the parameters associated with the fifth and eighth lowest residual sum of squares).
The four chosen parameter vectors are passed in turn as starting points to the nlminb() function.nlminb() is a very fast and efficient function and for well-behaved data it converges in general to the global optimum.In this case, the best out of the four solutions is also the solution to our optimization problem.To keep refining the solution in cases of problematic data we feed the current solution as a starting point to our own trust region implementation based on analytical gradient and Hessian.

General overview
The main function of drda is drda() with signature drda(formula, data, subset, weights, na.action, mean_function = "logistic4", lower_bound = NULL, upper_bound = NULL, start = NULL, max_iter = 1000) The first argument, formula, is a symbolic representation in the form y ~x of the model to be fitted, where y is the vector of responses and x is the vector of log-doses.
data is an optional argument, typically a 'data.frame'object, containing the variables in the model.When data is not specified the variables are taken from the environment where the function is being called.subset is a logical vector, or a vector of indices, specifying the portion of data to be used for model fitting.
weights is an optional argument that specifies the weights to be used for fitting.Usually weights are used in situations where observations are not equally informative, i.e., when it is known that some of the observations should have a smaller or larger impact on the fitting process.If the weights argument is not provided then the ordinary least squares method is applied.na.action defines a function for handling NAs found in data.The default option is to use na.omit(), i.e., to remove all data points associated with the missing values.
Arguments lower_bound and upper_bound are used for performing constrained optimization.They serve as the minimum and maximum values allowed for the model parameters.
They are vectors of length equal to the number of parameters of the model specified by the mean_function argument.Values -Inf and Inf are allowed.The parameters for the 5parameter (log-)logistic function are listed in the following order: α, δ, η, ϕ, ν.For the other models the order is preserved but some of the parameters are excluded.Obviously, values in upper_bound must be greater than or equal to the corresponding values in lower_bound.
start represents an optional vector of starting values for the parameters.
Finally, the max_iter argument sets the value for the maximum number of iterations in the optimization algorithm.
To evaluate the efficacy of the treatment it is also possible to compute the normalized area under or above the curve: The two-element vector xlim defines the interval of integration with respect to x.The twoelement vector ylim defines the theoretical minimum and maximum values for the response variable y.Therefore, xlim and ylim together define a rectangle that is partitioned into two regions by the dose-response curve.The normalized area under the curve (NAUC) is defined as the area of the "lower" rectangle region divided by the total area of the rectangle.The normalized area above the curve (NAAC) is simply its complement, i.e., 1− NAUC.Default value of xlim is c(-10, 10) for the logistic function and c(0, 1000) for the loglogistic function.The xlim default values were chosen on the basis of dose ranges that are commonly found in the literature.In the majority of real applications the response variable y is a relative measure against a control treatment, therefore the default value for ylim is chosen to be c(0, 1).
It is also possible to estimate effective doses at particular response levels: effective_dose(drda_object, y, type, level) The effective_dose() function estimates the doses at which the fitted curve is equal to the specified response values y.Response values can either be passed on a relative scale (type = "relative"), in which case the value of y is expected to be in the range (0, 1), or on an absolute scale (type = "absolute"), in which case the value of y can span the whole real line.Additionally, it is possible to specify the confidence interval level via the level parameter.Default values are y = 0.5, type = "relative", and level = 0.95.

Usage examples
To demonstrate how to use drda, we create the following example data: This example imitates an experiment where seven drug doses have been tested three times each.Relative viability measures have been obtained for each dose-replicate pair and, in this case, comprise 21 values in the (0, 1) interval.Note that any finite real number is accepted as a possible valid outcome.

Default fitting
After loading the package, the drda() function can be applied directly to the two variables with the mean_function selected to be one of the log-logistic functions.

R> summary(fit_l4)
Call: drda(formula = y ~x, data = test_data, mean_function = "l4") The summary() function provides information about the Pearson residuals, parameters' and residual standard error estimates, and their 95% confidence intervals.Together with the actual point estimates, the widths of confidence intervals are a good starting point for assessing the reliability of the model fit.The values of the log-likelihood function, AIC, and BIC are also provided.Finally, the summary() function warns the user if the algorithm converged and if so, in how many iterations.
Parameter estimates can be accessed using the coef() and sigma() functions.

Model comparison and selection
The anova() function can be used to compare competing models within the same logistic family of models.The constant model (δ = 0), i.e., a flat horizontal line, is always included by default in the comparisons.When the model being fitted is not the 5-parameter function, the latter is always included as the general reference model in the likelihood-ratio test.
R> fit_l2 <-drda(y ~x, data = test_data, mean_function = "logistic2") R> anova(fit_l2) Note that the p value refers here to the likelihood-ratio test with a χ 2 -distribution asymptotic approximation.In this particular case we are testing the null hypothesis that the complete 5parameter logistic function is equivalent, likelihood-wise, to our 2-parameter logistic function.The significant result indicates that the 2-parameter logistic function provides a worse fit for the observed data compared to a 5-parameter logistic function.

Analysis of Deviance
R> fit_gz <-drda(y ~x, data = test_data, mean_function = "gompertz") R> fit_l4 <-drda(y ~x, data = test_data, mean_function = "logistic4") R> anova(fit_l2, fit_gz, fit_l4) These results indicate the 4-parameter logistic function as the best fit for the data.Not only has the model the lowest AIC value, but the likelihood-ratio test (Model 5 vs. Model 4) is also not significant.Indeed, the data was generated from a 4-parameter logistic function with ψ = (0.86, −0.84, 1, −2) ⊤ and σ = 0.05.Note that the likelihood-ratio test between the 4-parameter logistic model and the Gompertz model could not be performed because of an equal number of parameters (they are not nested models).

Weighted fitting
In case when not all of the observations should be utilized equally in the model, the weights argument can be provided to the drda() function.All the generic functions described above also apply to a weighted fit object.

Constrained optimization
The drda() function allows the choice of admissible values for the parameters by setting the lower_bound and upper_bound arguments appropriately.Unconstrained parameters are set to -Inf and Inf respectively.While setting the constraints manually, one should be careful in choosing the values as the optimization problem might become very difficult to solve within a reasonable number of iterations.
In the next example α is fixed to 1, δ is fixed to −1, the growth rate η is allowed to vary in [0, 5], while the midpoint parameter ϕ is left unconstrained.

R> naac(fit_l4)
[1] 0.622 To allow the values to be comparable between different compounds and/or studies, the function sets a hard constraint on both the x and y variables (see Section 3.1).However, the intervals can be easily changed if needed.
R> effective_dose(fit_l4, y = c(0.The missing values are a consequence of the fact that the estimated upper bound of the curve is α = 0.88 and there is no dose such that µ(log(d); ψ) = 0.95.

Benchmarking
In this section we evaluate the performance and accuracy of drda.We want to compare its performance against the three packages described in Section 1.To do that we select the most widely used model for dose-response curve fitting, that is the 4-parameter log-logistic function.It is the only model implemented in all the four packages we aim to compare.As a control case, we also add to the comparison the nls() function from package stats.With nls() the doses are first log-transformed and the fit is done using the SSfpl() self-starting function.
Define, for each parameter, the following possible values: All the possible combinations of the previous values form a set of 162 parameter vectors.
For each parameter vector we simulate 100 datasets, with 7 doses and 3 replicates per dose, from a normal distribution centered on the corresponding log-logistic function.The seven doses are uniformly distributed on the log10 scale: 0.0001, 0.001, 0.01, 0.1, 1, 10, 100.Finally, we apply the fitting function of each package on each one of the 16200 simulated datasets.
For drm() from package drc we select the 4-parameter log-logistic model with argument fct = LL.4() and fix the maximum number of iterations to 10000, similar to drda().For nplr() from package nplr we set LPweight to 0 for the ordinary least squares method.We fix npars to four for the 4-parameter log-logistic model.For fitMod() from package DoseFinding we choose the 4-parameter log-logistic model by setting model = "sigEmax" (see Section 2.1).Since the fitMod() function requires the user to set constraints on the nonlinear parameters, we use the default values bnds = defBnds(max(dose))$sigEmax.For each dataset we store the residual sum of squares (RSS), convergence status, and the elapsed time of running the function.
Since all packages are solving the same optimization problem, i.e., minimization of the residual sum of squares, we define the absolute relative error of package k as For practical purposes small absolute relative errors can be considered equivalent to zero.
According to our simulation results (Table 1), drda provides the most accurate estimates: in its worst performance drda achieved a relative error of just 0.2% from the best solution.
Overall, drda converged to a solution 100% of the times, drc 83.8% of the times, while function nls() did so only 25.9% of the times.We cannot provide an accurate convergence rate for DoseFinding and nplr because the information is missing from their returned fit object.
The higher accuracy comes at a small computational cost as more steps are needed for exploring the parameter space.Our data analysis reveals that fitMOD() and nplr() are the fastest functions to complete the fit (mean elapsed time of 0.008s and 0.012s respectively).drda found the global optimum in 0.514 seconds on average.For completeness, drm() and nls() needed an average of 0.04 and 0.013 seconds respectively.
We further tested drda on a real large-scale drug sensitivity dataset downloaded from the Cancer Therapeutics Response Portal (CTRP, Rees et al. 2016;Seashore-Ludlow et al. 2015;Basu et al. 2013)  viability measures span the (0.0019, 2.864) interval.To speed up the benchmark process we sampled at random 15000 datasets from the full set.
Similarly to our simulation study, we fitted the 4-parameter log-logistic model with each one of the four packages.For each cell line-drug-package triple we performed the benchmark and summarized the results in Table 2.
drda is flagged as the absolute best fit in 96.97% of the cases.When we only consider relative errors greater than 1%, the percentage raises to 99.53% (70.87% for DoseFinding, 64.95% for drc, 76.14% for nplr).When compared directly against the other packages, drda outperforms DoseFinding in 29.07% of the cases (worse for 0.43%), drc in 34.81% of cases (worse for 0.02%), and nplr in 23.58% of the cases (worse for 0.06%).When nls() successfully converged, drda was better in 4.9% and worse in 0.02% of the cases.
Similar to the simulation study results, fitMod() and nplr() are again the fastest functions to complete the fit (mean elapsed time of 0.009s and 0.013s respectively).On average drda found the global optimum (or a very close solution) in 0.3 seconds, drm() in 0.092 seconds, while nls() in 0.019 seconds.
Results shown in this section were obtained with R version 4.2.1, running under Rocky Linux 8.6 with an Intel Xeon Gold 6230 x86_64 Processor, using the Intel oneAPI Math Kernel Library 2022.1 as the BLAS/LAPACK back-end.Because of system-specific numerical instabilities, executing the provided replication materials within a different environment will likely produce minor differences from what is shown here.

Summary and discussion
In this paper we introduced the drda package, aimed at evaluating dose-response relationship to advance our understanding of biological processes or pharmacological safety.These types of experiments are of high importance in drug discovery, as they establish an essential step for subsequent therapeutic advances.An appropriate interpretation of the experimental data is grounded on a reliable estimation of the dose-response relationship.Therefore, it is imperative to provide advanced optimization methods that allow more accurate estimation of dose-response parameters, and the assessment of their statistical significance.
One of the main limitations of most optimization procedures is their convergence to local solutions.The basic quasi-Newton methods applied to logistic curve fitting are sensitive to the selection of a starting point and to cases where data is non-informative.Our package effectively overcomes the convergence problem as we implement a comprehensive multi-step initialization procedure.In case of problematic data, we additionally rely on our own trust region Newton method implementation based on analytical gradient and Hessian.In addition to standard routines, the package allows a user to evaluate a model fitness via the assessment of confidence intervals for the whole dose-response curve, estimation of effective doses and advanced plot options.
We have compared our package with the three state-of-the-art packages -DoseFinding, drc, and nplr.Using simulations and a large-scale drug screening study, we have shown that drda has clearly outperformed the other three packages in terms of accuracy.Despite the fact that our package is on average slower than the other three packages, its gain in accuracy is a favorable compromise.For most, if not all, experimental applications, accuracy has a higher priority.
The current version of drda provides optimization tools for continuous data and relies on the family of logistic functions only.In most applications the dose-response relationship is fitted with a logistic function.However, in some specific scenarios, other models can be a better description of the dose-response relationship.For example, a linear no-threshold model might be preferred for fitting dose-response data in radiation protection studies.We plan to extend our package to cover such cases in the future.The package is currently completely implemented in base R, therefore there are still many opportunities for improving its performance, by, for example, refactoring core critical functions in C.
If a researcher is looking for a package providing improved accuracy at a relatively low speed cost, drda might provide a viable option.
Package drda is available on the Comprehensive R Archive Network (CRAN) at https: //CRAN.R-project.org/package=drda.Users are encouraged to contribute to package development at https://github.com/albertopessia/drda.

Figure 3 :Figure 4 :Figure 5 :
Figure 3: Basic plot of the data and the maximum likelihood curve for the logistic function with 5 parameters.
Compare the estimated values with those of the log-logistic fit:

Table Model
Full)Model 4 is the best model according to the Akaike Information Criterion.
75, 0.95))By default effective_dose() uses a relative scale for the response, that is the previous results are for the 75% and 95% response.Compare them to the actual values of 0.75 and 0.95:

Table 1 :
Summary statistics of the RSS relative error for simulated models.A total of 162 parameter vectors were analyzed.For each parameter, 100 datasets with 7 doses and 3 replicates per dose were simulated.

Table 2 :
Summary statistics of the RSS relative error as measured on the CTRPv2 data.