onls: Weighted Orthogonal Nonlinear Least-Squares Regression...
In onls: Orthogonal Nonlinear Least-Squares Regression
onls R Documentation
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, Rogers and Schnabel; 1987): 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)
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. 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.
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, spent in restarts of at most 1024 iterations each, see 'Details'); ftol and ptol, the relative-reduction and relative-step convergence tolerances passed to 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. At least one parameter must be free.
subset
Optional specification of a subset of 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 (including restarts) are printed.
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 and resid_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). See example 12.
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 scaling D. 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. ODRPACK exploits this structure to make each step cheap; onls delegates the step to the general dense solver of minpack.lm, so the cost per iteration grows faster with n than in ODRPACK.
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) Restarts and convergence. minpack.lm caps a single nls.lm call at 1024 iterations. The solver is restarted from its last iterate (which also resets the trust region) until it converges (nls.lm codes 1 to 4, giving convergence = 0), a restart no longer reduces the objective, or the budget is exhausted (convergence = 1, with a warning). 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 an angle or, for weighted fits, through the relative residual of this equation).
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. 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.
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 (residONLS) with those of the NLS fit (residNLS), 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, predictor corrections 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 Euclidean) 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
Value
An orthogonal fit of class onls with the following list items:
data
Original data argument.
call
Matched function call.
convInfo
Convergence information from the joint optimization: isConv (logical), finIter (total number of Levenberg-Marquardt iterations over all restarts), finTol (the ftol used), stopCode and stopMessage (termination code and message of the last nls.lm call).
na.action
Information on omitted observations.
dataClasses
Classes of model variables.
model
Model frame used for fitting.
formula
Model formula.
parNLS
Parameter estimates from the initial nonlinear least-squares fit (the starting values, if the warm start failed; fixed parameters at their starting values).
parONLS
Final orthogonal-distance parameter estimates.
x0, y0
Coordinates \hat\xi_i, f(\hat\xi_i,\hat\theta) of the fitted foot points on the model curve/surface, in sorted predictor order (see 'Details'; for p>1, sorting is a no-op and the original order is used).
fittedONLS
Model evaluated at the observed predictors, f(x_i, \hat\theta), using the final ODR parameter estimates.
fittedNLS
Model evaluated at the observed predictors using the initial NLS warm-start estimates.
residONLS
Vertical residuals y_i - f(x_i,\hat\theta) from the orthogonal model (observed x, not the foot point).
residNLS
Vertical residuals from the initial nonlinear least-squares fit.
resid_o
Weighted orthogonal distances d_i (see 'Details') at the final estimates, in sorted predictor order. Numerically equivalent to orth_dist, but reordered; compare to orth_dist for the original observation order.
pred, resp
Predictor and response values used internally during fitting (sorted order for p=1).
grad
Jacobian \partial f(\hat\xi_i,\theta)/\partial\theta of the model with respect to all model parameters (fixed and free), evaluated at the final foot points \hat\xi_i and estimates \hat\theta (not at the observed predictors), in the original observation order. Only its free-parameter columns are used to construct vcov.
QR
QR decomposition of grad.
weights
Observation weights used in fitting, in original observation order.
control
The resolved control settings actually passed to nls.lm (via nls.lm.control) for each restart, not the raw control argument.
coefficients
Named vector of final parameter estimates.
std_errors
Approximate standard errors of parameter estimates (0 for fixed parameters); see 'Details' for the covariance used.
vcov
Estimated parameter covariance matrix; see 'Details'.
xi
Estimated foot-point coordinates \hat\xi_i in the original observation order.
x_corrections
Estimated predictor corrections \hat\delta_i = \hat\xi_i - x_i.
orth_dist
Final weighted orthogonal distances d_i, in original observation order (see resid_o for the sorted-order equivalent).
objective
Minimized objective value S(\hat\theta, \hat\delta) = \sum_i d_i^2.
Q_x
Predictor precision structure used in fitting.
Sigma_x
Predictor covariance structure corresponding to Q_x.
Q_yy
Response precision vector.
is_diag_x
Logical indicating whether predictor precision is diagonal.
per_obs_x
Logical indicating whether predictor precision varies by observation.
sigma_x
User-supplied predictor error specification (after subset/na.action).
sigma_y
User-supplied response error specification (after subset/na.action).
known_sigma
Logical indicating whether measurement variances were treated as known.
reduced_chisq
Reduced chi-square statistic based on the orthogonal-distance criterion.
pred_names
Predictor variable names.
param_names
Model parameter names.
start
Starting parameter values.
fixed
Logical vector identifying fixed parameters.
lower, upper
Parameter bounds supplied to the optimizer.
y_obs, X_obs
Observed response and predictor values.
resp_name
Response variable name.
n, p, q
Numbers of observations, predictors and model parameters.
df_resid
Residual degrees of freedom, n minus the number of free (non-fixed) parameters.
convergence
Convergence code, 0 if the joint optimization converged (nls.lm termination codes 1 to 4), 1 otherwise.
iters
Total number of Levenberg-Marquardt iterations, summed over all restarts (NA if the joint optimization failed to run).
type
Model type, currently "nonlinear".
ortho
Orthogonality diagnostics returned by check_o (NA if they could not be computed).
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)
## Reference tests comparing to pivotal literature
## 4. Example from odrpack_guide.pdf, 2.C.i, pages 39ff.
x <- c(0, 0, 5, 7, 7.5, 10, 16, 26, 30, 34, 34.5, 100)
y <- c(1265, 1263.6, 1258, 1254, 1253, 1249.8, 1237, 1218, 1220.6,
1213.8, 1215.5, 1212)
DAT <- data.frame(x, y)
mod4 <- onls(y ~ b1 + b2 * (exp(b3 * x) -1)^2, data = DAT,
start = list(b1 = 1500, b2 = -50, b3 = -0.1))
deviance_o(mod4) # 21.445 as on page 47
summary(mod4) # 1264.65481 (1.03492) / -54.01838 (1.583992) / -0.08785 (6.33222E-3) as on page 48
## 5. Example from Algorithm 676: ODRPACK, page 355 + 356.
x <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 85, 90, 95, 100, 105)
y <- c(4.14, 8.52, 16.31, 32.18, 64.62, 98.76, 151.13, 224.74, 341.35,
423.36, 522.78, 674.32, 782.04, 920.01)
DAT <- data.frame(x, y)
mod5 <- onls(y ~ b1 * 10^(b2 * x/(b3 + x)), data = DAT,
start = list(b1 = 1, b2 = 5, b3 = 100))
deviance_o(mod5) # 15.263 as on page 363
summary(mod5) # 4.4879 (0.56876) / 7.1882 (0.69504) / 221.8383 (37.2313) as on page 363
## 6. Example with bounds from simple_example.f90
## in https://www.netlib.org/toms/869.zip.
x <- c(0.982, 1.998, 4.978, 6.01)
y <- c(2.7, 7.4, 148.0, 403.0)
DAT <- data.frame(x, y)
mod6 <- onls(y ~ b1 * exp(b2 * x), data = DAT,
start = list(b1 = 2, b2 = 0.5),
lower = c(0, 0), upper = c(10, 0.9))
coef(mod6) # 1.4376 / 0.9 ## Different to reference 1.6334 / 0.9
deviance_o(mod6) # 0.1919 => lower RSS than original ODRPACK with 0.2674!
## 7. Example similar to Deming regression
## Comparison to XLstat
## https://help.xlstat.com/6650-run-deming-regression-compare-methods-excel
x <- c(9.8, 9.7, 10.7, 10.9, 12.4, 12.5, 12.8, 12.8, 12.9, 13.3,
13.4, 13.5, 13.7, 14.9, 15.2, 15.5)
y <- c(10.1, 11.4, 10.8, 11.3, 11.8, 12.1, 12.3, 13.6, 14.2, 14.4,
14.6, 15.3, 15.5, 15.8, 16.2, 16.5)
DAT <- data.frame(x, y)
mod7 <- onls(y ~ a + b * x, data = DAT, start = list(a = 2, b = 3))
print(mod7) ## -1.909 / 1.208 as on webpage
plot(mod7)
## 8. Linear multivariate model, using the closed-form Total Least Squares
## (TLS) solution from Golub & Van Loan (1980)
tls_fit <- function(X, y) {
X <- as.matrix(X)
n <- nrow(X); p <- ncol(X)
Xc <- scale(X, center = TRUE, scale = FALSE)
yc <- y - mean(y)
xbar <- colMeans(X); ybar <- mean(y)
Z <- cbind(Xc, yc)
SVD <- svd(Z)
v <- SVD$v[, p + 1L]
v_x <- v[1:p]; v_y <- v[p + 1L]
slope <- -v_x / v_y; intercept <- ybar - sum(slope * xbar)
list(intercept = intercept, slope = setNames(slope, colnames(X)),
singular_values = SVD$d)
}
set.seed(11)
n <- 40
x1_true <- runif(n, 0, 10)
x2_true <- runif(n, 0, 10)
b0_true <- 3; b1_true <- 1.5; b2_true <- -0.8
y_true <- b0_true + b1_true * x1_true + b2_true * x2_true
x1 <- x1_true + rnorm(n, 0, 0.5)
x2 <- x2_true + rnorm(n, 0, 0.5)
y <- y_true + rnorm(n, 0, 0.5)
DAT <- data.frame(x1 = x1, x2 = x2, y = y)
TLS <- tls_fit(DAT[, c("x1", "x2")], DAT$y)
mod8 <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT, start = list(b0 = 1, b1 = 1, b2 = 1))
TLS_vec <- c(b0 = TLS$intercept, b1 = TLS$slope[["x1"]], b2 = TLS$slope[["x2"]])
ONLS_vec <- coef(mod8)[c("b0", "b1", "b2")]
print(data.frame(TLS_closed_form = TLS_vec, onls = ONLS_vec,
abs_diff = abs(TLS_vec - ONLS_vec))) # all equal
## 9. Pearson (1901) / York (1966) -> "Pearson's data with York's weights"
# Intercept: 5.47991 (SE 0.29497)
# Slope: -0.48053 (SE 0.05799)
x <- c(0.0, 0.9, 1.8, 2.6, 3.3, 4.4, 5.2, 6.1, 6.5, 7.4)
y <- c(5.9, 5.4, 4.4, 4.6, 3.5, 3.7, 2.8, 2.8, 2.4, 1.5)
sd_x <- 1/sqrt(c(1000.0, 1000.0, 500.0, 800.0, 200.0, 80.0, 60.0, 20.0, 1.8, 1.0))
sd_y <- 1/sqrt(c(1.0, 1.8, 4.0, 8.0, 20.0, 20.0, 70.0, 70.0, 100.0, 500.0))
DAT <- data.frame(x = x, y = y)
mod9 <- onls(y ~ b0 + b1*x, data = DAT,
start = list(b0 = 5, b1 = -0.5),
sigma_x = sd_x, sigma_y = sd_y)
summary(mod9) # 5.47991 (0.29497) / -0.48053 (0.05799) as in paper
## 10. Daeron & Vermeesch (2024), Table 3 / Figure 2C toy example.
x <- c(9, 19, 31, 41)
y <- c(21, 31, 39, 49)
DAT <- data.frame(x = x, y = y)
mod10 <- onls(y ~ a + b * x, data = DAT, start = list(a = 10, b = 1), sigma_x = 1, sigma_y = 1)
summary(mod10) # 13.71 / 0.851 as in Table 3 of paper
## 11. Full predictor covariance (correlated predictor errors), compared to the closed-form
## generalized Total Least Squares (TLS) solution. Whitening the predictors with the Cholesky
## factor L of the precision matrix (t(L) %*% L = solve(Sigma)) and scaling y by 1/sigma_y
## turns the problem into plain TLS, which is solved by an SVD.
set.seed(123)
n <- 40
Sigma <- matrix(c(0.25, 0.15, 0.15, 0.16), 2) # correlation of predictor errors = 0.75
sigma_y <- 0.3
xt <- cbind(runif(n, 0, 10), runif(n, 0, 10))
E <- matrix(rnorm(2 * n), n) %*% chol(Sigma)
DAT <- data.frame(x1 = xt[, 1] + E[, 1], x2 = xt[, 2] + E[, 2],
y = 3 + 1.5 * xt[, 1] - 0.8 * xt[, 2] + rnorm(n, 0, sigma_y))
mod12 <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT,
start = list(b0 = 1, b1 = 1, b2 = 1), sigma_x = Sigma, sigma_y = sigma_y,
control = list(ftol = 1e-13, ptol = 1e-13))
L <- chol(solve(Sigma))
U <- as.matrix(DAT[, c("x1", "x2")]) %*% t(L) # whitened predictors
v <- DAT$y / sigma_y
Z <- cbind(scale(U, scale = FALSE), v - mean(v))
V <- svd(Z)$v[, 3]
w <- -V[1:2] / V[3]
gTLS <- c(sigma_y * (mean(v) - sum(w * colMeans(U))), sigma_y * drop(t(L) %*% w))
print(data.frame(gen_TLS = gTLS, onls = coef(mod12),
abs_diff = abs(gTLS - coef(mod12)))) # all equal
## 12. ODRPACK's separate weights WE (response) and WD (predictor).
## ODRPACK takes both as precisions (inverse variances), per observation. In onls(), WE is
## passed as 'weights' (with the default sigma_y = 1) and WD through sigma_x = 1/sqrt(WD)
## (for p > 1: an n x p matrix 1/sqrt(WD)). known_sigma = FALSE gives ODRPACK's scaling of
## the standard errors by the residual variance.
set.seed(123)
n <- 30
xt <- seq(0.5, 10, length.out = n)
WD <- runif(n, 0.5, 4) # predictor weights
WE <- runif(n, 0.5, 4) # response weights
x <- xt + rnorm(n, 0, 0.3/sqrt(WD))
y <- 2 * exp(-0.3 * xt) + 0.5 + rnorm(n, 0, 0.03/sqrt(WE))
DAT <- data.frame(x, y)
mod12 <- onls(y ~ a * exp(-b * x) + c, data = DAT, start = list(a = 1.5, b = 0.2, c = 0.3),
weights = WE, sigma_x = 1/sqrt(WD), known_sigma = FALSE)
summary(mod12)
onls documentation built on Sept. 21, 2026, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| onls | R Documentation |
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, Rogers and Schnabel; 1987): 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)
Arguments
formula |
A two-sided nonlinear model |
data |
A data frame containing the variables appearing in |
start |
A named list (or named vector) of starting values for all model parameters. |
weights |
Optional non-negative observation weights. 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, |
sigma_x |
Specification of predictor measurement error, as a standard deviation (not a variance or a weight). Accepts |
sigma_y |
Response measurement error, as a standard deviation. May be either a single positive value applied to all observations or a length- |
known_sigma |
Logical indicating whether the supplied (or default) |
extend |
Optional numeric vector of length one or two. If |
window |
Optional integer window width, used only if |
control |
Optional list of optimization settings. Recognized elements are |
lower |
Optional vector of lower bounds for the model parameters, of length |
upper |
Optional vector of upper bounds for the model parameters, of length |
fixed |
Optional logical vector with the same length and ordering as |
subset |
Optional specification of a subset of observations to be used for fitting. Row-aligned |
na.action |
Function indicating how missing values should be handled. Row-aligned |
trace |
Logical. If |
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 and resid_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). See example 12.
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 scaling D. 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. ODRPACK exploits this structure to make each step cheap; onls delegates the step to the general dense solver of minpack.lm, so the cost per iteration grows faster with n than in ODRPACK.
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) Restarts and convergence. minpack.lm caps a single nls.lm call at 1024 iterations. The solver is restarted from its last iterate (which also resets the trust region) until it converges (nls.lm codes 1 to 4, giving convergence = 0), a restart no longer reduces the objective, or the budget is exhausted (convergence = 1, with a warning). 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 an angle or, for weighted fits, through the relative residual of this equation).
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. 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.
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 (residONLS) with those of the NLS fit (residNLS), 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, predictor corrections 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 Euclidean) 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
Value
An orthogonal fit of class onls with the following list items:
data |
Original |
call |
Matched function call. |
convInfo |
Convergence information from the joint optimization: |
na.action |
Information on omitted observations. |
dataClasses |
Classes of model variables. |
model |
Model frame used for fitting. |
formula |
Model formula. |
parNLS |
Parameter estimates from the initial nonlinear least-squares fit (the starting values, if the warm start failed; fixed parameters at their starting values). |
parONLS |
Final orthogonal-distance parameter estimates. |
x0, y0 |
Coordinates |
fittedONLS |
Model evaluated at the observed predictors, |
fittedNLS |
Model evaluated at the observed predictors using the initial NLS warm-start estimates. |
residONLS |
Vertical residuals |
residNLS |
Vertical residuals from the initial nonlinear least-squares fit. |
resid_o |
Weighted orthogonal distances |
pred, resp |
Predictor and response values used internally during fitting (sorted order for |
grad |
Jacobian |
QR |
QR decomposition of |
weights |
Observation weights used in fitting, in original observation order. |
control |
The resolved control settings actually passed to |
coefficients |
Named vector of final parameter estimates. |
std_errors |
Approximate standard errors of parameter estimates ( |
vcov |
Estimated parameter covariance matrix; see 'Details'. |
xi |
Estimated foot-point coordinates |
x_corrections |
Estimated predictor corrections |
orth_dist |
Final weighted orthogonal distances |
objective |
Minimized objective value |
Q_x |
Predictor precision structure used in fitting. |
Sigma_x |
Predictor covariance structure corresponding to |
Q_yy |
Response precision vector. |
is_diag_x |
Logical indicating whether predictor precision is diagonal. |
per_obs_x |
Logical indicating whether predictor precision varies by observation. |
sigma_x |
User-supplied predictor error specification (after |
sigma_y |
User-supplied response error specification (after |
known_sigma |
Logical indicating whether measurement variances were treated as known. |
reduced_chisq |
Reduced chi-square statistic based on the orthogonal-distance criterion. |
pred_names |
Predictor variable names. |
param_names |
Model parameter names. |
start |
Starting parameter values. |
fixed |
Logical vector identifying fixed parameters. |
lower, upper |
Parameter bounds supplied to the optimizer. |
y_obs, X_obs |
Observed response and predictor values. |
resp_name |
Response variable name. |
n, p, q |
Numbers of observations, predictors and model parameters. |
df_resid |
Residual degrees of freedom, |
convergence |
Convergence code, |
iters |
Total number of Levenberg-Marquardt iterations, summed over all restarts ( |
type |
Model type, currently |
ortho |
Orthogonality diagnostics returned by |
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)
## Reference tests comparing to pivotal literature
## 4. Example from odrpack_guide.pdf, 2.C.i, pages 39ff.
x <- c(0, 0, 5, 7, 7.5, 10, 16, 26, 30, 34, 34.5, 100)
y <- c(1265, 1263.6, 1258, 1254, 1253, 1249.8, 1237, 1218, 1220.6,
1213.8, 1215.5, 1212)
DAT <- data.frame(x, y)
mod4 <- onls(y ~ b1 + b2 * (exp(b3 * x) -1)^2, data = DAT,
start = list(b1 = 1500, b2 = -50, b3 = -0.1))
deviance_o(mod4) # 21.445 as on page 47
summary(mod4) # 1264.65481 (1.03492) / -54.01838 (1.583992) / -0.08785 (6.33222E-3) as on page 48
## 5. Example from Algorithm 676: ODRPACK, page 355 + 356.
x <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 85, 90, 95, 100, 105)
y <- c(4.14, 8.52, 16.31, 32.18, 64.62, 98.76, 151.13, 224.74, 341.35,
423.36, 522.78, 674.32, 782.04, 920.01)
DAT <- data.frame(x, y)
mod5 <- onls(y ~ b1 * 10^(b2 * x/(b3 + x)), data = DAT,
start = list(b1 = 1, b2 = 5, b3 = 100))
deviance_o(mod5) # 15.263 as on page 363
summary(mod5) # 4.4879 (0.56876) / 7.1882 (0.69504) / 221.8383 (37.2313) as on page 363
## 6. Example with bounds from simple_example.f90
## in https://www.netlib.org/toms/869.zip.
x <- c(0.982, 1.998, 4.978, 6.01)
y <- c(2.7, 7.4, 148.0, 403.0)
DAT <- data.frame(x, y)
mod6 <- onls(y ~ b1 * exp(b2 * x), data = DAT,
start = list(b1 = 2, b2 = 0.5),
lower = c(0, 0), upper = c(10, 0.9))
coef(mod6) # 1.4376 / 0.9 ## Different to reference 1.6334 / 0.9
deviance_o(mod6) # 0.1919 => lower RSS than original ODRPACK with 0.2674!
## 7. Example similar to Deming regression
## Comparison to XLstat
## https://help.xlstat.com/6650-run-deming-regression-compare-methods-excel
x <- c(9.8, 9.7, 10.7, 10.9, 12.4, 12.5, 12.8, 12.8, 12.9, 13.3,
13.4, 13.5, 13.7, 14.9, 15.2, 15.5)
y <- c(10.1, 11.4, 10.8, 11.3, 11.8, 12.1, 12.3, 13.6, 14.2, 14.4,
14.6, 15.3, 15.5, 15.8, 16.2, 16.5)
DAT <- data.frame(x, y)
mod7 <- onls(y ~ a + b * x, data = DAT, start = list(a = 2, b = 3))
print(mod7) ## -1.909 / 1.208 as on webpage
plot(mod7)
## 8. Linear multivariate model, using the closed-form Total Least Squares
## (TLS) solution from Golub & Van Loan (1980)
tls_fit <- function(X, y) {
X <- as.matrix(X)
n <- nrow(X); p <- ncol(X)
Xc <- scale(X, center = TRUE, scale = FALSE)
yc <- y - mean(y)
xbar <- colMeans(X); ybar <- mean(y)
Z <- cbind(Xc, yc)
SVD <- svd(Z)
v <- SVD$v[, p + 1L]
v_x <- v[1:p]; v_y <- v[p + 1L]
slope <- -v_x / v_y; intercept <- ybar - sum(slope * xbar)
list(intercept = intercept, slope = setNames(slope, colnames(X)),
singular_values = SVD$d)
}
set.seed(11)
n <- 40
x1_true <- runif(n, 0, 10)
x2_true <- runif(n, 0, 10)
b0_true <- 3; b1_true <- 1.5; b2_true <- -0.8
y_true <- b0_true + b1_true * x1_true + b2_true * x2_true
x1 <- x1_true + rnorm(n, 0, 0.5)
x2 <- x2_true + rnorm(n, 0, 0.5)
y <- y_true + rnorm(n, 0, 0.5)
DAT <- data.frame(x1 = x1, x2 = x2, y = y)
TLS <- tls_fit(DAT[, c("x1", "x2")], DAT$y)
mod8 <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT, start = list(b0 = 1, b1 = 1, b2 = 1))
TLS_vec <- c(b0 = TLS$intercept, b1 = TLS$slope[["x1"]], b2 = TLS$slope[["x2"]])
ONLS_vec <- coef(mod8)[c("b0", "b1", "b2")]
print(data.frame(TLS_closed_form = TLS_vec, onls = ONLS_vec,
abs_diff = abs(TLS_vec - ONLS_vec))) # all equal
## 9. Pearson (1901) / York (1966) -> "Pearson's data with York's weights"
# Intercept: 5.47991 (SE 0.29497)
# Slope: -0.48053 (SE 0.05799)
x <- c(0.0, 0.9, 1.8, 2.6, 3.3, 4.4, 5.2, 6.1, 6.5, 7.4)
y <- c(5.9, 5.4, 4.4, 4.6, 3.5, 3.7, 2.8, 2.8, 2.4, 1.5)
sd_x <- 1/sqrt(c(1000.0, 1000.0, 500.0, 800.0, 200.0, 80.0, 60.0, 20.0, 1.8, 1.0))
sd_y <- 1/sqrt(c(1.0, 1.8, 4.0, 8.0, 20.0, 20.0, 70.0, 70.0, 100.0, 500.0))
DAT <- data.frame(x = x, y = y)
mod9 <- onls(y ~ b0 + b1*x, data = DAT,
start = list(b0 = 5, b1 = -0.5),
sigma_x = sd_x, sigma_y = sd_y)
summary(mod9) # 5.47991 (0.29497) / -0.48053 (0.05799) as in paper
## 10. Daeron & Vermeesch (2024), Table 3 / Figure 2C toy example.
x <- c(9, 19, 31, 41)
y <- c(21, 31, 39, 49)
DAT <- data.frame(x = x, y = y)
mod10 <- onls(y ~ a + b * x, data = DAT, start = list(a = 10, b = 1), sigma_x = 1, sigma_y = 1)
summary(mod10) # 13.71 / 0.851 as in Table 3 of paper
## 11. Full predictor covariance (correlated predictor errors), compared to the closed-form
## generalized Total Least Squares (TLS) solution. Whitening the predictors with the Cholesky
## factor L of the precision matrix (t(L) %*% L = solve(Sigma)) and scaling y by 1/sigma_y
## turns the problem into plain TLS, which is solved by an SVD.
set.seed(123)
n <- 40
Sigma <- matrix(c(0.25, 0.15, 0.15, 0.16), 2) # correlation of predictor errors = 0.75
sigma_y <- 0.3
xt <- cbind(runif(n, 0, 10), runif(n, 0, 10))
E <- matrix(rnorm(2 * n), n) %*% chol(Sigma)
DAT <- data.frame(x1 = xt[, 1] + E[, 1], x2 = xt[, 2] + E[, 2],
y = 3 + 1.5 * xt[, 1] - 0.8 * xt[, 2] + rnorm(n, 0, sigma_y))
mod12 <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT,
start = list(b0 = 1, b1 = 1, b2 = 1), sigma_x = Sigma, sigma_y = sigma_y,
control = list(ftol = 1e-13, ptol = 1e-13))
L <- chol(solve(Sigma))
U <- as.matrix(DAT[, c("x1", "x2")]) %*% t(L) # whitened predictors
v <- DAT$y / sigma_y
Z <- cbind(scale(U, scale = FALSE), v - mean(v))
V <- svd(Z)$v[, 3]
w <- -V[1:2] / V[3]
gTLS <- c(sigma_y * (mean(v) - sum(w * colMeans(U))), sigma_y * drop(t(L) %*% w))
print(data.frame(gen_TLS = gTLS, onls = coef(mod12),
abs_diff = abs(gTLS - coef(mod12)))) # all equal
## 12. ODRPACK's separate weights WE (response) and WD (predictor).
## ODRPACK takes both as precisions (inverse variances), per observation. In onls(), WE is
## passed as 'weights' (with the default sigma_y = 1) and WD through sigma_x = 1/sqrt(WD)
## (for p > 1: an n x p matrix 1/sqrt(WD)). known_sigma = FALSE gives ODRPACK's scaling of
## the standard errors by the residual variance.
set.seed(123)
n <- 30
xt <- seq(0.5, 10, length.out = n)
WD <- runif(n, 0.5, 4) # predictor weights
WE <- runif(n, 0.5, 4) # response weights
x <- xt + rnorm(n, 0, 0.3/sqrt(WD))
y <- 2 * exp(-0.3 * xt) + 0.5 + rnorm(n, 0, 0.03/sqrt(WE))
DAT <- data.frame(x, y)
mod12 <- onls(y ~ a * exp(-b * x) + c, data = DAT, start = list(a = 1.5, b = 0.2, c = 0.3),
weights = WE, sigma_x = 1/sqrt(WD), known_sigma = FALSE)
summary(mod12)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.