Orthogonal Nonlinear Least-Squares Regression with onls"

In onls: Orthogonal Nonlinear Least-Squares Regression

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 6,
  fig.height = 5.5
)

library(onls)

1. Why orthogonal regression?

Ordinary nonlinear least squares (nls()) assumes that the predictor $x$ is known exactly and all the "noise" lives in the response $y$: it fits $y = f(x, \theta) + \varepsilon$ by minimizing the vertical distance between each point and the curve,

$$\min_{\theta} \sum_{i=1}^n \left[y_i - f(x_i, \theta)\right]^2.$$

In many real measurement situations -- calibration curves, instrument comparisons, physical or chemical assays -- both $x$ and $y$ carry measurement error. Orthogonal distance regression (ODR), also called errors-in-variables regression, accounts for this by allowing the model to also adjust the predictor: for each observation $i$, a foot point $\xi_i$ is found on the fitted curve/surface, and the model is fit by minimizing the distance from $(x_i, y_i)$ to $(\xi_i, f(\xi_i, \theta))$ instead of the purely vertical distance to $(x_i, f(x_i, \theta))$.

For a single, unweighted predictor this is literally the shortest (perpendicular) distance from each point to the curve -- hence "orthogonal" regression. onls() generalizes this to multiple predictors and to heteroscedastic/weighted measurement error in both $x$ and $y$, and it does so in the manner of ODRPACK (Boggs, Byrd, Rogers & Schnabel): same objective, same weights, same covariance estimate.

Deming regression (linear, known ratio of $x$/$y$ variances) and Total Least Squares (linear, unweighted, multivariate) are both special cases of what onls() fits; we'll see both later.

1.1 When is it worth it?

Classical least squares is biased when the predictors are measured with error. The coefficients of noisy predictors are attenuated towards zero; with correlated predictors the bias spills over to the other coefficients (even to those of nearly error-free predictors); and a nonlinear term such as $x_4^2$ is hit harder, because the error enters through the nonlinearity. Orthogonal regression puts the error where it belongs -- on the measured predictor -- and propagates it through the model function, which removes most of this bias. Two caveats matter in practice:

• the error standard deviations have to be supplied (sigma_x, sigma_y); they cannot be estimated from the data alone, and the result depends on them;
• the estimates have a larger variance than their least-squares counterparts, and they describe the relationship between the true quantities -- if the goal is merely to predict $y$ from future noisy $x$, ordinary least squares is the better predictor.

A sensible workflow is therefore to fit both and compare: if the estimates hardly differ, measurement error does not matter for your problem.

2. The mathematics

2.1 The joint criterion

For a nonlinear model $y = f(x, \theta)$ with predictor vector $x \in \mathbb{R}^p$ and parameter vector $\theta$, let $Qyy_i$ be the response precision (inverse variance) for observation $i$, and $Qx_i$ its $p \times p$ predictor precision matrix. onls() uses the formulation of ODRPACK's explicit ODR problem: the parameters and one correction $\delta_i = \xi_i - x_i$ per observation (so that $\xi_i = x_i + \delta_i$ is the foot point on the model surface) are estimated simultaneously by minimizing

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

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 the estimate $\hat\theta$ is the minimizer of $\sum_i d_i^2$. The important point is that onls() does not solve a separate foot-point problem for every trial value of $\theta$: all $q + np$ unknowns ($q$ free parameters plus $n$ corrections of length $p$) are optimized at once. This is the "combined" approach of ODRPACK.

2.2 Where the precisions come from

$Qyy_i$ and $Qx_i$ are built from the sigma_y, sigma_x, and weights arguments to onls():

• sigma_y (a standard deviation, not a variance): either one value shared by all observations, or a length-$n$ vector of observation-specific response error.
• weights: optional non-negative weights, combined multiplicatively, $Qyy_i = w_i / \sigma_{y,i}^2$. They act on the response side only.
• sigma_x: predictor measurement error, accepted as NULL (unit variance), a scalar (isotropic), a length-$p$ vector (diagonal, predictor-specific), a full $p \times p$ covariance matrix (allowing correlated predictor errors, $Qx = \Sigma_x^{-1}$), or an $n \times p$ matrix of observation-specific standard deviations.

When neither sigma_x, sigma_y, nor weights is supplied, $Qyy_i = 1$ and $Qx_i = I_p$ for every observation, and $d_i$ reduces to the plain Euclidean distance from $(x_i, y_i)$ to the curve -- the classical, unweighted orthogonal-regression case.

ODRPACK's own weights WE (response) and WD (predictor) are precisions, so they correspond to $Qyy_i = WE_i$ and $Qx_i = WD_i$; Section 3.4 shows how to pass them to onls().

2.3 The algorithm, briefly

1. An ordinary (vertical) nonlinear fit via Levenberg-Marquardt (minpack.lm::nlsLM()) gives warm-start parameter values.
2. The joint problem is written as a nonlinear least-squares problem in the unknowns $(\theta, \delta)$ with the residual vector of length $n + np$

$$r(\theta, \delta) = \Big(\big{Qyy_i^{1/2}\,[y_i - f(x_i + \delta_i, \theta)]\big}{i=1}^{n},\; \big{L_i\,\delta_i\big}{i=1}^{n}\Big), \qquad L_i^T L_i = Qx_i,$$

so that $S = r^T r$, and is solved by a single Levenberg-Marquardt run (minpack.lm::nls.lm()), starting at the warm-start values with all $\delta_i = 0$. The Jacobian of $r$ has the sparse "arrow" structure of ODRPACK, because each $\delta_i$ only affects observation $i$:

$$J = \begin{pmatrix} -W_y^{1/2} F_\theta & -W_y^{1/2} G \ 0 & L_x \end{pmatrix},$$

with $F_\theta$ the derivatives of $f$ with respect to the parameters and $G$ the (block diagonal) derivatives with respect to the predictors, both at the current foot points. onls() supplies them exactly, by symbolic differentiation of the model formula and, if the formula cannot be differentiated (for example because it calls a user-defined function), by central finite differences. 3. nls.lm() stops after 1024 iterations per call, so onls() restarts it from its last iterate until it converges, a restart brings no further improvement, or the total budget control$outer_max (default 5000 iterations) is used up. The convergence tolerances ftol and ptol are settable through control (default 1e-10). 4. The parameter covariance is the ODRPACK (Gauss-Newton) covariance, i.e. the parameter block of $(J^T J)^{-1}$, which can be written with effective-variance weights as

$$\widehat{\mathrm{Var}}(\hat\theta) = \hat\sigma^2\left(F_\theta^T W F_\theta\right)^{-1},\quad w_i = \left(Qyy_i^{-1} + g_i^T Qx_i^{-1} g_i\right)^{-1},\quad g_i = \nabla_x f(\hat\xi_i, \hat\theta).$$

Here $\hat\sigma^2$ is the reduced chi-square when known_sigma = FALSE, and $1$ when the supplied precisions are taken as known (known_sigma = TRUE, the default as soon as sigma_x or sigma_y is given). 5. If some free parameters have no measurable influence on the fitted model (for example because a sigmoid has collapsed to a constant), onls() issues a warning -- see Section 8.3.

2.4 Checking orthogonality

At the solution, $\partial S / \partial \delta_i = 0$ gives the stationarity (KKT) condition

$$Qx_i\,(\hat\xi_i - x_i) = Qyy_i\,\big(y_i - f(\hat\xi_i, \hat\theta)\big)\,\nabla_x f(\hat\xi_i, \hat\theta).$$

For unit precisions this says that the vector from the foot point to the observation is orthogonal to the model surface. check_o() verifies it after fitting -- either via the classical tangent-angle criterion (unweighted case, where a right angle is literally expected; the angle is computed with atan2() so that points whose foot point coincides with the observation are handled correctly) or via the relative residual of the stationarity condition (weighted case, where a plain right angle is no longer the correct geometric picture). The tolerances are tol_deg (default $0.05^\circ$) and tol_kkt (default $0.001$). We'll use it throughout.

3. A univariate example

We start with the classic DNase enzyme-linked immunosorbent assay data (from base R's nls documentation) and a 3-parameter logistic model -- a single predictor, no weighting, the simplest case where $d_i$ is literally the Euclidean distance to the curve.

DNase1 <- subset(DNase, Run == 1)
set.seed(1)
DNase1$density <- sapply(DNase1$density, function(x) rnorm(1, x, 0.1 * x))

mod_uni <- onls(density ~ Asym / (1 + exp((xmid - log(conc)) / scal)),
                data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1))
print(mod_uni)

summary(mod_uni)

The plot below shows the fitted orthogonal curve (red) alongside the ordinary vertical-least-squares warm-start curve (blue), with segments connecting each observation to its foot point. Because this fit is unweighted, plot.onls() uses a literal 1:1 axis aspect ratio by default -- the segments should genuinely look perpendicular to the red curve.

plot(mod_uni)

check_o() confirms this numerically: for every observation, the angle between the tangent to the curve and the line to the observation should be within $0.05^\circ$ of a right angle.

check_o(mod_uni, plot = FALSE)

3.1 Foot points and orthogonal residuals

The foot points $(\hat\xi_i, f(\hat\xi_i, \hat\theta))$ are available through x0() and y0() (for single-predictor models in the sorted predictor order that check_o() uses as well). Together with the observations they give the orthogonal distances by hand:

FP <- data.frame(x = mod_uni$pred, x0 = x0(mod_uni), y = mod_uni$resp, y0 = y0(mod_uni))
FP$dist <- sqrt((FP$x - FP$x0)^2 + (FP$y - FP$y0)^2)
head(FP)

## for unit precisions, the sum of squared distances is the minimized objective
all.equal(sum(FP$dist^2), deviance_o(mod_uni))

residuals_o() returns the fitted, precision-weighted orthogonal distances $\hat d_i$ -- the per-observation quantities that deviance_o() squares and sums -- directly, and in the original observation order (whereas x0()/y0() follow the sorted predictor order for a single predictor). For an unweighted single-predictor fit such as this one they are plain Euclidean distances, i.e. the dist column above in a different order:

d_o <- residuals_o(mod_uni)
all.equal(sort(d_o), sort(FP$dist), check.attributes = FALSE)
all.equal(sum(d_o^2), deviance_o(mod_uni))

## side by side with the vertical residuals of the same fit
head(data.frame(vertical = residuals(mod_uni), orthogonal = d_o))

If observations were dropped through na.action = na.exclude, residuals_o() re-inserts them as NA at their original positions, so that the result always lines up with the rows of the data:

DNase1_na <- DNase1
DNase1_na$density[c(3, 10)] <- NA
mod_na <- onls(density ~ Asym / (1 + exp((xmid - log(conc)) / scal)),
               data = DNase1_na, start = list(Asym = 3, xmid = 0, scal = 1),
               na.action = na.exclude)
residuals_o(mod_na)

3.2 Plot options

plot.onls() draws the observations, the orthogonal fit (red), the ordinary vertical-least-squares warm-start fit (blue) and the segments to the foot points. Each of these can be switched off (fitted.onls, fitted.nls, segments), and npoints sets the resolution of the drawn curves. The orthogonality is easiest to judge for a steep curve. Here is a noisy quadratic:

set.seed(123)
x <- 1:20
y <- 10 + 3 * x^2 + rnorm(20, 0, 50)
DAT_quad <- data.frame(x, y)
mod_quad <- onls(y ~ a + b * x^2, data = DAT_quad, start = list(a = 10, b = 3))
plot(mod_quad)

To zoom into a region, supplying xlim alone is enough: a matching ylim is chosen automatically. A strict 1:1 axis ratio is generally not possible for a zoomed region, so it has to be switched off with asp = FALSE:

plot(mod_quad, fitted.nls = FALSE, xlim = c(0, 10), asp = FALSE)

3.3 Adding measurement error: weighted regression

If the predictor and/or response measurement error is actually known (e.g. from instrument specifications or replicate measurements), supplying sigma_x/sigma_y gives a proper weighted orthogonal fit rather than treating all points as equally uncertain. A classic illustration is Pearson's (1901) data with York's (1966) weights -- ten points with wildly different, known per-point standard deviations in both $x$ and $y$:

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, 1000, 500, 800, 200, 80, 60, 20, 1.8, 1.0))
sd_y <- 1 / sqrt(c(1.0, 1.8, 4.0, 8.0, 20, 20, 70, 70, 100, 500))
DAT_py <- data.frame(x = x, y = y)

mod_w <- onls(y ~ b0 + b1 * x, data = DAT_py, start = list(b0 = 5, b1 = -0.5),
              sigma_x = sd_x, sigma_y = sd_y)
summary(mod_w)   # intercept 5.480 (0.295), slope -0.481 (0.058), matching York's published values

Since this fit is weighted, check_o() automatically switches from the tangent-angle criterion to the underlying KKT-residual criterion (a plain right angle is no longer the geometrically correct picture once the axes are rescaled by unequal precisions):

check_o(mod_w, plot = FALSE)

3.4 ODRPACK's weights WE and WD

ODRPACK takes a response weight WE and a predictor weight WD for every observation, both as precisions (inverse variances). In onls(), WE is passed as weights (with the default sigma_y = 1, or equivalently as sigma_y = 1/sqrt(WE)), and WD through sigma_x = 1/sqrt(WD) (for $p > 1$ an $n \times p$ matrix of $1/\sqrt{WD_{ij}}$). known_sigma = FALSE gives ODRPACK's scaling of the standard errors by the residual variance.

set.seed(7)
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_we <- data.frame(x, y)

mod_we <- onls(y ~ a * exp(-b * x) + c, data = DAT_we,
               start = list(a = 1.5, b = 0.2, c = 0.3),
               weights = WE, sigma_x = 1 / sqrt(WD), known_sigma = FALSE)
summary(mod_we)   # 1.98131 (0.03347) / 0.28984 (0.01468) / 0.48480 (0.02986), as scipy.odr with we = WE, wd = WD

The minimized objective is exactly the ODRPACK objective $\sum WE\,\varepsilon^2 + WD\,\delta^2$:

f <- function(x, b) b[1] * exp(-b[2] * x) + b[3]
xi <- mod_we$xi[, 1]
all.equal(sum(WE * (y - f(xi, coef(mod_we)))^2 + WD * (xi - x)^2), mod_we$objective)

4. Reference problems from the literature

onls() reproduces the published results of the standard ODRPACK test problems.

4.1 ODRPACK guide, example 2.C.i

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_guide <- data.frame(x, y)

mod_guide <- onls(y ~ b1 + b2 * (exp(b3 * x) - 1)^2, data = DAT_guide,
                  start = list(b1 = 1500, b2 = -50, b3 = -0.1))
deviance_o(mod_guide)   # 21.445, as on page 47 of the guide
summary(mod_guide)      # 1264.65481 (1.03492) / -54.01838 (1.583992) / -0.08785 (6.33222E-3), as on page 48

This example is also a good test for check_o(). The model has slope zero at $x = 0$ (the first two observations), and its last point lies far out on an almost flat part of the curve. For the two observations at $x = 0$ the foot point coincides with the observation, so the residual is purely vertical and the tangent-angle criterion has to be evaluated without dividing by $x_i - \xi_i = 0$:

check_o(mod_guide, plot = FALSE)

4.2 Algorithm 676 (ODRPACK), pages 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_676 <- data.frame(x, y)

mod_676 <- onls(y ~ b1 * 10^(b2 * x / (b3 + x)), data = DAT_676,
                start = list(b1 = 1, b2 = 5, b3 = 100))
deviance_o(mod_676)   # 15.263, as on page 363
summary(mod_676)      # 4.4879 (0.56876) / 7.1882 (0.69504) / 221.8383 (37.2313), as on page 363

4.3 Daeron & Vermeesch (2024), Table 3

A four-point toy example with unit measurement errors in both variables:

DAT_dv <- data.frame(x = c(9, 19, 31, 41), y = c(21, 31, 39, 49))
mod_dv <- onls(y ~ a + b * x, data = DAT_dv, start = list(a = 10, b = 1),
               sigma_x = 1, sigma_y = 1)
summary(mod_dv)   # 13.71 / 0.8516 (the exact TLS slope); Table 3 of the paper lists 13.71 / 0.851

5. A bivariate example, in 3D

With two predictors, the foot point $\xi_i$ is a point in a genuine 3D surface, not a 2D curve, and the orthogonal distance is the literal shortest distance from each observation to that surface. We fit a hyperboloid of one sheet,

$$z = c\sqrt{1 + (x_1/a)^2 + (x_2/b)^2},$$

with measurement error in both predictors. This surface is always real-valued (the argument under the square root is never less than 1), so it's a numerically well-behaved choice for a 3D example.

set.seed(2024)
n  <- 60
x1 <- runif(n, -5, 5)
x2 <- runif(n, -5, 5)

a_true <- 3; b_true <- 2; c_true <- 4
z <- c_true * sqrt(1 + (x1 / a_true)^2 + (x2 / b_true)^2) + rnorm(n, 0, 0.3)

x1 <- x1 + rnorm(n, 0, 0.2)
x2 <- x2 + rnorm(n, 0, 0.15)
DAT_hyp <- data.frame(x1 = x1, x2 = x2, z = z)

mod_hyp <- onls(z ~ c * sqrt(1 + (x1 / a)^2 + (x2 / b)^2), data = DAT_hyp,
                start = list(a = 2, b = 2, c = 3),
                sigma_x = c(0.2, 0.15), sigma_y = 0.3)
summary(mod_hyp)   # expect a, b, c close to 3, 2, 4

check_o(mod_hyp, plot = FALSE)

For exactly two predictors, plot.onls() produces an rgl-based 3D plot: the fitted surface (with nmesh x nmesh mesh divisions), the observations, and segments connecting each observation to its foot point on the surface. Its advantage is that it can be rotated and zoomed, which is the best way to inspect the orthogonality of the points by eye. This isn't run when this vignette is built (it opens an interactive graphics device), but works from an interactive R session:

plot(mod_hyp)

5.1 A half-dome with a parameter bound

Some models can only be evaluated on part of the parameter space. A half-dome of radius $r$,

$$z = \sqrt{r^2 - x_1^2 - x_2^2},$$

is only defined where $x_1^2 + x_2^2 \le r^2$, so $r$ must exceed the largest radius of any observation. A lower bound on the parameter guarantees this:

set.seed(123)
n <- 60
r_true <- 6
ang <- runif(n, 0, 2 * pi)
rad <- sqrt(runif(n, 0, 0.55)) * r_true
x1  <- rad * cos(ang)
x2  <- rad * sin(ang)
z <- sqrt(r_true^2 - x1^2 - x2^2) + rnorm(n, 0, 0.15)
x1 <- x1 + rnorm(n, 0, 0.1)
x2 <- x2 + rnorm(n, 0, 0.1)
DAT_dome <- data.frame(x1 = x1, x2 = x2, z = z)
maxrad <- max(sqrt(x1^2 + x2^2))

mod_dome <- onls(z ~ sqrt(r^2 - x1^2 - x2^2), data = DAT_dome,
                 start = list(r = r_true),
                 sigma_x = c(0.1, 0.1), sigma_y = 0.15,
                 lower = maxrad * 1.05, upper = 100)
summary(mod_dome)   # r close to 6
check_o(mod_dome, plot = FALSE)   # all orthogonal to the dome surface

plot(mod_dome)   # renders as a visibly round dome

6. A multivariate example -- and translating from lm()-style notation

onls() supports any number of predictors, but its formula must be written in fully explicit classical notation, with a distinct parameter symbol multiplying each term -- not the implicit lm()-style shorthand you may be used to. This matters because onls() has to identify each parameter by name (to differentiate with respect to it, to allow it to be fixed, bounded, etc.), which R's usual formula shorthand doesn't expose directly.

Suppose you wanted to relate a response $z$ to four predictors, each entering through a different functional form. In lm(), you might reach for something like:

## This is NOT valid onls() syntax -- shown only for comparison.
lm(z ~ x1 + I(x2^2) + sqrt(x3) + log(x4 + 1))

onls() needs the same relationship written out with an explicit coefficient on every term:

$$z = \beta_0 + \beta_1 x_1 + \beta_2 x_2^2 + \beta_3 \sqrt{x_3} + \beta_4 \log(x_4 + 1),$$

which translates directly into:

z ~ b0 + b1 * x1 + b2 * x2^2 + b3 * sqrt(x3) + b4 * log(x4 + 1)

Note that a transformed predictor such as $x_2^2$ is not a separate variable with its own error: the measurement error sits on $x_2$ and is propagated through the square by the model function, which is exactly what the foot-point formulation does.

Let's fit exactly that model, with measurement error on all four predictors (a diagonal sigma_x, i.e. no cross-predictor correlation -- see Section 6.1 for the fully correlated, matrix case):

set.seed(99)
n  <- 60
x1 <- runif(n, 0, 10)
x2 <- runif(n, 0, 5)
x3 <- runif(n, 2, 10)   # kept away from 0: sqrt() needs non-negative arguments
x4 <- runif(n, 2, 10)   # kept away from -1: log(x4 + 1) needs x4 + 1 > 0

b0 <- 2; b1 <- 0.8; b2 <- 0.5; b3 <- 2; b4 <- 3
z <- b0 + b1 * x1 + b2 * x2^2 + b3 * sqrt(x3) + b4 * log(x4 + 1) + rnorm(n, 0, 0.5)

sd_x <- c(0.3, 0.2, 0.3, 0.3)
x1 <- x1 + rnorm(n, 0, sd_x[1])
x2 <- x2 + rnorm(n, 0, sd_x[2])
x3 <- x3 + rnorm(n, 0, sd_x[3])
x4 <- x4 + rnorm(n, 0, sd_x[4])
DAT_mv <- data.frame(x1 = x1, x2 = x2, x3 = x3, x4 = x4, z = z)

mod_mv <- onls(z ~ b0 + b1 * x1 + b2 * x2^2 + b3 * sqrt(x3) + b4 * log(x4 + 1),
               data = DAT_mv,
               start = list(b0 = 1, b1 = 1, b2 = 1, b3 = 1, b4 = 1),
               sigma_x = sd_x, sigma_y = 0.5)
summary(mod_mv)   # expect b0..b4 close to 2, 0.8, 0.5, 2, 3

check_o(mod_mv, plot = FALSE)

With more than two predictors, plot.onls() draws a grid of partial-dependence panels by default -- one per predictor, with the others held at their mean foot-point value. These small panels are not a literally faithful orthogonality check: the drawn curve fixes the other predictors at a shared mean value, whereas each point's own segment ends at its own foot point in every dimension, so the segments should not be expected to look exactly perpendicular even for an unweighted fit. Use check_o() for a reliable per-observation, per-axis check:

plot(mod_mv)

Passing panel renders a single predictor as one full-size plot instead, with the same 1:1 aspect treatment as the univariate case -- by position or by name:

plot(mod_mv, panel = "x2")

6.1 Correlated predictor errors: a full covariance matrix

If the predictor errors are correlated, sigma_x can be a full $p \times p$ covariance matrix $\Sigma_x$. For a linear model this problem has a closed-form solution, the generalized Total Least Squares (TLS) estimate: whitening the predictors with the Cholesky factor $L$ of the precision matrix ($L^T L = \Sigma_x^{-1}$) and scaling $y$ by $1/\sigma_y$ turns it into plain TLS, which is solved by a singular value decomposition. This gives an exact reference for onls():

set.seed(2026)
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_cov <- 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))

mod_cov <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT_cov,
                start = list(b0 = 1, b1 = 1, b2 = 1), sigma_x = Sigma, sigma_y = sigma_y,
                control = list(ftol = 1e-13, ptol = 1e-13))

## closed-form generalized TLS
L <- chol(solve(Sigma))
U <- as.matrix(DAT_cov[, c("x1", "x2")]) %*% t(L)      # whitened predictors
v <- DAT_cov$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(mod_cov),
                 abs_diff = abs(gTLS - coef(mod_cov))))   # all equal

Ignoring the correlation (a diagonal sigma_x only) gives visibly different estimates:

mod_cov_d <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT_cov,
                  start = list(b0 = 1, b1 = 1, b2 = 1), sigma_x = sqrt(diag(Sigma)),
                  sigma_y = sigma_y)
coef(mod_cov_d)

7. Special cases

Two well-known linear methods are recovered exactly as special cases of onls().

7.1 Deming regression

Deming regression -- a linear model with known (or assumed) ratio of the predictor/response error variances -- is just onls() with a linear formula and appropriate sigma_x/sigma_y, as in the Pearson-York example above. With the default equal error variances it is orthogonal regression; the example below reproduces the XLSTAT Deming regression demonstration (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_dem <- data.frame(x, y)

mod_dem <- onls(y ~ a + b * x, data = DAT_dem, start = list(a = 2, b = 3))
print(mod_dem)   # -1.909 / 1.208 as on the webpage

7.2 Total Least Squares

Total Least Squares (unweighted, multivariate, linear) is recovered by an unweighted linear onls() fit. Here we compare it against the closed-form SVD solution of Golub & Van Loan (1980):

tls_fit <- function(X, y) {
  X <- as.matrix(X)
  p <- ncol(X)
  Xc <- scale(X, center = TRUE, scale = FALSE)
  yc <- y - mean(y)
  xbar <- colMeans(X); ybar <- mean(y)
  SVD <- svd(cbind(Xc, yc))
  v <- SVD$v[, p + 1L]
  slope <- -v[1:p] / v[p + 1L]
  list(intercept = ybar - sum(slope * xbar), slope = setNames(slope, colnames(X)))
}

set.seed(11)
n <- 40
x1_true <- runif(n, 0, 10)
x2_true <- runif(n, 0, 10)
y_true  <- 3 + 1.5 * x1_true - 0.8 * x2_true
DAT_tls <- data.frame(x1 = x1_true + rnorm(n, 0, 0.5),
                      x2 = x2_true + rnorm(n, 0, 0.5),
                      y  = y_true  + rnorm(n, 0, 0.5))

TLS <- tls_fit(DAT_tls[, c("x1", "x2")], DAT_tls$y)
mod_tls <- onls(y ~ b0 + b1 * x1 + b2 * x2, data = DAT_tls,
                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(mod_tls)[c("b0", "b1", "b2")]
print(data.frame(TLS_closed_form = TLS_vec, onls = ONLS_vec,
                 abs_diff = abs(TLS_vec - ONLS_vec)))   # equal to solver tolerance

8. Bounds, fixed parameters and control settings

8.1 Parameter bounds

lower and upper bound the model parameters (one value per parameter, in the order of start). This example with bounds comes from the simple_example.f90 of TOMS 869 (https://www.netlib.org/toms/869.zip); Section 5.1 shows a bound that keeps a model inside its domain:

DAT_bnd <- data.frame(x = c(0.982, 1.998, 4.978, 6.01),
                      y = c(2.7, 7.4, 148.0, 403.0))
mod_bnd <- onls(y ~ b1 * exp(b2 * x), data = DAT_bnd,
                start = list(b1 = 2, b2 = 0.5),
                lower = c(0, 0), upper = c(10, 0.9))
coef(mod_bnd)          # 1.4376 / 0.9, different to the reference 1.6334 / 0.9
deviance_o(mod_bnd)    # 0.1919, lower than the 0.2674 of the original ODRPACK

8.2 Fixed parameters

Parameters flagged in fixed are held at their starting values throughout. They are automatically excluded from the degrees of freedom, standard errors and correlation matrices (their standard error is reported as zero). Here the asymptote of the DNase model is fixed to 3:

mod_fix <- 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(mod_fix)

8.3 Control settings, starting values and degenerate solutions

Convergence is governed by control: the tolerances ftol/ptol (default 1e-10) and the total iteration budget outer_max (default 5000). The convergence details of a fit are stored in convInfo. Tighter tolerances sharpen the orthogonality angles reported by check_o(), which matters mostly for observations with very small residuals:

mod_ctrl <- onls(y ~ b1 + b2 * (exp(b3 * x) - 1)^2, data = DAT_guide,
                 start = list(b1 = 1500, b2 = -50, b3 = -0.1),
                 control = list(ftol = 1e-12, ptol = 1e-12, outer_max = 2000))
coef(mod_ctrl)
mod_ctrl$convInfo$isConv
mod_ctrl$convInfo$finIter   # total Levenberg-Marquardt iterations

Like all nonlinear fits, the result depends on the starting values, and an orthogonal fit can have several local minima. The worst case is a start at which the model is insensitive to some of its parameters. Consider the Richards-type growth curve $y = b_1 / (1 + \exp(b_2 - b_3 x))^{1/b_4}$: with the start $b_3 = 7$, $\exp(b_2 - b_3 x)$ is practically zero for every $x \ge 1$, the model is a constant, and the solver "converges" to that constant. onls() detects this and warns:

x <- 1:15
y <- c(16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53,
       590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41)
DAT_rich <- data.frame(x, y)

mod_flat <- withCallingHandlers(
  onls(y ~ b1 / (1 + exp(b2 - b3 * x))^(1 / b4), data = DAT_rich,
       start = list(b1 = 10, b2 = -1, b3 = 7, b4 = 9)),
  warning = function(w) {
    message("Warning: ", conditionMessage(w))
    invokeRestart("muffleWarning")
  })

Even with a sensible start (here read off the data: plateau near 750, inflection around $x = 7$) this data set has a peculiarity: its orthogonal optimum lies in the limit $b_4 \to 0$, where the Richards curve turns into the Gompertz curve $b_1 \exp(-\exp(c - b_3 x))$ with $c = b_2 - \log b_4$. Only this combination of $b_2$ and $b_4$ is identified, so the Richards fit slides along a flat valley and does not converge cleanly. Fitting the Gompertz form directly avoids the problem:

mod_gomp <- onls(y ~ b1 * exp(-exp(c - b3 * x)), data = DAT_rich,
                 start = list(b1 = 750, c = 2, b3 = 0.5))
summary(mod_gomp)
check_o(mod_gomp, plot = FALSE)

9. Diagnostics and further functions

Beyond check_o(), several functions distinguish between the classical vertical residuals (what an ordinary nls() fit would report) and the orthogonal residuals that onls() itself minimizes:

| uses vertical residuals | uses orthogonal residuals | |---|---| | deviance(), fitted(), residuals(), logLik() | deviance_o(), residuals_o(), logLik_o() |

residuals_o() is demonstrated in Section 3.1.

A low orthogonal residual sum of squares is not by itself evidence of a good fit: with unit precisions and a steep model, an orthogonal fit can lower its objective by shifting observations horizontally. It is therefore worth comparing the vertical residuals (residuals()) with those of the ordinary least-squares fit that onls() starts from, and choosing sigma_x/sigma_y to reflect the actual measurement errors -- in particular a small sigma_x for a predictor that is essentially error-free, such as a time index.

logLik_o() includes a precision-based normalizing correction, so AIC()/BIC() built from it are valid for comparing fits that used different weighting schemes -- as long as they were fit to the same data (the same response and number of observations). Here we fit the DNase model a second time, now with sigma_x/sigma_y supplied, and compare it against the original unweighted fit from Section 3:

mod_uni_w <- onls(density ~ Asym / (1 + exp((xmid - log(conc)) / scal)),
                  data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1),
                  sigma_x = 0.05, sigma_y = 0.1)
AIC(logLik_o(mod_uni))
AIC(logLik_o(mod_uni_w))

9.1 Confidence intervals: confint()

confint() computes bootstrap confidence intervals for all parameters. Unlike confint.nls(), which uses profile likelihoods, it refits the orthogonal model to nonparametric case resamples of the data, and it is therefore fully consistent with the criterion that onls() minimizes. In each replicate the rows of the data are resampled with replacement and the model is refitted with onls(); fits that fail to converge or violate the internal orthogonality checks are discarded, as are pathological solutions that deviate from the original estimate by more than twenty standard errors. After k successful fits, the limits are the empirical quantiles of the bootstrap distribution (level, default $0.95$). Further arguments are passed to update.onls().

Because k refits are needed, this takes a while. We use k = 100 here to keep the build time short -- in practice, take k >= 200:

set.seed(123)
confint(mod_uni, k = 100)

Further reading

See ?onls for the full mathematical details (construction of $Qyy_i$/$Qx_i$, the joint Levenberg-Marquardt algorithm, fixed parameters and bounds, and the approximate parameter covariance), ?check_o for the two orthogonality criteria in detail, and the reference list in ?onls for the underlying literature (Boggs, Byrd, Rogers & Schnabel's ODRPACK, York's weighted linear regression, and Daeron & Vermeesch's generalized least squares framing).

onls documentation built on Sept. 21, 2026, 9:06 a.m.