onls: Weighted Orthogonal Nonlinear Least-Squares Regression...

In onls: Orthogonal Nonlinear Least-Squares Regression

Weighted Orthogonal Nonlinear Least-Squares Regression (Orthogonal Distance Regression)

Description

Fits nonlinear orthogonal-distance regression (ODR, also called errors-in-variables) models in the manner of ODRPACK (Boggs, Byrd and Schnabel, 1987; Boggs et al., 1992): the model parameters and one foot-point correction per observation are estimated simultaneously by minimizing the same weighted sum of squared response and predictor residuals that ODRPACK minimizes, using a Levenberg-Marquardt algorithm on the joint problem. In contrast to ordinary nonlinear least squares, onls() allows measurement error in both predictors and response variables. Predictor uncertainty may be specified through predictor-specific standard deviations, observation-specific error structures, or full covariance matrices allowing correlated predictor errors. Response uncertainty may be supplied globally or on a per-observation basis and can be combined with observation weights. The algorithm supports single- and multi-predictor nonlinear models, parameter bounds, fixed parameters, and covariance estimation for the fitted parameters that is equivalent to the ODRPACK (Gauss-Newton) covariance.

Usage

  onls(formula, data, start = NULL, weights = NULL, sigma_x = NULL, sigma_y = 1, 
       known_sigma = NULL, extend = NULL, window = NULL, control = list(), 
       lower = NULL, upper = NULL, fixed = NULL, subset = NULL, na.action = NULL, 
       trace = FALSE, algorithm = c("auto", "structured", "dense"))

Arguments

formula	A two-sided nonlinear model formula describing the response as a nonlinear function of one or more predictors and unknown model parameters. The response must be a plain variable name (no transformation of the left-hand side).
data	A data frame containing the variables appearing in formula.
start	A named list (or named vector) of starting values for all model parameters.
weights	Optional non-negative observation weights. As in nls, weights, subset, sigma_x and sigma_y are evaluated in data first and then in the calling environment, so that columns of data can be named directly (e.g. weights = w, subset = conc > 0.1; NA in a logical subset drops the observation). The weights are incorporated into the orthogonal-distance criterion through the response precision matrix and therefore affect both the foot-point corrections and the parameter estimates. When supplied, weights must have length equal to the number of observations.
sigma_x	Specification of predictor measurement error, as a standard deviation (not a variance or a weight). Accepts NULL, a scalar, a length-p vector, a p \times p covariance matrix, or an n \times p matrix of observation-specific standard deviations (a length-n vector is accepted as shorthand for the latter when p = 1). See 'Details'.
sigma_y	Response measurement error, as a standard deviation. May be either a single positive value applied to all observations or a length-n vector of observation-specific standard deviations.
known_sigma	Logical indicating whether the supplied (or default) sigma_x/sigma_y should be treated as known measurement-error standard deviations rather than as estimated from the data. If NULL (default), it is set to TRUE if sigma_x is supplied or sigma_y is explicitly passed by the user (regardless of its value), and FALSE otherwise. Standard errors are scaled by the reduced chi-square statistic (see reduced_chisq in 'Value') only when known_sigma = FALSE; when known_sigma = TRUE the supplied variances are assumed to already be on the correct absolute scale and no such scaling is applied. If known_sigma is TRUE but sigma_x is not given, a warning points out that the predictor errors are then taken to have standard deviation 1 in the units of each predictor.
extend	Optional numeric vector of length one or two. If NULL (default), the foot-point corrections are unbounded, as in ODRPACK. If supplied (single-predictor models only), the foot points are confined to a range extended beyond the observed predictor range by the given fractions (lower, upper; a single value is used for both ends). See 'Details'.
window	Optional integer window width, used only if extend is supplied, for single-predictor models with n > 25: restricts the foot point of each observation to the neighborhood of the corresponding sorted observations. See 'Details'.
control	Optional list of optimization settings. Recognized elements are outer_max, the total budget of Levenberg-Marquardt iterations (default 5000; with algorithm = "dense" spent in restarts of at most 1024 iterations each, see 'Details'); ftol and ptol, the relative-reduction and relative-step convergence tolerances (as in nls.lm.control; default 1e-10 each); diag, an optional length-q vector of positive parameter scaling factors (foot-point corrections are given the scale 1; if omitted, MINPACK's adaptive scaling is used); and lower/upper, alternatives to the arguments of the same name.
lower	Optional vector of lower bounds for the model parameters, of length q (same ordering as start, fixed parameters included).
upper	Optional vector of upper bounds for the model parameters, of length q (same ordering as start, fixed parameters included).
fixed	Optional logical vector with the same length and ordering as start. Parameters marked TRUE are held fixed at their supplied starting values throughout optimization. If all parameters are fixed, only the foot points are optimized (as used by profile.onls for one-parameter models).
subset	Optional logical or index expression, evaluated in data, selecting the observations to be used for fitting. Row-aligned weights, sigma_x and sigma_y are subset accordingly.
na.action	Function indicating how missing values should be handled. Row-aligned weights, sigma_x and sigma_y follow the removal of observations.
trace	Logical. If TRUE, progress messages from the NLS warm start and the orthogonal-distance optimization (iterations or restarts) are printed.
algorithm	The solver for the joint problem; both minimize the same objective. "auto" (default) uses "structured" for more than 100 observations or when parameter bounds (lower/upper) are given, and "dense" otherwise. "structured" uses the block structure of the problem, as ODRPACK does, and its cost per iteration grows linearly with n. "dense" passes the full joint problem to nls.lm (the solver of onls 0.2), with a cost per iteration growing with n^3. Fits with the foot-point bounds extend/window always use "dense". See 'Details'.

Details

Model and objective. Assume a nonlinear model y = f(x, \theta), where x is a vector of p predictors and \theta is a vector of q unknown model parameters. Unlike ordinary nonlinear least squares, which assumes error only in the response, ODR assumes that both the response and the predictors are observed with error. For observation i, let x_i be the observed predictor vector, y_i the observed response, \delta_i the unknown correction to the predictors (so that \xi_i = x_i + \delta_i is the foot point on the model surface), Qyy_i the response precision and Qx_i the predictor precision matrix (both constructed below). The parameters and all n corrections are estimated jointly by minimizing

S(\theta, \delta_1, \dots, \delta_n) = \sum_{i=1}^{n} \left[ Qyy_i \left(y_i - f(x_i + \delta_i, \theta)\right)^2 + \delta_i^T Qx_i \, \delta_i \right]

This is the explicit ODR problem of ODRPACK, with Qyy_i playing the role of the response weights (WE) and Qx_i the role of the predictor weights (WD). It has q_{free} + np unknowns, where q_{free} is the number of free (non-fixed) parameters. For fixed \theta, minimizing S over \delta_i alone gives the weighted squared orthogonal distance of observation i to the model surface,

d_i^2 = Qyy_i \left[y_i - f(\hat\xi_i, \theta) \right]^2 + (\hat\xi_i - x_i)^T Qx_i (\hat\xi_i - x_i)

so that the estimate \hat\theta is the minimizer of \sum_i d_i^2; the d_i are returned as orth_dist (residuals_o), and \sum_i d_i^2 as objective.

Predictor precision (sigma_x):
The predictor precision Qx_i used above is constructed from sigma_x, which is always supplied on the standard-deviation scale:
NULL: unit variance for every predictor, Qx_i = I_p for all i
scalar \sigma: isotropic, Qx_i = I_p/\sigma^2 for all i
length-p vector (\sigma_{x1}, \dots, \sigma_{xp}): diagonal, predictor-specific, but identical across observations,

Qx_i = \mathrm{diag}\left(1/\sigma_{x1}^2, \dots, 1/\sigma_{xp}^2\right) \text{ for all } i

p \times p matrix \Sigma_x: global covariance allowing correlated predictor errors, identical across observations,

Qx_i = \Sigma_x^{-1} \text{ for all } i

n \times p matrix (or, when p = 1, a length-n vector as shorthand): observation-specific standard deviations \sigma_{x,i1}, \dots, \sigma_{x,ip}, giving a genuinely observation-varying, diagonal-only precision

Qx_i = \mathrm{diag}\left(1/\sigma_{x,i1}^2, \dots, 1/\sigma_{x,ip}^2\right)

Cross-predictor correlation is not supported at the per-observation level (only globally, via the p \times p form above), since that would require a separate p \times p matrix per observation. In every case except the last, Qx_i is the same fixed matrix/vector for all observations and is simply written Qx. If \sigma_x is (nearly) zero for a predictor that is in fact measured without error (for example a time index), supplying a small sigma_x makes the fit approach ordinary nonlinear least squares in that direction. With the default unit precisions, the relative scale of the predictors and the response directly determines how much of the misfit is attributed to x rather than y.

Response precision (sigma_y, weights): The baseline response precision is 1/\sigma_{y,i}^2, using either the common sigma_y or, if a length-n vector is supplied, the observation-specific value. If observation weights w_i are supplied,

Qyy_i = w_i / \sigma_{y,i}^2

so that observation weights act multiplicatively on the response precision and thus enter the objective S directly (affecting both the foot-point corrections and the parameter estimates), rather than acting solely as conventional nls-style regression weights. Weights act on the response side only; the predictor precision Qx_i is not scaled by them.

Relation to ODRPACK weights (WE, WD). ODRPACK takes a response weight WE and a predictor weight WD for every observation, both as precisions (inverse variances). In onls they correspond to Qyy_i = WE_i and Qx_i = WD_i (a diagonal matrix if p > 1), and are specified as weights = WE (together with the default sigma_y = 1; alternatively sigma_y = 1/sqrt(WE) without weights) and sigma_x = 1/sqrt(WD) (for p > 1 an n \times p matrix with entries 1/\sqrt{WD_{ij}}). ODRPACK always scales the standard errors by the residual variance, which corresponds to known_sigma = FALSE (supplying sigma_x or sigma_y otherwise sets known_sigma = TRUE). A full weight matrix that differs between observations (a three-dimensional WD in ODRPACK) is not supported, only a global covariance matrix (see sigma_x). The demo gauge (doGauge) fits the ODRPACK example with separate WE and WD weights.

Algorithm.

1) Warm start. A conventional nonlinear least-squares model is fitted using nlsLM (optionally weighted by weights; fixed parameters are substituted into the model formula) to obtain starting values \theta^{(0)}. If this fit fails, a warning is issued and the raw start values are used instead.
2) Joint Levenberg-Marquardt problem. With L_i the factor satisfying L_i^T L_i = Qx_i (elementwise square root for a diagonal Qx_i, Cholesky factor for a correlated global \Sigma_x), define the residual vector of length n + np

r(\theta, \delta) = \left( \left\{ Qyy_i^{1/2}\left(y_i - f(x_i + \delta_i, \theta)\right) \right\}_{i=1}^{n}, \; \left\{ L_i \delta_i \right\}_{i=1}^{n} \right)

so that S = r^T r. Starting at (\theta^{(0)}, \delta = 0), the Levenberg-Marquardt iteration (Moré 1978; nls.lm) repeatedly solves

\left(J^T J + \mu D^2\right) \Delta z = -J^T r \qquad z = (\theta_{free}, \delta_1, \dots, \delta_n)

with a trust-region controlled damping \mu and a diagonal scaling D (the running maximum of the Jacobian column norms, or control$diag). The Jacobian of r has the sparse 'arrow' structure of ODRPACK,

J = \left( \begin{array}{cc} -W_y^{1/2} F_\theta & -W_y^{1/2} G \\ 0 & L_x \end{array} \right)

where W_y = \mathrm{diag}(Qyy_i), F_\theta is the n \times q_{free} matrix with rows \partial f(\xi_i,\theta)/\partial\theta^T, G is the n \times np block-diagonal matrix whose i-th row holds g_i^T = \partial f(\xi_i,\theta)/\partial x^T in the columns belonging to \delta_i, and L_x = \mathrm{blockdiag}(L_1,\dots,L_n). Each correction \delta_i thus affects only observation i. The two algorithms differ only in how the step is computed:
algorithm = "structured" (default for n > 100 or with parameter bounds). In the normal equations, the \delta-part is block diagonal, with one p \times p block D_i = Qyy_i g_i g_i^T + Qx_i (+ damping) per observation. These blocks are eliminated (Schur complement, as in ODRPACK), so that each step needs one q_{free} \times q_{free} solve and n small block solves (Sherman-Morrison for diagonal Qx_i, a Cholesky factorization vectorized over the observations for a correlated \Sigma_x):

\Big(A + \mu D_\theta^2 - \sum_i B_i \tilde D_i^{-1} B_i^T\Big)\Delta\theta = -g_\theta + \sum_i B_i \tilde D_i^{-1} g_{\delta,i} \qquad \Delta\delta_i = -\tilde D_i^{-1}\big(g_{\delta,i} + B_i^T \Delta\theta\big)

with A = F_\theta^T W_y F_\theta, B_i = Qyy_i \, \partial f_i/\partial\theta \, g_i^T, \tilde D_i = D_i + \mu D_{\delta,i}^2, and g_\theta, g_{\delta,i} the parts of J^T r. The damping \mu is chosen by MINPACK's trust-region rule (Moré 1978): \mu = 0 (a Gauss-Newton step) when that step lies inside the trust region, otherwise \mu such that \|D \Delta z\| matches its radius. Parameter bounds lower/upper are handled by holding a parameter whose step would leave its bounds at the bound and re-solving for the others. The cost per iteration is O(n p^3 + n q^2 + q^3), linear in n; a fit with n = 100\,000 takes seconds. For small n its per-iteration overhead in R is higher than that of the compiled minpack.lm code, so near-degenerate fits that run thousands of iterations (for example some bootstrap resamples in confint.onls) can take longer than with "dense".
algorithm = "dense". The joint problem is passed as one nonlinear least-squares problem with the full (n + np) \times (q + np) Jacobian to nls.lm, whose dense QR factorization costs O((n + q)^3) per iteration (about 15 s per fit at n = 1000). It is kept for comparison and is the only algorithm that supports the foot-point bounds extend/window.
Both algorithms give the same estimates up to the convergence tolerance on well-posed problems. For models with several local minima (for example periodic models with free periods) they can end in different minima, because their iterates differ slightly.
3) Analytic derivatives. The blocks F_\theta and G are supplied to the solver exactly, by symbolic differentiation of the model formula (deriv) when possible, and by central finite differences otherwise (for example if the formula calls a user-defined function). This matters: MINPACK's own forward-difference Jacobian uses a step proportional to the size of each unknown, which becomes smaller than the floating-point resolution of f(x_i + \delta_i) once a correction \delta_i is small, and would then drive that correction to exactly zero, i.e. to a foot point that is not orthogonal.
4) Iterations and convergence. The structured algorithm iterates until convergence (codes 1 to 3, as in nls.lm, giving convergence = 0) or until control$outer_max iterations are used (convergence = 1, with a warning). minpack.lm caps a single nls.lm call at 1024 iterations, so the dense algorithm is restarted from its last iterate (which also resets the trust region) until it converges (nls.lm codes 1 to 4), a restart no longer reduces the objective, or the budget is exhausted. Convergence is declared with the tolerances ftol and ptol (default 1e-10). If the joint optimization fails to run at all, the NLS estimates with zero corrections are returned, with a warning.

Stationarity (first-order) conditions. At a solution, \partial S/\partial \delta_i = 0 and \partial S/\partial \theta = 0 give

Qx_i (\hat\xi_i - x_i) = Qyy_i \left(y_i - f(\hat\xi_i, \hat\theta)\right) \nabla_x f(\hat\xi_i, \hat\theta) \qquad \sum_{i=1}^{n} Qyy_i \left(y_i - f(\hat\xi_i, \hat\theta)\right) \frac{\partial f(\hat\xi_i, \hat\theta)}{\partial \theta} = 0

For unit precisions the first condition states that the vector from the foot point to the observation is orthogonal to the model surface at the foot point, which is what check_o verifies, as the angle between the residual vector and the model surface after each axis has been rescaled by its precision.

Optional foot-point bounds (extend, window). By default the corrections \delta_i are unbounded, as in ODRPACK. For single-predictor models, supplying extend confines the foot points to [\text{XLOW}, \text{XUPP}],

\text{XLOW} = \min(x) - \text{extend[1]} \times \mathrm{range}(x) \qquad \text{XUPP} = \max(x) + \text{extend[2]} \times \mathrm{range}(x)

by bounding \delta_i. If, in addition, window is given and n > 25, the foot point of observation i (in sorted predictor order) is confined to [x_{(i-\text{window}+1)},\, x_{(i+\text{window}-1)}], clipped at the two ends of the data to XLOW and XUPP. With bounds, the bounded variant of the Levenberg-Marquardt algorithm is used, and a foot point sitting at a bound is generally not orthogonal. The bounded variant can also stop before the optimum, in particular with window; onls then warns when foot points are neither orthogonal nor at their bound. window has no effect unless extend is supplied; for multivariate models (p>1) neither argument is used. These bounds are rarely needed and are kept for backward compatibility and for models whose foot points must be kept on one branch of the curve.

Fixed parameters: If fixed is supplied, parameters flagged TRUE are held, throughout both the NLS warm start and the joint optimization, at their user-supplied start value (not at the NLS warm-start estimate). Fixed parameters contribute neither to the residual degrees of freedom (df_resid = n - q_{free}, with q_{free} the number of free parameters) nor to the Jacobian columns used for vcov/std_errors; their standard error is reported as 0. Observations with zero weight do not count towards df_resid.

Covariance of the parameter estimates. The covariance is the ODRPACK covariance: the \theta-block of the inverse Gauss-Newton matrix (J^T J)^{-1} of the joint problem (second-derivative terms are neglected, as in ODRPACK). Eliminating the \delta-blocks of J^T J exactly (Schur complement, using the Woodbury identity) yields the compact form used by onls, valid for any number of predictors p:

\widehat{\mathrm{Var}}(\hat\theta_{free}) = \left(F_\theta^T W F_\theta\right)^{-1} \qquad \widehat{\mathrm{Var}}(\hat\theta_{free}) = \hat\sigma^2 \left(F_\theta^T W F_\theta\right)^{-1} \quad

for known_sigma = TRUE and known_sigma = FALSE, respectively, and with F_\theta evaluated at the final foot points \hat\xi_i and estimates \hat\theta, W = \mathrm{diag}(w_1,\dots,w_n) and the effective (“effective-variance”) observation weights

w_i = \left(Qyy_i^{-1} + g_i^T \, Qx_i^{-1} \, g_i\right)^{-1} \qquad g_i = \nabla_x f(\hat\xi_i, \hat\theta)

which combine the response error and the predictor error propagated through the local model gradient (for p = 1, w_i = (1/Qyy_i + s_i^2/Qx_i)^{-1} with slope s_i). Here \hat\sigma^2 = \left(\sum_i d_i^2\right)/\text{df\_resid} is the reduced chi-square (reduced_chisq in 'Value'). Standard errors are the square roots of the diagonal.

Starting values, local minima and degenerate solutions. The objective can have several local minima (for example for periodic models with free periods, or for models that contain a limiting case such as the Richards curve, whose limit is the Gompertz curve), and the solution found depends on start. Starting values that make the model insensitive to some parameters (for example a sigmoid whose exponent underflows for all observations, so that the model is a constant) create a degenerate stationary point at which the solver may report convergence; onls then issues a warning naming the parameters that have no measurable influence on the fit, and the results should be discarded. A low orthogonal residual sum of squares is not by itself evidence of a good fit: with unit precisions and a steep model, an ODR fit can lower the objective by shifting observations horizontally. Compare the vertical residuals (residuals) with those of the NLS fit at parNLS, and choose sigma_x/sigma_y to reflect the actual measurement errors.

The returned object contains both classical vertical residual information and orthogonal-distance quantities such as foot points (xi, yi; the predictor corrections are xi - X_obs) and weighted orthogonal distances.

IMPORTANT: If not all points are orthogonal to the fitted curve, print.onls gives a “FAILED: Only X out of Y fitted points are orthogonal” message. In this case, it is suggested to conduct a more detailed analysis using check_o. Because orthogonality of the foot points is the stationarity condition \partial S/\partial \delta_i = 0, the most common cause is a convergence tolerance that is too loose for observations with very small residuals, for which the orthogonality angle is very sensitive to the foot-point position: tightening control = list(ftol = 1e-12, ptol = 1e-12) (or increasing control$outer_max if the iteration budget was exhausted) will normally resolve it. Foot points held at a bound set through extend/window are not expected to be orthogonal.

ALSO IMPORTANT: Regular R-like notation such as lm(y ~ x1 + x2) or lm(y ~ x1 * x2) needs to be changed to classical notation y = b_0 + b_1x1 + b_2x2 or y = b_0 + b_1x1 + b_2x2 + b_3x1x2, respectively.

The resulting orthogonal model houses information in respect to the (classical) vertical residuals as well as the (minimized, precision-weighted) orthogonal residuals.

The following functions use the vertical residuals:
deviance
fitted
residuals
logLik

The following functions use the orthogonal residuals:
deviance_o
residuals_o
logLik_o anova.onls

Value

An orthogonal fit of class onls (also inheriting from nls), a list with the following elements:

call, data, formula	The matched call, the data argument as given (unevaluated) and the model formula.
algorithm	The algorithm used, "structured" or "dense" (see algorithm).
coefficients	Named vector of the final parameter estimates.
std_errors, vcov	Standard errors (0 for fixed parameters) and covariance matrix of the estimates; see 'Details'.
parNLS	Parameter estimates of the initial nonlinear least-squares fit (the starting values, if the warm start failed; fixed parameters at their starting values).
convergence, convInfo	Convergence code (0 if the joint optimization converged, termination codes 1 to 4; 1 otherwise) and details: isConv, finIter (total number of Levenberg-Marquardt iterations), finTol (the ftol used), stopCode and stopMessage.
objective	The minimized objective S = \sum_i d_i^2.
reduced_chisq	Reduced chi-square S/\mathrm{df\_resid}.
known_sigma	Logical: were the variances treated as known?
df_resid	Residual degrees of freedom: the number of observations with non-zero weight minus the number of free parameters.
n	Number of observations used (after subset and na.action).
X_obs, y_obs	Observed predictors (n \times p matrix) and responses.
xi, yi	Foot points \hat\xi_i (n \times p matrix) and foot-point responses f(\hat\xi_i, \hat\theta); see x0.
fittedONLS	Model evaluated at the observed predictors, f(x_i, \hat\theta); the vertical residuals are y_obs - fittedONLS (residuals).
orth_dist	Weighted orthogonal distances d_i.
weights	Observation weights.
Q_yy, Lx	Response precisions Qyy_i (length n) and predictor precisions as an n \times p \times p array of upper Cholesky factors, L_i^T L_i = Qx_i: one representation for every form of sigma_x.
sigma_x, sigma_y, per_obs_x	sigma_x and sigma_y as used (after subset/na.action), and whether sigma_x was given per observation.
fixed, lower, upper, control, extend, window	The arguments of the same names, used when the model is refitted by profile.onls and confint.onls.
na.action, dataClasses	Information on omitted observations, and the classes of the model variables (used by predict).
ortho	Logical vector, one value per observation: the Ortho column of check_o (NA if it could not be computed).

All per-observation elements are in the original observation order. Version 0.3-0 removed the duplicated elements of earlier versions (parONLS, x0/y0, resid_o, pred/resp, residONLS, residNLS, fittedNLS, grad, QR, Q_x, Sigma_x, model, ...); use coef, x0/y0, residuals_o, residuals and fitted instead.

Author(s)

Andrej-Nikolai Spiess

References

A stable and efficient algorithm for nonlinear orthogonal distance regression.
Boggs PT, Byrd RH and Schnabel RB.
SIAM J Sci Stat Comput (1987), 8: 1052-1078.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1137/0908085")}.

The Levenberg-Marquardt algorithm: implementation and theory.
More JJ.
In: Watson GA (ed.), Numerical Analysis, Lecture Notes in Mathematics 630, Springer (1978): 105-116.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/BFb0067700")}.

Orthogonal Distance Regression.
Boggs PT and Rogers JE.
NISTIR (1990), 89-4197: 1-15.
https://static.scipy.org/doc/external/odr_ams.pdf.

User's Reference Guide for ODRPACK Version 2.01
Software for Weighted Orthogonal Distance Regression.
Boggs PT, Byrd RH, Rogers JE and Schnabel RB.
NISTIR (1992), 4834: 1-113.
https://static.scipy.org/doc/external/odrpack_guide.pdf.

ALGORITHM 676 ODRPACK: Software for Weighted Orthogonal Distance Regression.
Boggs PT, Donaldson JR, Byrd RH and Schnabel RB.
ACM Trans Math Soft (1989), 15, 348-364.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1145/76909.76913")}.

Omnivariant generalized least squares regression.
Daeron M and Vermeesch P, Chemical Geology (2024), 647: 121881.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.chemgeo.2023.121881")}.

An analysis of the Total Least Squares problem.
Golub GH and Van Loan CF.
SIAM Journal on Numerical Analysis (1980), 17, 883-893.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1137/0717073")}.

Least squares fitting of a straight line with correlated errors.
York D.
Earth and Planetary Science Letters (1968), 5, 320-324.
\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/S0012-821X(68)80059-7")}.

Examples

## 1. The DNase data from 'nls', use all generic functions.
DNase1 <- subset(DNase, Run == 1)
set.seed(123)
DNase1$density <- sapply(DNase1$density, function(x) rnorm(1, x, 0.1 * x))
mod1 <- onls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), 
             data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1))
print(mod1)
plot(mod1)
summary(mod1)
predict(mod1, newdata = data.frame(conc = 6))
logLik(mod1)
deviance(mod1)
formula(mod1)
weights(mod1)
df.residual(mod1)
fitted(mod1)
residuals(mod1)
vcov(mod1)
coef(mod1)

## 2a. Update model
DNase2 <- DNase1
DNase2$conc <- DNase2$conc * 2
mod2a <- update(mod1, data = DNase2)
print(mod2a)

## 2b. Example with a fixed parameter
## => Asym = 3.
mod2b <- onls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), 
             data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1), 
             fixed = c(TRUE, FALSE, FALSE))
print(mod2b)

## 3. Multivariate example: matched curvature, 
## low noise, decorrelated predictors 
set.seed(123)
n  <- 25
x1 <- runif(n, 1, 5)     
x2 <- runif(n, 1, 5)
b1_true <- 5
b2_true <- 2
b3_true <- 1.5
z_true <- b1_true + b2_true * x1 + b3_true * x2^2
z  <- z_true + rnorm(n, 0, 2)      
x1 <- x1 + rnorm(n, 0, 0.5)
x2 <- x2 + rnorm(n, 0, 0.5)
DAT <- data.frame(x1 = x1, x2 = x2, z = z)
mod3 <- onls(z ~ b1 + b2 * x1 + b3 * x2^2, data = DAT, 
             start = list(b1 = 1, b2 = 1, b3 = 1), trace = TRUE)
print(mod3)

## The reference tests against published ODR results, closed-form solutions, ODRPACK and
## the NIST certified values are collected in one demo: doGauge(), or
## demo("gauge", package = "onls").


## 4. The two algorithms: same fit, very different cost for large n.
set.seed(1)
n <- 400
xt <- seq(0, 10, length.out = n)
DAT <- data.frame(x = xt + rnorm(n, 0, 0.2), y = 5/(1 + exp(-(xt - 5))) + rnorm(n, 0, 0.2))
system.time(m_s <- onls(y ~ b1/(1 + exp(-(x - b2)/b3)), data = DAT,
                        start = list(b1 = 4, b2 = 4, b3 = 1.5)))   # auto: structured for n > 100
system.time(m_d <- onls(y ~ b1/(1 + exp(-(x - b2)/b3)), data = DAT,
                        start = list(b1 = 4, b2 = 4, b3 = 1.5), algorithm = "dense"))
rbind(structured = coef(m_s), dense = coef(m_d))

onls documentation built on Nov. 9, 2026, 5:06 p.m.