Orthogonal Nonlinear Least-Squares Regression with onls

Andrej-Nikolai Spiess

2026-09-21

library(onls)
#> Loading required package: minpack.lm
#> Loading required package: rgl

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:

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

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)
#> Nonlinear orthogonal regression model
#>   model: density ~ Asym/(1 + exp((xmid - log(conc))/scal))
#>    data: DNase1
#>  Asym  xmid  scal 
#> 3.303 2.352 1.315 
#>  vertical residual sum-of-squares: 0.1832
#>  orthogonal residual sum-of-squares: 0.1782
#>  PASSED: 16 out of 16 fitted points are orthogonal.
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 1e-10
summary(mod_uni)
#> 
#> Formula: density ~ Asym/(1 + exp((xmid - log(conc))/scal))
#> 
#> Parameters:
#>      Estimate Std. Error t value Pr(>|t|)    
#> Asym   3.3034     1.4820   2.229 0.044079 *  
#> xmid   2.3520     1.1587   2.030 0.063351 .  
#> scal   1.3150     0.3012   4.366 0.000764 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 0.1187 on 13 degrees of freedom
#> Residual standard error of orthogonal distances: 0.1171 on 13 degrees of freedom
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 1e-10

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)

Orthogonal vs vertical fit with foot-point segments

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)
#>           conc          x0    density         y0    alpha      df/dx Ortho
#> 1   0.04882812  0.02879479 0.01593503 0.03678759 89.99984 0.96071836  TRUE
#> 2   0.04882812  0.03010727 0.01833056 0.03804126 89.99986 0.94978583  TRUE
#> 3   0.19531250  0.17729284 0.11088889 0.14183707 89.99999 0.58225289  TRUE
#> 4   0.19531250  0.19170494 0.14378148 0.15012768 90.00000 0.56846042  TRUE
#> 5   0.39062500  0.37665149 0.21278786 0.24347577 89.99999 0.45534239  TRUE
#> 6   0.39062500  0.37074460 0.19735993 0.24077865 89.99999 0.45787592  TRUE
#> 7   0.78125000  0.77919269 0.39537608 0.40135468 89.99998 0.34411148  TRUE
#> 8   0.78125000  0.78111200 0.40161334 0.40201480 89.99998 0.34375231  TRUE
#> 9   1.56250000  1.56746128 0.64935297 0.62927803 89.99999 0.24713755  TRUE
#> 10  1.56250000  1.55369409 0.59040185 0.62586771 89.99999 0.24829245  TRUE
#> 11  3.12500000  3.16187573 1.17305050 0.94592294 90.00000 0.16235687  TRUE
#> 12  3.12500000  3.14090312 1.04002331 0.94251025 90.00000 0.16308705  TRUE
#> 13  6.25000000  6.24247076 1.25112651 1.32893758 89.99999 0.09676313  TRUE
#> 14  6.25000000  6.22426885 1.06191494 1.32717412 89.99999 0.09700388  TRUE
#> 15 12.50000000 12.50817620 1.92461305 1.76105267 89.99999 0.04998904  TRUE
#> 16 12.50000000 12.49708885 1.70231635 1.76049817 89.99999 0.05003562  TRUE

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)
#>            x         x0          y         y0        dist
#> 1 0.04882812 0.02879479 0.01593503 0.03678759 0.028916488
#> 2 0.04882812 0.03010727 0.01833056 0.03804126 0.027184228
#> 3 0.19531250 0.17729284 0.11088889 0.14183707 0.035811982
#> 4 0.19531250 0.19170494 0.14378148 0.15012768 0.007299918
#> 5 0.39062500 0.37665149 0.21278786 0.24347577 0.033719534
#> 6 0.39062500 0.37074460 0.19735993 0.24077865 0.047753701

## for unit precisions, the sum of squared distances is the minimized objective
all.equal(sum(FP$dist^2), deviance_o(mod_uni))
#> [1] "Attributes: < target is NULL, current is list >"

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)
#> [1] TRUE
all.equal(sum(d_o^2), deviance_o(mod_uni))
#> [1] "Attributes: < target is NULL, current is list >"

## side by side with the vertical residuals of the same fit
head(data.frame(vertical = residuals(mod_uni), orthogonal = d_o))
#>       vertical  orthogonal
#> 1 -0.038732828 0.028916488
#> 2 -0.036337299 0.027184228
#> 3 -0.041283634 0.035811982
#> 4 -0.008391046 0.007299918
#> 5 -0.037009708 0.033719534
#> 6 -0.052437638 0.047753701

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)
#> Warning in X - X0: longer object length is not a multiple of shorter object
#> length
#> Warning in Y - Y0: longer object length is not a multiple of shorter object
#> length
#> Warning in sqrt(dx^2 + dy^2) <= sqrt(.Machine$double.eps) * sc: longer object
#> length is not a multiple of shorter object length
residuals_o(mod_na)
#>          1          2          3          4          5          6          7 
#> 0.03139515 0.02973997         NA 0.01436914 0.04256587 0.05656476 0.01604731 
#>          8          9         10         11         12         13         14 
#> 0.01014890 0.01222736         NA 0.22537196 0.09402879 0.07872927 0.26707291 
#>         15         16 
#> 0.16397977 0.05803631 
#> attr(,"label")
#> [1] "Orthogonal residuals from orthogonal model"

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)

Noisy quadratic with orthogonal and vertical fit and foot-point segments

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)

Zoom into the left half of the quadratic fit, without the NLS curve

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
#> 
#> Formula: y ~ b0 + b1 * x
#> 
#> Parameters:
#>    Estimate Std. Error t value Pr(>|t|)    
#> b0  5.47991    0.29497  18.578 7.27e-08 ***
#> b1 -0.48053    0.05799  -8.287 3.38e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 0.3585 on 8 degrees of freedom
#> Residual standard error of orthogonal distances: 1.218 on 8 degrees of freedom
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 1e-10

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)
#>      x            x0   y       y0    rel_resid      df/dx Ortho
#> 1  0.0 -0.0002018208 5.9 5.480007 6.358921e-07 -0.4805333  TRUE
#> 2  0.9  0.8996951674 5.4 5.047576 6.358920e-07 -0.4805333  TRUE
#> 3  1.8  1.8008248014 4.4 4.614554 6.358529e-07 -0.4805333  TRUE
#> 4  2.6  2.5982286299 4.6 4.231374 6.358920e-07 -0.4805333  TRUE
#> 5  3.3  3.3185127438 3.5 3.885254 6.358724e-07 -0.4805333  TRUE
#> 6  4.4  4.3620157262 3.7 3.383816 6.358934e-07 -0.4805333  TRUE
#> 7  5.2  5.2799979510 2.8 2.942695 6.358734e-07 -0.4805333  TRUE
#> 8  6.1  5.8662161319 2.8 2.660997 6.358586e-07 -0.4805333  TRUE
#> 9  6.5  6.4159121599 2.4 2.396850 6.359340e-07 -0.4805333  TRUE
#> 10 7.4  8.2747003350 1.5 1.503641 6.356767e-07 -0.4805333  TRUE

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
#> 
#> Formula: y ~ a * exp(-b * x) + c
#> 
#> Parameters:
#>   Estimate Std. Error t value Pr(>|t|)    
#> a  1.98131    0.03347   59.20  < 2e-16 ***
#> b  0.28984    0.01468   19.74  < 2e-16 ***
#> c  0.48480    0.02986   16.23 1.87e-15 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 0.06119 on 27 degrees of freedom
#> Residual standard error of orthogonal distances: 0.05895 on 27 degrees of freedom
#> 
#> Number of iterations to convergence: 4 
#> Achieved convergence tolerance: 1e-10

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)
#> [1] TRUE

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
#> [1] 21.4455
#> attr(,"label")
#> [1] "Deviance (RSS) of orthogonal residuals from orthogonal model"
summary(mod_guide)      # 1264.65481 (1.03492) / -54.01838 (1.583992) / -0.08785 (6.33222E-3), as on page 48
#> 
#> Formula: y ~ b1 + b2 * (exp(b3 * x) - 1)^2
#> 
#> Parameters:
#>      Estimate Std. Error t value Pr(>|t|)    
#> b1  1.265e+03  1.035e+00 1221.97  < 2e-16 ***
#> b2 -5.402e+01  1.584e+00  -34.10 7.91e-11 ***
#> b3 -8.785e-02  6.332e-03  -13.87 2.22e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 2.3 on 9 degrees of freedom
#> Residual standard error of orthogonal distances: 1.544 on 9 degrees of freedom
#> 
#> Number of iterations to convergence: 8 
#> Achieved convergence tolerance: 1e-10

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)
#>        x        x0      y       y0    alpha         df/dx Ortho
#> 1    0.0  0.000000 1265.0 1264.655 90.00000  8.411821e-10  TRUE
#> 2    0.0  0.000000 1263.6 1264.655 90.00000  8.411821e-10  TRUE
#> 3    5.0  4.934914 1258.0 1257.970 89.99999 -2.164255e+00  TRUE
#> 4    7.0  6.732794 1254.0 1253.886 89.99999 -2.345577e+00  TRUE
#> 5    7.5  7.168638 1253.0 1252.860 89.99999 -2.362594e+00  TRUE
#> 6   10.0  8.693582 1249.8 1249.247 89.99999 -2.361733e+00  TRUE
#> 7   16.0 14.674734 1237.0 1236.300 90.00000 -1.894398e+00  TRUE
#> 8   26.0 27.458864 1218.0 1219.884 89.99996 -7.742995e-01  TRUE
#> 9   30.0 28.811955 1220.6 1218.891 89.99998 -6.951041e-01  TRUE
#> 10  34.0 34.771259 1213.8 1215.609 89.99996 -4.263092e-01  TRUE
#> 11  34.5 34.582416 1215.5 1215.690 89.99997 -4.330819e-01  TRUE
#> 12 100.0 99.998044 1212.0 1210.653 90.00000 -1.452283e-03  TRUE

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
#> [1] 15.26281
#> attr(,"label")
#> [1] "Deviance (RSS) of orthogonal residuals from orthogonal model"
summary(mod_676)      # 4.4879 (0.56876) / 7.1882 (0.69504) / 221.8383 (37.2313), as on page 363
#> 
#> Formula: y ~ b1 * 10^(b2 * x/(b3 + x))
#> 
#> Parameters:
#>    Estimate Std. Error t value Pr(>|t|)    
#> b1   4.4879     0.5688   7.891 7.44e-06 ***
#> b2   7.1882     0.6951  10.342 5.28e-07 ***
#> b3 221.8378    37.2323   5.958 9.48e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 14.45 on 11 degrees of freedom
#> Residual standard error of orthogonal distances: 1.178 on 11 degrees of freedom
#> 
#> Number of iterations to convergence: 11 
#> Achieved convergence tolerance: 1e-10

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
#> 
#> Formula: y ~ a + b * x
#> 
#> Parameters:
#>   Estimate Std. Error t value Pr(>|t|)   
#> a 13.70916    1.50991   9.079  0.01191 * 
#> b  0.85163    0.05438  15.660  0.00405 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 1.171 on 2 degrees of freedom
#> Residual standard error of orthogonal distances: 0.8916 on 2 degrees of freedom
#> 
#> Number of iterations to convergence: 5 
#> Achieved convergence tolerance: 1e-10

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
#> 
#> Formula: z ~ c * sqrt(1 + (x1/a)^2 + (x2/b)^2)
#> 
#> Parameters:
#>   Estimate Std. Error t value Pr(>|t|)    
#> a   3.0448     0.1458   20.89   <2e-16 ***
#> b   2.1083     0.0751   28.07   <2e-16 ***
#> c   4.1001     0.1077   38.07   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 0.3611 on 57 degrees of freedom
#> Residual standard error of orthogonal distances: 0.946 on 57 degrees of freedom
#> 
#> Number of iterations to convergence: 6 
#> Achieved convergence tolerance: 1e-10
check_o(mod_hyp, plot = FALSE)
#>             x1         x2       x0_x1      x0_x2         z        y0
#> 1   3.31637169 -4.1460487  3.43691527 -4.3239707 10.891614 10.437372
#> 2  -1.85703817 -0.6958970 -1.89347425 -0.7119703  5.142744  5.022811
#> 3   1.54289801 -1.1748560  1.60755020 -1.2330349  5.480296  5.219804
#> 4   1.94227335  1.2943736  1.96762909  1.3142428  5.598371  5.510261
#> 5  -0.57219586 -1.9894797 -0.59445663 -2.0809000  6.086398  5.816142
#> 6   2.01997134  3.2610579  2.00344632  3.2298035  7.889616  7.971200
#> 7  -1.05758904  4.0070166 -1.01602597  3.8235161  8.164071  8.600642
#> 8  -2.18626713  2.1342758 -2.21479987  2.1670283  6.698198  6.592807
#> 9   3.65207878 -2.0115378  3.68037959 -2.0298500  7.618809  7.546800
#> 10 -3.70560605 -2.0558336 -3.58587768 -1.9783386  7.104299  7.411357
#> 11  3.82479352  1.4275194  3.97372464  1.4931747  7.681311  7.339959
#> 12  4.55778304 -3.1723957  4.43433401 -3.0720589  9.064866  9.389213
#> 13  1.90558734  1.4039676  1.91407271  1.4113077  5.597177  5.566556
#> 14  0.16520093 -0.4000851  0.16527916 -0.4003074  4.181704  4.179250
#> 15 -0.16727403 -1.6835699 -0.17358522 -1.7585823  5.585369  5.344263
#> 16  3.25371695  2.1726698  3.23814695  2.1604823  7.269131  7.312762
#> 17 -4.13704362 -1.9696061 -4.23366411 -2.0237925  8.277878  8.049914
#> 18  3.73776762 -4.5677915  3.80061458 -4.6581599 11.412668 11.183203
#> 19 -3.62407938  2.0025521 -3.47168902  1.9044759  6.843349  7.237561
#> 20 -2.10571389 -3.0681032 -2.05671754 -2.9846853  7.401555  7.627013
#> 21 -4.41695680 -3.0848892 -4.36163026 -3.0396536  9.140841  9.287019
#> 22 -4.12937601 -2.4751312 -4.16891090 -2.5029787  8.586210  8.486348
#> 23  1.64506784 -0.1353757  1.70866225 -0.1415568  4.927094  4.709590
#> 24  3.84865476 -4.7929777  3.77090966 -4.6797835 10.912596 11.199100
#> 25  3.14630150  2.0517050  3.08287763  2.0033519  6.836771  7.015873
#> 26 -1.93068754 -0.3223832 -2.01284513 -0.3385973  5.210057  4.958901
#> 27  1.49950507  4.6446288  1.51670612  4.7072604 10.380373 10.236321
#> 28  3.18350800 -1.8993748  3.15946760 -1.8825693  6.885225  6.950853
#> 29 -4.82029199 -3.7111148 -4.88688768 -3.7714109 10.853170 10.672698
#> 30 -4.57868102  3.6523600 -4.65886260  3.7276259 10.649456 10.426784
#> 31 -0.81621090  0.4152604 -0.75270131  0.3778566  3.838105  4.286939
#> 32 -2.10583416  4.1751985 -2.08995261  4.1383050  9.371130  9.460334
#> 33 -2.75927995  1.2048784 -2.81184629  1.2318969  6.214308  6.073421
#> 34  4.98268250 -5.0163720  4.90025084 -4.9192876 12.066370 12.323607
#> 35  1.59818459  2.5762397  1.61342261  2.6051051  6.950510  6.869998
#> 36  1.21924456 -3.1075141  1.20068376 -3.0521604  7.251134  7.392944
#> 37 -4.39290044  2.4906694 -4.29119454  2.4232873  8.259108  8.509365
#> 38  4.60004654  0.1462656  4.47422505  0.1415941  7.038371  7.292851
#> 39 -4.92112837  4.6488972 -4.94992012  4.6808394 12.090856 12.004215
#> 40  1.53212445  0.8674578  1.55845304  0.8849985  5.020051  4.916976
#> 41  4.07630822 -2.2017604  4.13263443 -2.2375392  8.305951  8.167814
#> 42  4.10405827  1.9304122  4.07340010  1.9135159  7.721119  7.793908
#> 43  1.23796719 -2.7605990  1.27380053 -2.8548174  7.359861  7.111619
#> 44 -0.76817734  0.2928545 -0.77723172  0.2969125  4.332489  4.270754
#> 45  3.43843031 -2.4069100  3.37012449 -2.3510067  7.443967  7.636010
#> 46 -0.64477913  3.2811056 -0.65021511  3.3136063  7.767560  7.687807
#> 47 -0.81589312 -3.4918418 -0.81765795 -3.5007064  8.044568  8.023080
#> 48 -2.55273660 -4.2882912 -2.60766982 -4.3969604 10.376560 10.112227
#> 49  1.78811396 -2.5587037  1.74937764 -2.4939167  6.587493  6.773605
#> 50  1.94403860  3.1223372  1.87292386  2.9891805  7.191804  7.547401
#> 51 -4.46034199 -4.0395153 -4.45275961 -4.0314613 10.665185 10.687769
#> 52  2.59608285  1.2978642  2.60814161  1.3049426  5.999562  5.965338
#> 53 -0.56834841  3.1144261 -0.56534040  3.0951057  7.274187  7.322533
#> 54 -0.06358762 -4.1752716 -0.06376183 -4.1886981  9.150792  9.119873
#> 55  2.61195264 -0.9495082  2.52580506 -0.9129761  5.377598  5.615245
#> 56  2.38980202  4.0464821  2.32733267  3.9229465  8.903749  9.210520
#> 57 -4.55252729  3.4048911 -4.73175397  3.5632331 10.750311 10.267724
#> 58 -3.04758905 -3.3953457 -3.07039357 -3.4251916  8.928788  8.847250
#> 59 -0.29262549 -0.1182806 -0.28413021 -0.1142722  3.970875  4.123873
#> 60 -2.29364057 -2.5303912 -2.30402250 -2.5438390  7.174981  7.135086
#>    rel_resid_x1 rel_resid_x2    df/dx_x1    df/dx_x2 Ortho
#> 1  2.921382e-06 6.569744e-07  0.59709059 -1.56675662  TRUE
#> 2  2.855104e-06 5.656511e-07 -0.68355762 -0.53607430  TRUE
#> 3  3.179163e-06 4.220068e-07  0.55843538 -0.89336953  TRUE
#> 4  2.848862e-06 5.338269e-07  0.64749081  0.90201417  TRUE
#> 5  3.464880e-06 7.585572e-08 -0.18533087 -1.35308847  TRUE
#> 6  2.926646e-06 4.047274e-07  0.45573936  1.53236537  TRUE
#> 7  3.242477e-06 2.590067e-07 -0.21420840  1.68128760  TRUE
#> 8  2.838459e-06 5.454163e-07 -0.60915339  1.24309453  TRUE
#> 9  2.359256e-06 1.039031e-06  0.88428519 -1.01721108  TRUE
#> 10 2.417225e-06 1.091000e-06 -0.87732454 -1.00951510  TRUE
#> 11 2.230472e-06 1.455867e-06  0.98167269  0.76935530  TRUE
#> 12 2.430313e-06 1.027793e-06  0.85637065 -1.23740055  TRUE
#> 13 2.861667e-06 4.859987e-07  0.62349711  0.95883777  TRUE
#> 14 3.734622e-06 4.627480e-07  0.07171043 -0.36224734  TRUE
#> 15 3.562770e-06 5.030779e-08 -0.05889623 -1.24447119  TRUE
#> 16 2.497271e-06 8.657401e-07  0.80292997  1.11732424  TRUE
#> 17 2.254394e-06 1.252447e-06 -0.95364713 -0.95079016  TRUE
#> 18 2.764202e-06 6.479305e-07  0.61624045 -1.57528143  TRUE
#> 19 2.481013e-06 1.146893e-06 -0.86978349  0.99516071  TRUE
#> 20 2.946846e-06 4.353272e-07 -0.48897060 -1.47997385  TRUE
#> 21 2.400884e-06 9.801481e-07 -0.85159886 -1.23782063  TRUE
#> 22 2.345081e-06 1.061998e-06 -0.89076746 -1.11544018  TRUE
#> 23 3.043973e-06 5.214577e-07  0.65786338 -0.11367308  TRUE
#> 24 2.749206e-06 6.335051e-07  0.61055610 -1.58034754  TRUE
#> 25 2.531740e-06 8.687507e-07  0.79677764  1.07990483  TRUE
#> 26 2.845011e-06 7.961701e-07 -0.73601644 -0.25823077  TRUE
#> 27 3.042383e-06 3.172390e-07  0.26867058  1.73913899  TRUE
#> 28 2.467298e-06 8.994126e-07  0.82421091 -1.02428975  TRUE
#> 29 2.452834e-06 9.671785e-07 -0.83027258 -1.33640960  TRUE
#> 30 2.499006e-06 9.404083e-07 -0.81019972  1.35204737  TRUE
#> 31 5.429726e-06 2.517281e-06 -0.31837422  0.33334172  TRUE
#> 32 2.939108e-06 3.907178e-07 -0.40058301  1.65434459  TRUE
#> 33 2.487794e-06 9.628024e-07 -0.83950009  0.76709817  TRUE
#> 34 2.593326e-06 7.857198e-07  0.72101288 -1.50964322  TRUE
#> 35 3.017245e-06 3.357712e-07  0.42584747  1.43409713  TRUE
#> 36 3.079425e-06 2.544060e-07  0.29449219 -1.56134865  TRUE
#> 37 2.311537e-06 1.127226e-06 -0.91441562  1.07700505  TRUE
#> 38 1.856168e-06 1.658431e-06  1.11245647  0.07342734  TRUE
#> 39 2.553432e-06 8.226151e-07 -0.74769934  1.47468732  TRUE
#> 40 3.051064e-06 3.368799e-07  0.57472245  0.68069784  TRUE
#> 41 2.304124e-06 1.128201e-06  0.91745274 -1.03603578  TRUE
#> 42 2.226528e-06 1.158993e-06  0.94768581  0.92851039  TRUE
#> 43 3.194779e-06 2.738220e-07  0.32478465 -1.51816788  TRUE
#> 44 3.557312e-06 2.407024e-07 -0.32999589  0.26292621  TRUE
#> 45 2.515569e-06 8.850194e-07  0.80028009 -1.16438697  TRUE
#> 46 3.137681e-06 2.012251e-07 -0.15336178  1.63007800  TRUE
#> 47 3.098797e-06 2.249915e-07 -0.18479622 -1.65015411  TRUE
#> 48 2.954917e-06 4.734662e-07 -0.46759337 -1.64443114  TRUE
#> 49 2.994135e-06 3.739480e-07  0.46830224 -1.39242561  TRUE
#> 50 3.078362e-06 4.071671e-07  0.44997179  1.49783747  TRUE
#> 51 2.539326e-06 8.218989e-07 -0.75544834 -1.42654480  TRUE
#> 52 2.555190e-06 8.170503e-07  0.79279092  0.82730624  TRUE
#> 53 3.134717e-06 1.789249e-07 -0.13999460  1.59854231  TRUE
#> 54 3.102021e-06 2.223990e-07 -0.01267753 -1.73700166  TRUE
#> 55 2.565508e-06 9.285080e-07  0.81563094 -0.61489459  TRUE
#> 56 2.955491e-06 4.458629e-07  0.45818081  1.61078734  TRUE
#> 57 2.572818e-06 1.083207e-06 -0.83562332  1.31244172  TRUE
#> 58 2.738480e-06 6.287500e-07 -0.62928632 -1.46415314  TRUE
#> 59 3.961514e-06 6.726174e-07 -0.12493228 -0.10479610  TRUE
#> 60 2.824634e-06 5.262076e-07 -0.58553122 -1.34834291  TRUE

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
#> 
#> Formula: z ~ sqrt(r^2 - x1^2 - x2^2)
#> 
#> Parameters:
#>   Estimate Std. Error t value Pr(>|t|)    
#> r  5.99775    0.01768   339.3   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 0.1502 on 59 degrees of freedom
#> Residual standard error of orthogonal distances: 0.9079 on 59 degrees of freedom
#> 
#> Number of iterations to convergence: 5 
#> Achieved convergence tolerance: 1e-10
check_o(mod_dome, plot = FALSE)   # all orthogonal to the dome surface
#>             x1         x2      x0_x1       x0_x2        z       y0 rel_resid_x1
#> 1  -0.83711561  3.4219365 -0.8382999  3.42677761 4.835096 4.850514 2.743587e-10
#> 2   0.23188986 -1.2045283  0.2337108 -1.21398705 5.766069 5.868957 7.828490e-10
#> 3  -2.36755994  1.4573945 -2.3738515  1.46126741 5.278940 5.310609 1.930333e-10
#> 4   1.70344010 -1.6496203  1.7205199 -1.66616040 5.375975 5.498796 1.521644e-10
#> 5   3.92288043 -1.4911203  3.9194902 -1.48983159 4.296823 4.288477 2.366742e-10
#> 6   2.79359165  0.8216875  2.7915709  0.82109313 5.253156 5.244614 1.898113e-10
#> 7  -3.91908376 -0.5925722 -3.9094501 -0.59111562 4.534974 4.509969 2.962689e-10
#> 8   3.13655325 -2.5007474  3.1338681 -2.49860656 4.470538 4.461936 2.489055e-10
#> 9  -3.85673158 -1.1840677 -3.8003838 -1.16676817 4.640786 4.490966 4.032239e-09
#> 10 -2.84920792  0.7445927 -2.7856436  0.72798125 5.531623 5.261489 1.409376e-08
#> 11  3.86821263 -1.0140447  3.8553437 -1.01067113 4.515615 4.481953 3.192389e-10
#> 12 -3.33387788  0.9877011 -3.4313109  1.01656683 4.505575 4.813080 5.977463e-08
#> 13 -1.64410177 -3.3587705 -1.6266304 -3.32307784 4.834704 4.720621 1.157038e-09
#> 14 -0.14208978 -0.1385513 -0.1431707 -0.13960532 5.892589 5.994418 1.062170e-09
#> 15  2.24300737  1.7174160  2.2835033  1.74842277 5.053225 5.263237 3.836369e-09
#> 16  1.80074137 -1.0292353  1.7806943 -1.01777713 5.778926 5.636160 7.721389e-10
#> 17 -0.07866027  2.8015712 -0.0787115  2.80339609 5.293915 5.301680 1.250800e-09
#> 18  3.43629743  0.7847006  3.4751285  0.79356797 4.702292 4.823564 6.404091e-10
#> 19 -1.05023115  2.2680821 -1.0531043  2.27428689 5.415553 5.449003 1.135632e-10
#> 20  1.27881463 -0.5368445  1.2828458 -0.53853679 5.792902 5.834151 2.055676e-10
#> 21  1.75838535 -1.1849759  1.7617965 -1.18727463 5.584442 5.608876 1.590955e-10
#> 22 -1.30530662 -3.2733787 -1.3075791 -3.27907754 4.829846 4.848807 1.761964e-10
#> 23 -1.98316244 -2.2480725 -1.9808665 -2.24546988 5.210386 5.196833 2.909993e-10
#> 24  3.79645116 -0.0878831  3.8048323 -0.08807711 4.612597 4.635572 1.568031e-10
#> 25 -0.95708123 -1.2253370 -0.9569381 -1.22515381 5.794729 5.792780 2.530215e-10
#> 26 -0.80906547 -2.8830044 -0.8038899 -2.86456211 5.283229 5.207791 4.845086e-10
#> 27 -4.39410831 -1.2860121 -4.2873795 -1.25477608 4.226281 4.002118 7.036706e-08
#> 28 -3.42181133 -2.4042651 -3.4104205 -2.39626154 4.345182 4.312772 3.410953e-10
#> 29 -0.81044244  4.2284345 -0.8063347  4.20700256 4.246214 4.198094 6.484351e-10
#> 30  0.99272872  1.4806770  0.9818129  1.46439587 5.876176 5.732768 7.635409e-10
#> 31  1.64421123 -0.3584545  1.6233866 -0.35391453 5.929356 5.763020 9.375575e-10
#> 32  3.01638046 -2.0470794  2.9892979 -2.02869978 4.885236 4.787642 6.621229e-10
#> 33 -0.91609104 -2.3058334 -0.9175364 -2.30947138 5.439362 5.458709 1.266740e-10
#> 34  0.91539224 -3.5115324  0.9204483 -3.53092795 4.701248 4.760080 2.815216e-11
#> 35  2.47660429  0.2887072  2.4411483  0.28457397 5.649886 5.471092 1.645271e-09
#> 36 -1.93707592  0.4356447 -1.9468995  0.43785404 5.591837 5.656050 3.016842e-12
#> 37  0.26539203 -3.9742190  0.2583425 -3.86865240 4.856960 4.576005 6.630385e-08
#> 38  0.24795389  1.2587993  0.2440664  1.23906360 6.073420 5.863292 1.971302e-09
#> 39 -1.16519292  2.6417539 -1.1756826  2.66553627 5.137437 5.242683 2.682202e-10
#> 40  0.32913242  3.0327654  0.3349175  3.08607204 4.932508 5.131961 4.093186e-09
#> 41  2.25522411  2.6365744  2.2700172  2.65386895 4.804673 4.876171 2.615958e-11
#> 42 -2.31077860  1.3748423 -2.2936202  1.36463354 5.461640 5.371231 4.706038e-10
#> 43 -2.79048397  1.7157110 -2.7761076  1.70687175 5.093831 5.035162 3.763056e-10
#> 44 -2.62872648  3.2612720 -2.6694452  3.31178868 4.083255 4.228375 5.087465e-09
#> 45  1.73720194  2.4928685  1.7597852  2.52527534 4.999099 5.147735 9.025687e-10
#> 46  2.73035728  3.2209474  2.7264267  3.21631059 4.279394 4.265558 2.753358e-10
#> 47  0.51639161  4.1604868  0.5231234  4.21472386 4.112409 4.235030 2.043220e-09
#> 48 -3.44102149  0.6651114 -3.4992473  0.67636581 4.643381 4.823985 4.346589e-09
#> 49 -0.23401814  2.9257226 -0.2341668  2.92758071 5.222014 5.229481 1.783943e-10
#> 50  1.10667419 -1.4312399  1.0927226 -1.41319657 5.890025 5.725545 7.312351e-10
#> 51  4.10462829  1.4176944  4.1171958  1.42203515 4.094734 4.123051 1.119465e-10
#> 52 -2.27637961  0.8585273 -2.2612950  0.85283823 5.571675 5.489285 3.859923e-10
#> 53  0.32837178 -1.0236252  0.3344388 -1.04253774 5.656278 5.896974 2.539265e-09
#> 54  3.33502290  2.9289792  3.3091133  2.90622409 4.143172 4.071446 7.233621e-10
#> 55 -3.57781045 -1.4688406 -3.5373409 -1.45222617 4.739691 4.620746 1.342147e-09
#> 56  0.34316146  1.4846011  0.3430432  1.48408966 5.805605 5.801107 1.766254e-11
#> 57  2.29832612  2.3504390  2.2963338  2.34840150 5.028252 5.018456 2.521492e-10
#> 58  0.12137637 -4.3035416  0.1225793 -4.34619273 4.040196 4.131419 8.067484e-10
#> 59  2.73455602 -2.0536203  2.7580005 -2.07122688 4.812932 4.906780 1.337483e-10
#> 60 -2.04023287  1.9297524 -2.0653027  1.95346466 5.137004 5.281244 6.072608e-10
#>    rel_resid_x2    df/dx_x1    df/dx_x2 Ortho
#> 1  1.750621e-10  0.17282702 -0.70647720  TRUE
#> 2  5.479148e-11 -0.03982153  0.20684887  TRUE
#> 3  1.195966e-10  0.44700176 -0.27516006  TRUE
#> 4  1.854918e-10 -0.31289028  0.30300458  TRUE
#> 5  2.465476e-10 -0.91395850  0.34740341  TRUE
#> 6  2.311165e-10 -0.53227385 -0.15655930  TRUE
#> 7  4.908074e-11  0.86684633  0.13106867  TRUE
#> 8  2.656057e-10 -0.70235609  0.55998259  TRUE
#> 9  4.083306e-09  0.84622864  0.25980340  TRUE
#> 10 1.441766e-08  0.52944012 -0.13836030  TRUE
#> 11 2.795949e-10 -0.86019272  0.22549791  TRUE
#> 12 5.981678e-08  0.71291380 -0.21120923  TRUE
#> 13 1.065934e-09  0.34457972  0.70394925  TRUE
#> 14 1.241237e-10  0.02388400  0.02328922  TRUE
#> 15 3.829526e-09 -0.43385907 -0.33219530  TRUE
#> 16 7.801753e-10 -0.31594107  0.18057990  TRUE
#> 17 1.970961e-10  0.01484652 -0.52877503  TRUE
#> 18 4.806018e-10 -0.72044830 -0.16451901  TRUE
#> 19 1.643319e-10  0.19326552 -0.41737673  TRUE
#> 20 1.177110e-10 -0.21988559  0.09230765  TRUE
#> 21 8.654853e-11 -0.31410865  0.21167782  TRUE
#> 22 1.496159e-10  0.26967029  0.67626486  TRUE
#> 23 2.669343e-10  0.38116802  0.43208430  TRUE
#> 24 2.975913e-11 -0.82079019  0.01900027  TRUE
#> 25 6.395074e-11  0.16519496  0.21149668  TRUE
#> 26 3.550829e-10  0.15436295  0.55005324  TRUE
#> 27 7.035606e-08  1.07127750  0.31352797  TRUE
#> 28 3.826967e-10  0.79077230  0.55561983  TRUE
#> 29 4.186799e-10  0.19207162 -1.00212208  TRUE
#> 30 6.170088e-10 -0.17126334 -0.25544308  TRUE
#> 31 1.328156e-09 -0.28169026  0.06141130  TRUE
#> 32 6.969092e-10 -0.62437795  0.42373676  TRUE
#> 33 1.850433e-10  0.16808670  0.42308012  TRUE
#> 34 4.203274e-11 -0.19336825  0.74177913  TRUE
#> 35 1.635370e-09 -0.44619029 -0.05201411  TRUE
#> 36 1.876978e-10  0.34421539 -0.07741339  TRUE
#> 37 6.606752e-08 -0.05645589  0.84542128  TRUE
#> 38 1.573997e-09 -0.04162618 -0.21132559  TRUE
#> 39 1.123256e-10  0.22425208 -0.50842980  TRUE
#> 40 3.920749e-09 -0.06526112 -0.60134358  TRUE
#> 41 1.270378e-11 -0.46553277 -0.54425269  TRUE
#> 42 4.651252e-10  0.42701948 -0.25406347  TRUE
#> 43 3.011810e-10  0.55134428 -0.33899045  TRUE
#> 44 5.079793e-09  0.63131701 -0.78322962  TRUE
#> 45 9.044868e-10 -0.34185620 -0.49056045  TRUE
#> 46 2.601600e-10 -0.63917240 -0.75401879  TRUE
#> 47 2.284746e-09 -0.12352294 -0.99520509  TRUE
#> 48 4.511472e-09  0.72538511 -0.14020892  TRUE
#> 49 1.438813e-10  0.04477820 -0.55982238  TRUE
#> 50 7.631705e-10 -0.19085040  0.24682305  TRUE
#> 51 2.017697e-10 -0.99857980 -0.34489872  TRUE
#> 52 3.793675e-10  0.41194709 -0.15536417  TRUE
#> 53 2.298267e-09 -0.05671363  0.17679198  TRUE
#> 54 7.368962e-10 -0.81276127 -0.71380645  TRUE
#> 55 1.255914e-09  0.76553455  0.31428390  TRUE
#> 56 1.727594e-10 -0.05913410 -0.25582872  TRUE
#> 57 2.278856e-10 -0.45757776 -0.46795300  TRUE
#> 58 7.039467e-10 -0.02967003  1.05198551  TRUE
#> 59 1.248539e-10 -0.56207950  0.42211529  TRUE
#> 60 6.224540e-10  0.39106366 -0.36988721  TRUE
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
#> 
#> Formula: z ~ b0 + b1 * x1 + b2 * x2^2 + b3 * sqrt(x3) + b4 * log(x4 + 
#>     1)
#> 
#> Parameters:
#>    Estimate Std. Error t value Pr(>|t|)    
#> b0  1.07302    0.78335    1.37    0.176    
#> b1  0.82247    0.03223   25.52  < 2e-16 ***
#> b2  0.46701    0.01638   28.51  < 2e-16 ***
#> b3  2.32536    0.19975   11.64  < 2e-16 ***
#> b4  3.20383    0.28101   11.40 4.14e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 0.8931 on 55 degrees of freedom
#> Residual standard error of orthogonal distances: 1.126 on 55 degrees of freedom
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 1e-10
check_o(mod_mv, plot = FALSE)
#>             x1        x2       x3        x4      x0_x1     x0_x2    x0_x3
#> 1   5.61879200 5.0879067 5.696754  6.000746  5.5139512 4.8322092 5.634316
#> 2   0.55590250 1.2412896 9.845152  3.728927  0.6997911 1.3384957 9.909768
#> 3   6.71043421 3.3948488 5.495541  6.119715  6.6833892 3.3491327 5.479208
#> 4   9.68187179 3.6687397 7.242301  6.076300  9.6197380 3.5571857 7.209589
#> 5   4.98501015 1.9310335 4.442625 10.022502  5.0850008 2.0336671 4.509191
#> 6   9.17784184 4.1172199 6.582825  7.958369  9.2865649 4.3562686 6.642460
#> 7   6.56693204 2.4121412 7.545158  9.290169  6.6044786 2.4587354 7.564456
#> 8   3.50647092 4.8166631 9.765502  8.255345  3.4673873 4.7234862 9.747806
#> 9   3.51752218 0.4932512 3.335610  3.074175  3.6517989 0.5291100 3.437984
#> 10  2.07702062 4.3433237 3.634255  8.837310  1.9900066 4.1605997 3.569145
#> 11  5.64915975 0.3746755 2.487899  6.144391  5.6469107 0.3742507 2.485883
#> 12  5.40688005 3.1492629 4.198677  7.026766  5.4067274 3.1490203 4.198572
#> 13  1.42742273 1.9264733 7.215666  4.423846  1.5726828 2.0788882 7.291711
#> 14  6.04503928 2.1284702 9.843158  6.387454  6.0891390 2.1769242 9.863009
#> 15  6.60265818 2.8633080 6.265470  6.446699  6.5096518 2.7349247 6.212721
#> 16  6.07610882 0.6156038 8.009476  8.462834  5.9639551 0.5826237 7.953257
#> 17  3.76559716 2.5396780 4.712928  8.323285  3.8255473 2.6189214 4.751806
#> 18  0.41333479 0.8153768 6.248696  9.472516  0.4821666 0.8447231 6.287501
#> 19  1.03711901 3.0905850 7.062698  7.758364  0.9467932 2.9558309 7.014486
#> 20  1.65507783 3.6847121 6.931380  8.075119  1.6502251 3.6757092 6.928774
#> 21  2.03668064 2.3540597 7.819214  2.083140  2.2449565 2.6305898 7.923809
#> 22  0.73980792 3.2296605 9.757762  6.617533  0.7943709 3.3211209 9.782423
#> 23  7.58186389 0.4694483 7.545026  3.993201  7.6548648 0.4874068 7.582503
#> 24  6.03966681 2.7172389 2.628439  5.461921  6.1594625 2.8921049 2.730916
#> 25  7.39136825 1.7667416 3.682776  3.193242  7.4028439 1.7770342 3.691219
#> 26  3.24980207 0.7996270 6.450513  7.576476  3.1322836 0.7548537 6.384766
#> 27 -0.01251893 0.3363646 7.022145  4.027751 -0.0758849 0.3259403 6.988260
#> 28  7.73904717 1.3706710 5.752881  6.696195  7.7805446 1.3999932 5.777287
#> 29 -0.14474396 0.3676947 3.288390  3.098810 -0.2043757 0.3569514 3.241569
#> 30  1.71239868 3.8896428 9.824505  3.460625  1.8284535 4.1316551 9.876708
#> 31  4.88658002 0.2846455 6.243691  5.909736  4.6350546 0.2525804 6.099723
#> 32  2.98948421 4.3691548 3.907278  6.724598  2.9596760 4.3043961 3.885901
#> 33  3.77864768 0.2812746 8.183964  9.215723  3.8488815 0.2916118 8.218597
#> 34  5.65312230 2.0697625 8.477986  8.028960  5.7758350 2.2064181 8.537356
#> 35  0.11135357 0.7532173 9.602470  5.182275  0.1544243 0.7699551 9.622099
#> 36  4.04901777 1.3360559 7.846513  2.614834  4.0257806 1.3205680 7.834777
#> 37  6.72463411 4.6750540 4.353793  4.205896  6.7854245 4.8230349 4.394785
#> 38  8.71183531 2.0443263 6.238425  4.090302  8.7945922 2.1334379 6.285090
#> 39  5.51944100 2.0420687 6.601037  7.379839  5.7546694 2.3171737 6.729225
#> 40  9.60319258 0.5298407 5.713619  6.724472  9.6054891 0.5304556 5.714977
#> 41  8.42845695 0.3359622 6.200162  5.437532  8.3279961 0.3197494 6.142863
#> 42  5.43122008 4.2453708 8.312230  9.758304  5.2795711 3.9435319 8.237537
#> 43  7.59676577 3.6094921 6.675345  9.790735  7.6794511 3.7666869 6.720434
#> 44  9.70528901 0.5638968 5.938802  4.144041  9.4668653 0.5033277 5.798837
#> 45  5.41384481 3.2776243 2.802756  5.138898  5.2775329 3.0666414 2.685161
#> 46  9.29178534 3.5748901 9.916818  3.703802  9.1280235 3.3019696 9.843029
#> 47  3.66400029 2.1383352 8.876380  2.390644  3.5028512 1.9774953 8.799584
#> 48  4.16454901 3.2609965 8.736659  4.424804  4.0900533 3.1428280 8.700958
#> 49  6.07011463 0.7119269 8.955285  2.653649  6.0905729 0.7193547 8.964944
#> 50  0.99972528 0.9955491 2.195873  6.654212  0.9224809 0.9581923 2.120893
#> 51  6.85116330 2.3542185 7.034289  8.891751  6.9098614 2.4260941 7.065506
#> 52  8.17465242 4.7005051 6.430949  7.835172  8.1134259 4.5596031 6.396727
#> 53  7.89091778 1.1678032 6.447807  2.935066  7.8981881 1.1721042 6.451853
#> 54 10.10215160 5.0977765 1.417108  8.145280 10.1381714 5.1921696 1.459260
#> 55  2.61067139 3.0488042 8.692222  3.530999  2.6802192 3.1597171 8.725505
#> 56  8.00782718 1.5467286 1.815612  4.164790  8.2017851 1.7145758 2.009055
#> 57  0.76112534 1.4426118 6.310713  5.088857  0.5415697 1.2986978 6.185922
#> 58  5.70557378 0.9355302 4.027586  5.771934  5.7677730 0.9658514 4.071164
#> 59  9.04393154 1.3446890 7.877682  2.239446  9.0252986 1.3321608 7.868291
#> 60  4.36886049 1.9879074 2.722036  8.180615  4.3245599 1.9444311 2.683809
#>        x0_x4         z        y0 rel_resid_x1 rel_resid_x2 rel_resid_x3
#> 1   5.941916 27.886043 28.240127 2.825776e-07 3.805126e-07 3.730699e-07
#> 2   3.844622 15.346631 14.860667 2.825983e-07 9.508336e-07 3.980808e-07
#> 3   6.104887 23.441998 23.533338 2.825743e-07 4.977199e-07 3.867745e-07
#> 4   6.041930 27.181682 27.391530 2.825911e-07 5.052127e-07 3.898465e-07
#> 5  10.057727 20.161528 19.823824 2.826054e-07 8.966964e-07 4.220094e-07
#> 6   8.005399 30.975175 30.607978 2.825770e-07 5.806356e-07 3.912858e-07
#> 7   9.304363 23.323775 23.196967 2.825801e-07 7.495913e-07 4.037639e-07
#> 8   8.238866 28.595970 28.727970 2.826060e-07 4.463273e-07 3.873244e-07
#> 9   3.198750 13.569212 13.115711 2.825977e-07 1.072853e-06 4.122663e-07
#> 10  8.802733 22.206409 22.500286 2.826012e-07 4.029861e-07 3.596121e-07
#> 11  6.143165 15.740806 15.748402 2.826057e-07 6.680174e-07 4.320478e-07
#> 12  7.026692 21.587987 21.588503 2.828160e-07 5.250338e-07 3.871941e-07
#> 13  4.526238 16.631585 16.140990 2.825989e-07 8.037771e-07 3.956559e-07
#> 14  6.410635 22.163145 22.014204 2.825274e-07 7.439856e-07 3.995848e-07
#> 15  6.397725 21.813485 22.127601 2.825471e-07 4.956807e-07 3.891426e-07
#> 16  8.416438 19.500265 19.879048 2.825543e-07 6.987495e-07 4.044976e-07
#> 17  8.348266 19.855123 19.652650 2.825991e-07 6.596048e-07 4.003598e-07
#> 18  9.498057 15.398929 15.166459 2.825620e-07 9.256119e-07 4.171764e-07
#> 19  7.718005 18.723114 19.028177 2.825986e-07 4.582591e-07 3.849871e-07
#> 20  8.073035 21.909928 21.926317 2.824829e-07 4.856522e-07 3.870967e-07
#> 21  2.326998 17.251493 16.548070 2.825993e-07 7.934770e-07 3.847029e-07
#> 22  6.645334 20.851561 20.667282 2.825971e-07 5.416659e-07 3.915765e-07
#> 23  4.049517 19.317549 19.070999 2.826124e-07 9.623908e-07 4.101007e-07
#> 24  5.533346 20.305866 19.901274 2.825584e-07 6.654566e-07 3.883116e-07
#> 25  3.203875 17.743470 17.704712 2.825849e-07 5.423904e-07 3.864389e-07
#> 26  7.522764 16.259156 16.656058 2.825977e-07 5.425398e-07 3.964176e-07
#> 27  3.978168 12.135724 12.349734 2.825969e-07 4.116410e-07 3.837416e-07
#> 28  6.717142 20.663850 20.523698 2.825611e-07 8.516031e-07 4.163417e-07
#> 29  3.041332  9.424084  9.625482 2.825602e-07 2.749800e-07 3.361478e-07
#> 30  3.559770 23.109956 22.717998 2.825963e-07 4.736745e-07 3.876902e-07
#> 31  5.764902 15.933523 16.783014 2.825974e-07 8.845502e-07 3.773844e-07
#> 32  6.709537 23.186823 23.287496 2.825382e-07 4.085010e-07 3.650501e-07
#> 33  9.242434 18.635737 18.398532 2.825961e-07 1.112366e-06 4.187530e-07
#> 34  8.081595 22.374308 21.959863 2.825881e-07 9.263420e-07 4.045629e-07
#> 35  5.209295 14.685855 14.540390 2.825970e-07 7.271107e-07 3.989130e-07
#> 36  2.589617 15.723523 15.802003 2.825966e-07 4.665371e-07 3.877718e-07
#> 37  4.250993 27.910620 27.705309 2.825682e-07 4.301747e-07 3.716319e-07
#> 38  4.152863 21.793969 21.514469 2.826407e-07 7.582907e-07 4.006867e-07
#> 39  7.487795 21.991983 21.197532 2.826062e-07 1.277554e-06 4.014027e-07
#> 40  6.725630 21.221799 21.214043 2.825101e-07 8.833802e-07 4.254355e-07
#> 41  5.376158 19.329673 19.668965 2.826005e-07 5.765755e-07 4.023507e-07
#> 42  9.703111 26.434610 26.946783 2.825682e-07 4.625356e-07 3.858907e-07
#> 43  9.820502 27.952224 27.672966 2.825722e-07 6.320873e-07 3.957156e-07
#> 44  3.956667 18.900434 19.705676 2.825956e-07 4.373876e-07 3.791003e-07
#> 45  5.051148 18.923265 19.383639 2.825977e-07 3.484288e-07 3.398694e-07
#> 46  3.564032 25.278818 25.831901 2.825786e-07 4.036844e-07 3.857016e-07
#> 47  2.194115 15.854592 16.398851 2.825979e-07 2.493508e-07 3.783626e-07
#> 48  4.370773 21.042927 21.294526 2.825994e-07 4.302523e-07 3.859530e-07
#> 49  2.675332 17.525823 17.456728 2.825921e-07 6.469428e-07 3.961188e-07
#> 50  6.614697 11.890151 12.151032 2.825973e-07 4.081734e-07 3.577754e-07
#> 51  8.914812 23.233925 23.035681 2.826493e-07 7.785892e-07 4.053600e-07
#> 52  7.808095 30.100086 30.306870 2.826601e-07 4.567782e-07 3.839847e-07
#> 53  2.942249 18.536580 18.512026 2.825776e-07 6.114475e-07 3.963728e-07
#> 54  8.160596 32.028120 31.906469 2.826292e-07 4.544736e-07 3.105425e-07
#> 55  3.590022 19.925966 19.691078 2.825268e-07 5.107019e-07 3.885035e-07
#> 56  4.307153 18.490069 17.835004 2.825824e-07 1.317072e-06 4.341193e-07
#> 57  4.944996 13.059103 13.800621 2.825976e-07 1.962620e-08 3.715888e-07
#> 58  5.807526 17.299519 17.089450 2.825994e-07 8.542884e-07 4.237461e-07
#> 59  2.216883 19.528051 19.590981 2.825565e-07 5.363079e-07 3.923241e-07
#> 60  8.161779 17.151990 17.301609 2.825978e-07 5.496267e-07 4.002574e-07
#>    rel_resid_x4  df/dx_x1  df/dx_x2  df/dx_x3  df/dx_x4 Ortho
#> 1  9.556521e-07 0.8224703 4.5133518 0.4898238 0.4615200  TRUE
#> 2  9.598100e-07 0.8224703 1.2501739 0.3693420 0.6613174  TRUE
#> 3  9.613770e-07 0.8224703 3.1281373 0.4967085 0.4509337  TRUE
#> 4  9.701475e-07 0.8224703 3.3224618 0.4330173 0.4549652  TRUE
#> 5  9.696831e-07 0.8224703 1.8994738 0.5475341 0.2897370  TRUE
#> 6  9.508743e-07 0.8224703 4.0688165 0.4511241 0.3557680  TRUE
#> 7  9.732895e-07 0.8224703 2.2964937 0.4227383 0.3109200  TRUE
#> 8  9.581042e-07 0.8224703 4.4118029 0.3723977 0.3467777  TRUE
#> 9  9.376461e-07 0.8224703 0.4941963 0.6270592 0.7630444  TRUE
#> 10 9.606661e-07 0.8224703 3.8860590 0.6154296 0.3268306  TRUE
#> 11 9.817922e-07 0.8224703 0.3495554 0.7374289 0.4485173  TRUE
#> 12 9.606388e-07 0.8224703 2.9412295 0.5674266 0.3991474  TRUE
#> 13 9.387550e-07 0.8224703 1.9417110 0.4305719 0.5797493  TRUE
#> 14 9.786320e-07 0.8224703 2.0332780 0.3702165 0.4323291  TRUE
#> 15 9.755419e-07 0.8224703 2.5544583 0.4664655 0.4330835  TRUE
#> 16 9.933642e-07 0.8224703 0.5441788 0.4122759 0.3402383  TRUE
#> 17 9.632010e-07 0.8224703 2.4461096 0.5333731 0.3427195  TRUE
#> 18 9.782077e-07 0.8224703 0.7889833 0.4636833 0.3051834  TRUE
#> 19 9.656104e-07 0.8224703 2.7607878 0.4389980 0.3674961  TRUE
#> 20 9.582904e-07 0.8224703 3.4331644 0.4417050 0.3531159  TRUE
#> 21 8.828432e-07 0.8224703 2.4570081 0.4130413 0.9629802  TRUE
#> 22 9.526360e-07 0.8224703 3.1019739 0.3717382 0.4190573  TRUE
#> 23 1.014783e-06 0.8224703 0.4552448 0.4222349 0.6344831  TRUE
#> 24 9.356664e-07 0.8224703 2.7012669 0.7035683 0.4903816  TRUE
#> 25 9.440937e-07 0.8224703 1.6597751 0.6051674 0.7621142  TRUE
#> 26 9.849380e-07 0.8224703 0.7050440 0.4601378 0.3759148  TRUE
#> 27 9.533417e-07 0.8224703 0.3044329 0.4398210 0.6435768  TRUE
#> 28 9.883730e-07 0.8224703 1.3076135 0.4837249 0.4151580  TRUE
#> 29 8.922822e-07 0.8224703 0.3333976 0.6457774 0.7927666  TRUE
#> 30 8.906221e-07 0.8224703 3.8590244 0.3699596 0.7026304  TRUE
#> 31 9.954380e-07 0.8224703 0.2359137 0.4707663 0.4735964  TRUE
#> 32 9.497790e-07 0.8224703 4.0203668 0.5898133 0.4155675  TRUE
#> 33 9.906334e-07 0.8224703 0.2723695 0.4055661 0.3128000  TRUE
#> 34 9.710424e-07 0.8224703 2.0608257 0.3979227 0.3527831  TRUE
#> 35 9.805302e-07 0.8224703 0.7191490 0.3748224 0.5159737  TRUE
#> 36 9.356933e-07 0.8224703 1.2334292 0.4153815 0.8925277  TRUE
#> 37 9.074252e-07 0.8224703 4.5047828 0.5546150 0.6101385  TRUE
#> 38 9.634527e-07 0.8224703 1.9926612 0.4637722 0.6217578  TRUE
#> 39 9.484489e-07 0.8224703 2.1642730 0.4482063 0.3774635  TRUE
#> 40 1.005449e-06 0.8224703 0.4954530 0.4863548 0.4147019  TRUE
#> 41 1.010050e-06 0.8224703 0.2986505 0.4691104 0.5024708  TRUE
#> 42 9.695104e-07 0.8224703 3.6833146 0.4050996 0.2993366  TRUE
#> 43 9.604552e-07 0.8224703 3.5181389 0.4484993 0.2960891  TRUE
#> 44 1.042918e-06 0.8224703 0.4701152 0.4828253 0.6463684  TRUE
#> 45 9.666573e-07 0.8224703 2.8642864 0.7095374 0.5294587  TRUE
#> 46 9.948756e-07 0.8224703 3.0840863 0.3705920 0.7019743  TRUE
#> 47 9.425678e-07 0.8224703 1.8470086 0.3919488 1.0030424  TRUE
#> 48 9.559271e-07 0.8224703 2.9354458 0.3941640 0.5965311  TRUE
#> 49 9.846201e-07 0.8224703 0.6718875 0.3883172 0.8717125  TRUE
#> 50 9.656037e-07 0.8224703 0.8949652 0.7983643 0.4207433  TRUE
#> 51 9.723799e-07 0.8224703 2.2660062 0.4374101 0.3231360  TRUE
#> 52 9.620202e-07 0.8224703 4.2587338 0.4597074 0.3637373  TRUE
#> 53 9.802542e-07 0.8224703 1.0947619 0.4577393 0.8126916  TRUE
#> 54 9.501369e-07 0.8224703 4.8495599 0.9624854 0.3497407  TRUE
#> 55 9.115716e-07 0.8224703 2.9512205 0.3936091 0.6979995  TRUE
#> 56 9.225429e-07 0.8224703 1.6014381 0.8202846 0.6036820  TRUE
#> 57 9.788368e-07 0.8224703 1.2130021 0.4674748 0.5389126  TRUE
#> 58 9.827643e-07 0.8224703 0.9021189 0.5762370 0.4706310  TRUE
#> 59 9.799560e-07 0.8224703 1.2442570 0.4144959 0.9959434  TRUE
#> 60 9.714962e-07 0.8224703 1.8161261 0.7097161 0.3496955  TRUE

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)

Grid of four partial-dependence panels

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")

Single full-size panel for one predictor

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
#>       gen_TLS       onls     abs_diff
#> b0  3.2925307  3.2925307 1.146137e-08
#> b1  1.5044330  1.5044330 1.035817e-09
#> b2 -0.8344564 -0.8344564 1.412917e-09

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)
#>         b0         b1         b2 
#>  3.3853065  1.5012131 -0.8494335

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
#> Nonlinear orthogonal regression model
#>   model: y ~ a + b * x
#>    data: DAT_dem
#>      a      b 
#> -1.909  1.208 
#>  vertical residual sum-of-squares: 9.224
#>  orthogonal residual sum-of-squares: 3.751
#>  PASSED: 16 out of 16 fitted points are orthogonal.
#> 
#> Number of iterations to convergence: 9 
#> Achieved convergence tolerance: 1e-10

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
#>    TLS_closed_form       onls     abs_diff
#> b0        3.681744  3.6817440 4.529555e-07
#> b1        1.468923  1.4689227 2.333856e-09
#> b2       -0.923844 -0.9238439 1.007331e-07

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
#>       b1       b2 
#> 1.437562 0.900000
deviance_o(mod_bnd)    # 0.1919, lower than the 0.2674 of the original ODRPACK
#> [1] 0.1918817
#> attr(,"label")
#> [1] "Deviance (RSS) of orthogonal residuals from orthogonal model"

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)
#> Nonlinear orthogonal regression model
#>   model: density ~ Asym/(1 + exp((xmid - log(conc))/scal))
#>    data: DNase1
#>  Asym  xmid  scal 
#> 3.000 2.106 1.253 
#>  vertical residual sum-of-squares: 0.183
#>  orthogonal residual sum-of-squares: 0.1788
#>  PASSED: 16 out of 16 fitted points are orthogonal.
#> 
#> Number of iterations to convergence: 5 
#> Achieved convergence tolerance: 1e-10

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)
#>            b1            b2            b3 
#> 1264.65481020  -54.01838894   -0.08784984
mod_ctrl$convInfo$isConv
#> [1] TRUE
mod_ctrl$convInfo$finIter   # total Levenberg-Marquardt iterations
#> [1] 9

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")
  })
#> Warning: NLS warm-start failed (singular gradient matrix at initial parameter estimates); using raw start values instead.
#> Warning: Parameter(s) b2, b3, b4 have no measurable influence on the fitted model: the solution is a degenerate stationary point (e.g. a model collapsed to a constant). Results are not meaningful -- try different `start` values.
#> Warning: Covariance matrix singular; SEs unavailable.

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)
#> 
#> Formula: y ~ b1 * exp(-exp(c - b3 * x))
#> 
#> Parameters:
#>     Estimate Std. Error t value Pr(>|t|)    
#> b1 837.36088  101.92414   8.216 2.86e-06 ***
#> c    1.71631    0.27825   6.168 4.81e-05 ***
#> b3   0.27560    0.06802   4.052  0.00161 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error of vertical distances: 56.16 on 12 degrees of freedom
#> Residual standard error of orthogonal distances: 1.397 on 12 degrees of freedom
#> 
#> Number of iterations to convergence: 35 
#> Achieved convergence tolerance: 1e-10
check_o(mod_gomp, plot = FALSE)
#>     x        x0      y        y0    alpha    df/dx Ortho
#> 1   1  1.239877  16.08  16.06630 89.99999 17.50543  TRUE
#> 2   2  1.997074  33.83  33.83010 89.99999 29.91795  TRUE
#> 3   3  2.840179  65.80  65.80346 90.00000 46.12819  TRUE
#> 4   4  3.444444  97.20  97.20963 90.00000 57.69030  TRUE
#> 5   5  4.817179 191.55 191.55235 90.00000 77.87157  TRUE
#> 6   6  6.441535 326.20 326.19479 90.00000 84.75203  TRUE
#> 7   7  7.165805 386.87 386.86799 90.00000 82.32816  TRUE
#> 8   8  8.925557 520.53 520.51643 90.00000 68.20185  TRUE
#> 9   9 10.035782 590.03 590.01181 90.00000 56.92987  TRUE
#> 10 10 11.252389 651.92 651.89216 90.00000 44.98215  TRUE
#> 11 11 13.252254 724.93 724.85186 89.99998 28.82377  TRUE
#> 12 12 12.453406 699.56 699.54692 89.99999 34.66826  TRUE
#> 13 13 12.185728 689.96 689.98212 89.99999 36.81216  TRUE
#> 14 14 10.945056 637.56 637.62380 90.00000 47.88654  TRUE
#> 15 15 13.003734 717.41 717.47534 89.99998 30.55300  TRUE

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))
#> [1] -18.55403
AIC(logLik_o(mod_uni_w))
#> [1] -114.1184

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)
#> ......
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> .
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ..10.........20..
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> .
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ......30..
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> .......40.........50
#> ...
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ..
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ...
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> .60..
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ....
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ...
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> 70..
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> .
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ...
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ...80....
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> .....90.
#> Warning in value[[3L]](cond): Covariance matrix singular; SEs unavailable.
#> ........100
#>           2.5%     97.5%
#> Asym 1.2646580 12.938978
#> xmid 0.2811913  5.527804
#> scal 0.7237133  1.688956

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).