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
An ordinary (vertical) nonlinear fit via Levenberg-Marquardt
(minpack.lm::nlsLM()) gives warm-start parameter
values.
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 Levenberg-Marquardt, 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.
Two algorithms compute the Levenberg-Marquardt steps
(algorithm); the default "auto" uses
"structured" for more than 100 observations or when
parameter bounds are given, and "dense" otherwise.
"structured" eliminates the \(n\) small \(p
\times p\) blocks of the foot-point corrections from the normal
equations, as ODRPACK does, and chooses the damping by MINPACK’s
trust-region rule; its cost grows linearly with \(n\). "dense" passes the full
joint problem to minpack.lm::nls.lm(), at a cost growing
with \(n^3\); it is needed only for the
foot-point bounds extend/window (used
automatically then). Both stop when the convergence tolerances
ftol and ptol (default 1e-10) are
met or after control$outer_max iterations (default
5000).
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).
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. With precisions, the
same holds after each axis is rescaled by its precision (“whitened”
coordinates \(L_i x\) and \(Qyy_i^{1/2} y\), with \(L_i^T L_i = Qx_i\)), where the weighted
distance becomes an ordinary Euclidean one. check_o()
therefore computes, for every observation, the angle \(\alpha_i\) between the residual vector and
the tangent space of the model surface in these coordinates: \(90^\circ\) at a stationary foot point, for
any number of predictors and any form of sigma_x, including
correlated predictor errors. For one predictor and unit precisions this
is the classical tangent angle. The tolerance is tol_deg
(default \(0.05^\circ\)); observations
that lie on the fitted surface have no angle and count as orthogonal.
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
#> Convergence tolerance (ftol): 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
#> Convergence tolerance (ftol): 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.
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.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
#>
#> Formula: y ~ b0 + b1 * x
#>
#> Parameters (variances known: z values):
#> Estimate Std. Error z value Pr(>|z|)
#> b0 5.47991 0.29497 18.578 <2e-16 ***
#> b1 -0.48053 0.05799 -8.287 <2e-16 ***
#> ---
#> 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
#> Convergence tolerance (ftol): 1e-10
For this weighted fit, check_o() measures the angle
after rescaling each axis by its precision (in these whitened
coordinates the weighted distance is Euclidean again), so 90 degrees is
still the criterion:
check_o(mod_w, plot = FALSE)
#> x x0 y y0 alpha df/dx Ortho
#> 1 0.0 -0.0002018208 5.9 5.480007 90.00000 -0.4805333 TRUE
#> 2 0.9 0.8996951674 5.4 5.047576 90.00000 -0.4805333 TRUE
#> 3 1.8 1.8008248014 4.4 4.614554 90.00000 -0.4805333 TRUE
#> 4 2.6 2.5982286299 4.6 4.231374 90.00000 -0.4805333 TRUE
#> 5 3.3 3.3185127438 3.5 3.885254 89.99999 -0.4805333 TRUE
#> 6 4.4 4.3620157262 3.7 3.383816 89.99999 -0.4805333 TRUE
#> 7 5.2 5.2799979510 2.8 2.942695 89.99999 -0.4805333 TRUE
#> 8 6.1 5.8662161319 2.8 2.660997 89.99998 -0.4805333 TRUE
#> 9 6.5 6.4159121599 2.4 2.396850 89.99999 -0.4805333 TRUE
#> 10 7.4 8.2747003350 1.5 1.503641 90.00000 -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
#> Convergence tolerance (ftol): 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
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 (variances known: z values):
#> Estimate Std. Error z value Pr(>|z|)
#> 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
#> Convergence tolerance (ftol): 1e-10
check_o(mod_hyp, plot = FALSE)
#> x1 x2 x0_x1 x0_x2 z y0 alpha
#> 1 3.31637169 -4.1460487 3.43691527 -4.3239707 10.891614 10.437372 89.99995
#> 2 -1.85703817 -0.6958970 -1.89347425 -0.7119703 5.142744 5.022811 89.99994
#> 3 1.54289801 -1.1748560 1.60755020 -1.2330349 5.480296 5.219804 89.99994
#> 4 1.94227335 1.2943736 1.96762909 1.3142428 5.598371 5.510261 89.99994
#> 5 -0.57219586 -1.9894797 -0.59445663 -2.0809000 6.086398 5.816142 89.99998
#> 6 2.01997134 3.2610579 2.00344632 3.2298035 7.889616 7.971200 89.99996
#> 7 -1.05758904 4.0070166 -1.01602597 3.8235161 8.164071 8.600642 89.99998
#> 8 -2.18626713 2.1342758 -2.21479987 2.1670283 6.698198 6.592807 89.99995
#> 9 3.65207878 -2.0115378 3.68037959 -2.0298500 7.618809 7.546800 89.99994
#> 10 -3.70560605 -2.0558336 -3.58587768 -1.9783386 7.104299 7.411357 89.99993
#> 11 3.82479352 1.4275194 3.97372464 1.4931747 7.681311 7.339959 89.99993
#> 12 4.55778304 -3.1723957 4.43433401 -3.0720589 9.064866 9.389213 89.99993
#> 13 1.90558734 1.4039676 1.91407271 1.4113077 5.597177 5.566556 89.99994
#> 14 0.16520093 -0.4000851 0.16527916 -0.4003074 4.181704 4.179250 89.99999
#> 15 -0.16727403 -1.6835699 -0.17358522 -1.7585823 5.585369 5.344263 89.99999
#> 16 3.25371695 2.1726698 3.23814695 2.1604823 7.269131 7.312762 89.99994
#> 17 -4.13704362 -1.9696061 -4.23366411 -2.0237925 8.277878 8.049914 89.99993
#> 18 3.73776762 -4.5677915 3.80061458 -4.6581599 11.412668 11.183203 89.99995
#> 19 -3.62407938 2.0025521 -3.47168902 1.9044759 6.843349 7.237561 89.99993
#> 20 -2.10571389 -3.0681032 -2.05671754 -2.9846853 7.401555 7.627013 89.99995
#> 21 -4.41695680 -3.0848892 -4.36163026 -3.0396536 9.140841 9.287019 89.99994
#> 22 -4.12937601 -2.4751312 -4.16891090 -2.5029787 8.586210 8.486348 89.99993
#> 23 1.64506784 -0.1353757 1.70866225 -0.1415568 4.927094 4.709590 89.99994
#> 24 3.84865476 -4.7929777 3.77090966 -4.6797835 10.912596 11.199100 89.99995
#> 25 3.14630150 2.0517050 3.08287763 2.0033519 6.836771 7.015873 89.99994
#> 26 -1.93068754 -0.3223832 -2.01284513 -0.3385973 5.210057 4.958901 89.99994
#> 27 1.49950507 4.6446288 1.51670612 4.7072604 10.380373 10.236321 89.99997
#> 28 3.18350800 -1.8993748 3.15946760 -1.8825693 6.885225 6.950853 89.99994
#> 29 -4.82029199 -3.7111148 -4.88688768 -3.7714109 10.853170 10.672698 89.99994
#> 30 -4.57868102 3.6523600 -4.65886260 3.7276259 10.649456 10.426784 89.99994
#> 31 -0.81621090 0.4152604 -0.75270131 0.3778566 3.838105 4.286939 89.99993
#> 32 -2.10583416 4.1751985 -2.08995261 4.1383050 9.371130 9.460334 89.99996
#> 33 -2.75927995 1.2048784 -2.81184629 1.2318969 6.214308 6.073421 89.99994
#> 34 4.98268250 -5.0163720 4.90025084 -4.9192876 12.066370 12.323607 89.99994
#> 35 1.59818459 2.5762397 1.61342261 2.6051051 6.950510 6.869998 89.99996
#> 36 1.21924456 -3.1075141 1.20068376 -3.0521604 7.251134 7.392944 89.99997
#> 37 -4.39290044 2.4906694 -4.29119454 2.4232873 8.259108 8.509365 89.99993
#> 38 4.60004654 0.1462656 4.47422505 0.1415941 7.038371 7.292851 89.99995
#> 39 -4.92112837 4.6488972 -4.94992012 4.6808394 12.090856 12.004215 89.99994
#> 40 1.53212445 0.8674578 1.55845304 0.8849985 5.020051 4.916976 89.99994
#> 41 4.07630822 -2.2017604 4.13263443 -2.2375392 8.305951 8.167814 89.99993
#> 42 4.10405827 1.9304122 4.07340010 1.9135159 7.721119 7.793908 89.99994
#> 43 1.23796719 -2.7605990 1.27380053 -2.8548174 7.359861 7.111619 89.99997
#> 44 -0.76817734 0.2928545 -0.77723172 0.2969125 4.332489 4.270754 89.99996
#> 45 3.43843031 -2.4069100 3.37012449 -2.3510067 7.443967 7.636010 89.99994
#> 46 -0.64477913 3.2811056 -0.65021511 3.3136063 7.767560 7.687807 89.99998
#> 47 -0.81589312 -3.4918418 -0.81765795 -3.5007064 8.044568 8.023080 89.99998
#> 48 -2.55273660 -4.2882912 -2.60766982 -4.3969604 10.376560 10.112227 89.99996
#> 49 1.78811396 -2.5587037 1.74937764 -2.4939167 6.587493 6.773605 89.99996
#> 50 1.94403860 3.1223372 1.87292386 2.9891805 7.191804 7.547401 89.99996
#> 51 -4.46034199 -4.0395153 -4.45275961 -4.0314613 10.665185 10.687769 89.99994
#> 52 2.59608285 1.2978642 2.60814161 1.3049426 5.999562 5.965338 89.99994
#> 53 -0.56834841 3.1144261 -0.56534040 3.0951057 7.274187 7.322533 89.99999
#> 54 -0.06358762 -4.1752716 -0.06376183 -4.1886981 9.150792 9.119873 89.99999
#> 55 2.61195264 -0.9495082 2.52580506 -0.9129761 5.377598 5.615245 89.99994
#> 56 2.38980202 4.0464821 2.32733267 3.9229465 8.903749 9.210520 89.99996
#> 57 -4.55252729 3.4048911 -4.73175397 3.5632331 10.750311 10.267724 89.99993
#> 58 -3.04758905 -3.3953457 -3.07039357 -3.4251916 8.928788 8.847250 89.99995
#> 59 -0.29262549 -0.1182806 -0.28413021 -0.1142722 3.970875 4.123873 89.99998
#> 60 -2.29364057 -2.5303912 -2.30402250 -2.5438390 7.174981 7.135086 89.99995
#> df/dx_x1 df/dx_x2 Ortho
#> 1 0.59709059 -1.56675662 TRUE
#> 2 -0.68355762 -0.53607430 TRUE
#> 3 0.55843538 -0.89336953 TRUE
#> 4 0.64749081 0.90201417 TRUE
#> 5 -0.18533087 -1.35308847 TRUE
#> 6 0.45573936 1.53236537 TRUE
#> 7 -0.21420840 1.68128760 TRUE
#> 8 -0.60915339 1.24309453 TRUE
#> 9 0.88428519 -1.01721108 TRUE
#> 10 -0.87732454 -1.00951510 TRUE
#> 11 0.98167269 0.76935530 TRUE
#> 12 0.85637065 -1.23740055 TRUE
#> 13 0.62349711 0.95883777 TRUE
#> 14 0.07171043 -0.36224734 TRUE
#> 15 -0.05889623 -1.24447119 TRUE
#> 16 0.80292997 1.11732424 TRUE
#> 17 -0.95364713 -0.95079016 TRUE
#> 18 0.61624045 -1.57528143 TRUE
#> 19 -0.86978349 0.99516071 TRUE
#> 20 -0.48897060 -1.47997385 TRUE
#> 21 -0.85159886 -1.23782063 TRUE
#> 22 -0.89076746 -1.11544018 TRUE
#> 23 0.65786338 -0.11367308 TRUE
#> 24 0.61055610 -1.58034754 TRUE
#> 25 0.79677764 1.07990483 TRUE
#> 26 -0.73601644 -0.25823077 TRUE
#> 27 0.26867058 1.73913899 TRUE
#> 28 0.82421091 -1.02428975 TRUE
#> 29 -0.83027258 -1.33640960 TRUE
#> 30 -0.81019972 1.35204737 TRUE
#> 31 -0.31837422 0.33334172 TRUE
#> 32 -0.40058301 1.65434459 TRUE
#> 33 -0.83950009 0.76709817 TRUE
#> 34 0.72101288 -1.50964322 TRUE
#> 35 0.42584747 1.43409713 TRUE
#> 36 0.29449219 -1.56134865 TRUE
#> 37 -0.91441562 1.07700505 TRUE
#> 38 1.11245647 0.07342734 TRUE
#> 39 -0.74769934 1.47468732 TRUE
#> 40 0.57472245 0.68069784 TRUE
#> 41 0.91745274 -1.03603578 TRUE
#> 42 0.94768581 0.92851039 TRUE
#> 43 0.32478465 -1.51816788 TRUE
#> 44 -0.32999589 0.26292621 TRUE
#> 45 0.80028009 -1.16438697 TRUE
#> 46 -0.15336178 1.63007800 TRUE
#> 47 -0.18479622 -1.65015411 TRUE
#> 48 -0.46759337 -1.64443114 TRUE
#> 49 0.46830224 -1.39242561 TRUE
#> 50 0.44997179 1.49783747 TRUE
#> 51 -0.75544834 -1.42654480 TRUE
#> 52 0.79279092 0.82730624 TRUE
#> 53 -0.13999460 1.59854231 TRUE
#> 54 -0.01267753 -1.73700166 TRUE
#> 55 0.81563094 -0.61489459 TRUE
#> 56 0.45818081 1.61078734 TRUE
#> 57 -0.83562332 1.31244172 TRUE
#> 58 -0.62928632 -1.46415314 TRUE
#> 59 -0.12493228 -0.10479610 TRUE
#> 60 -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:
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 (variances known: z values):
#> Estimate Std. Error z value Pr(>|z|)
#> 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
#> Convergence tolerance (ftol): 1e-10
check_o(mod_dome, plot = FALSE) # all orthogonal to the dome surface
#> x1 x2 x0_x1 x0_x2 z y0 alpha
#> 1 -0.83711561 3.4219365 -0.8382999 3.42677761 4.835096 4.850514 90
#> 2 0.23188986 -1.2045283 0.2337108 -1.21398705 5.766069 5.868957 90
#> 3 -2.36755994 1.4573945 -2.3738515 1.46126741 5.278940 5.310609 90
#> 4 1.70344010 -1.6496203 1.7205199 -1.66616040 5.375975 5.498796 90
#> 5 3.92288043 -1.4911203 3.9194902 -1.48983159 4.296823 4.288477 90
#> 6 2.79359165 0.8216875 2.7915709 0.82109313 5.253156 5.244614 90
#> 7 -3.91908376 -0.5925722 -3.9094501 -0.59111562 4.534974 4.509969 90
#> 8 3.13655325 -2.5007474 3.1338681 -2.49860656 4.470538 4.461936 90
#> 9 -3.85673158 -1.1840677 -3.8003838 -1.16676817 4.640786 4.490966 90
#> 10 -2.84920792 0.7445927 -2.7856436 0.72798125 5.531623 5.261489 90
#> 11 3.86821263 -1.0140447 3.8553437 -1.01067113 4.515615 4.481953 90
#> 12 -3.33387788 0.9877011 -3.4313109 1.01656683 4.505575 4.813080 90
#> 13 -1.64410177 -3.3587705 -1.6266304 -3.32307784 4.834704 4.720621 90
#> 14 -0.14208978 -0.1385513 -0.1431707 -0.13960532 5.892589 5.994418 90
#> 15 2.24300737 1.7174160 2.2835033 1.74842277 5.053225 5.263237 90
#> 16 1.80074137 -1.0292353 1.7806943 -1.01777713 5.778926 5.636160 90
#> 17 -0.07866027 2.8015712 -0.0787115 2.80339609 5.293915 5.301680 90
#> 18 3.43629743 0.7847006 3.4751285 0.79356797 4.702292 4.823564 90
#> 19 -1.05023115 2.2680821 -1.0531043 2.27428689 5.415553 5.449003 90
#> 20 1.27881463 -0.5368445 1.2828458 -0.53853679 5.792902 5.834151 90
#> 21 1.75838535 -1.1849759 1.7617965 -1.18727463 5.584442 5.608876 90
#> 22 -1.30530662 -3.2733787 -1.3075791 -3.27907754 4.829846 4.848807 90
#> 23 -1.98316244 -2.2480725 -1.9808665 -2.24546988 5.210386 5.196833 90
#> 24 3.79645116 -0.0878831 3.8048323 -0.08807711 4.612597 4.635572 90
#> 25 -0.95708123 -1.2253370 -0.9569381 -1.22515381 5.794729 5.792780 90
#> 26 -0.80906547 -2.8830044 -0.8038899 -2.86456211 5.283229 5.207791 90
#> 27 -4.39410831 -1.2860121 -4.2873795 -1.25477608 4.226281 4.002118 90
#> 28 -3.42181133 -2.4042651 -3.4104205 -2.39626154 4.345182 4.312772 90
#> 29 -0.81044244 4.2284345 -0.8063347 4.20700256 4.246214 4.198094 90
#> 30 0.99272872 1.4806770 0.9818129 1.46439587 5.876176 5.732768 90
#> 31 1.64421123 -0.3584545 1.6233866 -0.35391453 5.929356 5.763020 90
#> 32 3.01638046 -2.0470794 2.9892979 -2.02869978 4.885236 4.787642 90
#> 33 -0.91609104 -2.3058334 -0.9175364 -2.30947138 5.439362 5.458709 90
#> 34 0.91539224 -3.5115324 0.9204483 -3.53092795 4.701248 4.760080 90
#> 35 2.47660429 0.2887072 2.4411483 0.28457397 5.649886 5.471092 90
#> 36 -1.93707592 0.4356447 -1.9468995 0.43785404 5.591837 5.656050 90
#> 37 0.26539203 -3.9742190 0.2583425 -3.86865240 4.856960 4.576005 90
#> 38 0.24795389 1.2587993 0.2440664 1.23906360 6.073420 5.863292 90
#> 39 -1.16519292 2.6417539 -1.1756826 2.66553627 5.137437 5.242683 90
#> 40 0.32913242 3.0327654 0.3349175 3.08607204 4.932508 5.131961 90
#> 41 2.25522411 2.6365744 2.2700172 2.65386895 4.804673 4.876171 90
#> 42 -2.31077860 1.3748423 -2.2936202 1.36463354 5.461640 5.371231 90
#> 43 -2.79048397 1.7157110 -2.7761076 1.70687175 5.093831 5.035162 90
#> 44 -2.62872648 3.2612720 -2.6694452 3.31178868 4.083255 4.228375 90
#> 45 1.73720194 2.4928685 1.7597852 2.52527534 4.999099 5.147735 90
#> 46 2.73035728 3.2209474 2.7264267 3.21631059 4.279394 4.265558 90
#> 47 0.51639161 4.1604868 0.5231234 4.21472386 4.112409 4.235030 90
#> 48 -3.44102149 0.6651114 -3.4992473 0.67636581 4.643381 4.823985 90
#> 49 -0.23401814 2.9257226 -0.2341668 2.92758071 5.222014 5.229481 90
#> 50 1.10667419 -1.4312399 1.0927226 -1.41319657 5.890025 5.725545 90
#> 51 4.10462829 1.4176944 4.1171958 1.42203515 4.094734 4.123051 90
#> 52 -2.27637961 0.8585273 -2.2612950 0.85283823 5.571675 5.489285 90
#> 53 0.32837178 -1.0236252 0.3344388 -1.04253774 5.656278 5.896974 90
#> 54 3.33502290 2.9289792 3.3091133 2.90622409 4.143172 4.071446 90
#> 55 -3.57781045 -1.4688406 -3.5373409 -1.45222617 4.739691 4.620746 90
#> 56 0.34316146 1.4846011 0.3430432 1.48408966 5.805605 5.801107 90
#> 57 2.29832612 2.3504390 2.2963338 2.34840150 5.028252 5.018456 90
#> 58 0.12137637 -4.3035416 0.1225793 -4.34619273 4.040196 4.131419 90
#> 59 2.73455602 -2.0536203 2.7580005 -2.07122688 4.812932 4.906780 90
#> 60 -2.04023287 1.9297524 -2.0653027 1.95346466 5.137004 5.281244 90
#> df/dx_x1 df/dx_x2 Ortho
#> 1 0.17282702 -0.70647720 TRUE
#> 2 -0.03982153 0.20684887 TRUE
#> 3 0.44700176 -0.27516006 TRUE
#> 4 -0.31289028 0.30300458 TRUE
#> 5 -0.91395850 0.34740341 TRUE
#> 6 -0.53227385 -0.15655930 TRUE
#> 7 0.86684633 0.13106867 TRUE
#> 8 -0.70235609 0.55998259 TRUE
#> 9 0.84622864 0.25980340 TRUE
#> 10 0.52944012 -0.13836030 TRUE
#> 11 -0.86019272 0.22549791 TRUE
#> 12 0.71291380 -0.21120923 TRUE
#> 13 0.34457972 0.70394925 TRUE
#> 14 0.02388400 0.02328922 TRUE
#> 15 -0.43385907 -0.33219530 TRUE
#> 16 -0.31594107 0.18057990 TRUE
#> 17 0.01484652 -0.52877503 TRUE
#> 18 -0.72044830 -0.16451901 TRUE
#> 19 0.19326552 -0.41737673 TRUE
#> 20 -0.21988559 0.09230765 TRUE
#> 21 -0.31410865 0.21167782 TRUE
#> 22 0.26967029 0.67626486 TRUE
#> 23 0.38116802 0.43208430 TRUE
#> 24 -0.82079019 0.01900027 TRUE
#> 25 0.16519496 0.21149668 TRUE
#> 26 0.15436295 0.55005324 TRUE
#> 27 1.07127750 0.31352797 TRUE
#> 28 0.79077230 0.55561983 TRUE
#> 29 0.19207162 -1.00212208 TRUE
#> 30 -0.17126334 -0.25544308 TRUE
#> 31 -0.28169026 0.06141130 TRUE
#> 32 -0.62437795 0.42373676 TRUE
#> 33 0.16808670 0.42308012 TRUE
#> 34 -0.19336825 0.74177913 TRUE
#> 35 -0.44619029 -0.05201411 TRUE
#> 36 0.34421539 -0.07741339 TRUE
#> 37 -0.05645589 0.84542128 TRUE
#> 38 -0.04162618 -0.21132559 TRUE
#> 39 0.22425208 -0.50842980 TRUE
#> 40 -0.06526112 -0.60134358 TRUE
#> 41 -0.46553277 -0.54425269 TRUE
#> 42 0.42701948 -0.25406347 TRUE
#> 43 0.55134428 -0.33899045 TRUE
#> 44 0.63131701 -0.78322962 TRUE
#> 45 -0.34185620 -0.49056045 TRUE
#> 46 -0.63917240 -0.75401879 TRUE
#> 47 -0.12352294 -0.99520509 TRUE
#> 48 0.72538511 -0.14020892 TRUE
#> 49 0.04477820 -0.55982238 TRUE
#> 50 -0.19085040 0.24682305 TRUE
#> 51 -0.99857980 -0.34489872 TRUE
#> 52 0.41194709 -0.15536417 TRUE
#> 53 -0.05671363 0.17679198 TRUE
#> 54 -0.81276127 -0.71380645 TRUE
#> 55 0.76553455 0.31428390 TRUE
#> 56 -0.05913410 -0.25582872 TRUE
#> 57 -0.45757776 -0.46795300 TRUE
#> 58 -0.02967003 1.05198551 TRUE
#> 59 -0.56207950 0.42211529 TRUE
#> 60 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 (variances known: z values):
#> Estimate Std. Error z value Pr(>|z|)
#> b0 1.07302 0.78335 1.37 0.171
#> 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 <2e-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
#> Convergence tolerance (ftol): 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 alpha df/dx_x1 df/dx_x2 df/dx_x3
#> 1 5.941916 27.886043 28.240127 89.99998 0.8224703 4.5133518 0.4898238
#> 2 3.844622 15.346631 14.860667 89.99997 0.8224703 1.2501739 0.3693420
#> 3 6.104887 23.441998 23.533338 89.99998 0.8224703 3.1281373 0.4967085
#> 4 6.041930 27.181682 27.391530 89.99998 0.8224703 3.3224618 0.4330173
#> 5 10.057727 20.161528 19.823824 89.99997 0.8224703 1.8994738 0.5475341
#> 6 8.005399 30.975175 30.607978 89.99998 0.8224703 4.0688165 0.4511241
#> 7 9.304363 23.323775 23.196967 89.99997 0.8224703 2.2964937 0.4227383
#> 8 8.238866 28.595970 28.727970 89.99998 0.8224703 4.4118029 0.3723977
#> 9 3.198750 13.569212 13.115711 89.99998 0.8224703 0.4941963 0.6270592
#> 10 8.802733 22.206409 22.500286 89.99998 0.8224703 3.8860590 0.6154296
#> 11 6.143165 15.740806 15.748402 89.99998 0.8224703 0.3495554 0.7374289
#> 12 7.026692 21.587987 21.588503 89.99998 0.8224703 2.9412295 0.5674266
#> 13 4.526238 16.631585 16.140990 89.99997 0.8224703 1.9417110 0.4305719
#> 14 6.410635 22.163145 22.014204 89.99997 0.8224703 2.0332780 0.3702165
#> 15 6.397725 21.813485 22.127601 89.99998 0.8224703 2.5544583 0.4664655
#> 16 8.416438 19.500265 19.879048 89.99999 0.8224703 0.5441788 0.4122759
#> 17 8.348266 19.855123 19.652650 89.99998 0.8224703 2.4461096 0.5333731
#> 18 9.498057 15.398929 15.166459 89.99998 0.8224703 0.7889833 0.4636833
#> 19 7.718005 18.723114 19.028177 89.99998 0.8224703 2.7607878 0.4389980
#> 20 8.073035 21.909928 21.926317 89.99998 0.8224703 3.4331644 0.4417050
#> 21 2.326998 17.251493 16.548070 89.99997 0.8224703 2.4570081 0.4130413
#> 22 6.645334 20.851561 20.667282 89.99998 0.8224703 3.1019739 0.3717382
#> 23 4.049517 19.317549 19.070999 89.99998 0.8224703 0.4552448 0.4222349
#> 24 5.533346 20.305866 19.901274 89.99997 0.8224703 2.7012669 0.7035683
#> 25 3.203875 17.743470 17.704712 89.99998 0.8224703 1.6597751 0.6051674
#> 26 7.522764 16.259156 16.656058 89.99999 0.8224703 0.7050440 0.4601378
#> 27 3.978168 12.135724 12.349734 89.99998 0.8224703 0.3044329 0.4398210
#> 28 6.717142 20.663850 20.523698 89.99998 0.8224703 1.3076135 0.4837249
#> 29 3.041332 9.424084 9.625482 89.99998 0.8224703 0.3333976 0.6457774
#> 30 3.559770 23.109956 22.717998 89.99998 0.8224703 3.8590244 0.3699596
#> 31 5.764902 15.933523 16.783014 89.99999 0.8224703 0.2359137 0.4707663
#> 32 6.709537 23.186823 23.287496 89.99998 0.8224703 4.0203668 0.5898133
#> 33 9.242434 18.635737 18.398532 89.99999 0.8224703 0.2723695 0.4055661
#> 34 8.081595 22.374308 21.959863 89.99997 0.8224703 2.0608257 0.3979227
#> 35 5.209295 14.685855 14.540390 89.99998 0.8224703 0.7191490 0.3748224
#> 36 2.589617 15.723523 15.802003 89.99998 0.8224703 1.2334292 0.4153815
#> 37 4.250993 27.910620 27.705309 89.99998 0.8224703 4.5047828 0.5546150
#> 38 4.152863 21.793969 21.514469 89.99997 0.8224703 1.9926612 0.4637722
#> 39 7.487795 21.991983 21.197532 89.99996 0.8224703 2.1642730 0.4482063
#> 40 6.725630 21.221799 21.214043 89.99998 0.8224703 0.4954530 0.4863548
#> 41 5.376158 19.329673 19.668965 89.99998 0.8224703 0.2986505 0.4691104
#> 42 9.703111 26.434610 26.946783 89.99998 0.8224703 3.6833146 0.4050996
#> 43 9.820502 27.952224 27.672966 89.99998 0.8224703 3.5181389 0.4484993
#> 44 3.956667 18.900434 19.705676 89.99998 0.8224703 0.4701152 0.4828253
#> 45 5.051148 18.923265 19.383639 89.99998 0.8224703 2.8642864 0.7095374
#> 46 3.564032 25.278818 25.831901 89.99998 0.8224703 3.0840863 0.3705920
#> 47 2.194115 15.854592 16.398851 89.99998 0.8224703 1.8470086 0.3919488
#> 48 4.370773 21.042927 21.294526 89.99998 0.8224703 2.9354458 0.3941640
#> 49 2.675332 17.525823 17.456728 89.99998 0.8224703 0.6718875 0.3883172
#> 50 6.614697 11.890151 12.151032 89.99998 0.8224703 0.8949652 0.7983643
#> 51 8.914812 23.233925 23.035681 89.99997 0.8224703 2.2660062 0.4374101
#> 52 7.808095 30.100086 30.306870 89.99998 0.8224703 4.2587338 0.4597074
#> 53 2.942249 18.536580 18.512026 89.99998 0.8224703 1.0947619 0.4577393
#> 54 8.160596 32.028120 31.906469 89.99998 0.8224703 4.8495599 0.9624854
#> 55 3.590022 19.925966 19.691078 89.99998 0.8224703 2.9512205 0.3936091
#> 56 4.307153 18.490069 17.835004 89.99996 0.8224703 1.6014381 0.8202846
#> 57 4.944996 13.059103 13.800621 89.99999 0.8224703 1.2130021 0.4674748
#> 58 5.807526 17.299519 17.089450 89.99998 0.8224703 0.9021189 0.5762370
#> 59 2.216883 19.528051 19.590981 89.99997 0.8224703 1.2442570 0.4144959
#> 60 8.161779 17.151990 17.301609 89.99998 0.8224703 1.8161261 0.7097161
#> df/dx_x4 Ortho
#> 1 0.4615200 TRUE
#> 2 0.6613174 TRUE
#> 3 0.4509337 TRUE
#> 4 0.4549652 TRUE
#> 5 0.2897370 TRUE
#> 6 0.3557680 TRUE
#> 7 0.3109200 TRUE
#> 8 0.3467777 TRUE
#> 9 0.7630444 TRUE
#> 10 0.3268306 TRUE
#> 11 0.4485173 TRUE
#> 12 0.3991474 TRUE
#> 13 0.5797493 TRUE
#> 14 0.4323291 TRUE
#> 15 0.4330835 TRUE
#> 16 0.3402383 TRUE
#> 17 0.3427195 TRUE
#> 18 0.3051834 TRUE
#> 19 0.3674961 TRUE
#> 20 0.3531159 TRUE
#> 21 0.9629802 TRUE
#> 22 0.4190573 TRUE
#> 23 0.6344831 TRUE
#> 24 0.4903816 TRUE
#> 25 0.7621142 TRUE
#> 26 0.3759148 TRUE
#> 27 0.6435768 TRUE
#> 28 0.4151580 TRUE
#> 29 0.7927666 TRUE
#> 30 0.7026304 TRUE
#> 31 0.4735964 TRUE
#> 32 0.4155675 TRUE
#> 33 0.3128000 TRUE
#> 34 0.3527831 TRUE
#> 35 0.5159737 TRUE
#> 36 0.8925277 TRUE
#> 37 0.6101385 TRUE
#> 38 0.6217578 TRUE
#> 39 0.3774635 TRUE
#> 40 0.4147019 TRUE
#> 41 0.5024708 TRUE
#> 42 0.2993366 TRUE
#> 43 0.2960891 TRUE
#> 44 0.6463684 TRUE
#> 45 0.5294587 TRUE
#> 46 0.7019743 TRUE
#> 47 1.0030424 TRUE
#> 48 0.5965311 TRUE
#> 49 0.8717125 TRUE
#> 50 0.4207433 TRUE
#> 51 0.3231360 TRUE
#> 52 0.3637373 TRUE
#> 53 0.8126916 TRUE
#> 54 0.3497407 TRUE
#> 55 0.6979995 TRUE
#> 56 0.6036820 TRUE
#> 57 0.5389126 TRUE
#> 58 0.4706310 TRUE
#> 59 0.9959434 TRUE
#> 60 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 check:
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")
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://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.439982 0.900000
deviance_o(mod_bnd) # 0.1919, lower than the 0.2674 of the original ODRPACK
#> [1] 0.1918681
#> 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
#> Convergence tolerance (ftol): 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
#> Convergence tolerance (ftol): 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:
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() is the log-likelihood of the
errors-in-variables model: each observation contributes \(1 + p\) Gaussian error components (response
and predictors) with the precisions given by weights,
sigma_x and sigma_y. With known variances the
scale is 1, otherwise it is estimated, and the foot points are maximized
out without being counted as parameters.
AIC()/BIC() built from it compare
onls() fits to the same data with the
same error variances, e.g. two model functions. (The ratio of
predictor to response variance cannot be chosen this way: the likelihood
is unbounded as the ratio tends to zero.) Here we compare the
three-parameter logistic model with a four-parameter one, both with
known error variances:
mod_3 <- 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)
mod_4 <- onls(density ~ y0 + (Asym - y0) / (1 + exp((xmid - log(conc)) / scal)),
data = DNase1, start = list(y0 = 0, Asym = 3, xmid = 0, scal = 1),
sigma_x = 0.05, sigma_y = 0.1)
c(AIC(logLik_o(mod_3)), AIC(logLik_o(mod_4)))
#> [1] -86.57912 -85.70925
9.1 Confidence intervals: confint()
confint() offers three kinds of intervals through
method:
"profile" (default): each parameter is held fixed on a
grid around its estimate while all other parameters and all foot points
are re-optimized (profile()); the limits are where the
increase of the orthogonal objective reaches the \(t\) (estimated variances) or normal (known
variances) quantile. Profile intervals follow the shape of the
objective, so they can be asymmetric."wald": estimate \(\pm\) quantile \(\times\) standard error, from
vcov(). Instant, and accurate for well-determined
parameters, but symmetric by construction."boot": nonparametric case bootstrap with
k refits. The rows are resampled together with their
weights, sigma_x and sigma_y; at
most max_tries refits are attempted.
For the DNase fit of Section 3, the three methods tell different
stories about Asym and xmid. With
Asym held fixed at 10, 20, 50 or 100, the other parameters
compensate and the objective rises only to \(\tau = 1.07, 1.25, 1.36, 1.39\), below the
95 % quantile \(t_{0.975}(13) = 2.16\):
the data do not bound Asym from above. The profile reports
this as NA (with a warning), the bootstrap gives an upper
quantile far beyond the data, and the Wald interval hides it, because it
is symmetric by construction. scal is well determined, and
there the three methods are much closer:
confint(mod_uni) # profile
#> Warning: Some limits are NA:
#> Asym, upper limit: the profile ran 10 cut-offs (in standard errors) from the estimate without reaching the level, so the data hardly bound the parameter on this side
#> xmid, upper limit: the profile flattens out, so the data do not bound the parameter on this side
#> 2.5 % 97.5 %
#> Asym 1.895539 NA
#> xmid 0.970884 NA
#> scal 0.831269 2.179495
confint(mod_uni, method = "wald")
#> 2.5 % 97.5 %
#> Asym 0.1017637 6.505078
#> xmid -0.1511966 4.855241
#> scal 0.6643038 1.965708
set.seed(123)
confint(mod_uni, method = "boot", k = 200)
#> 2.5 % 97.5 %
#> Asym 1.3038939 7.112539e+06
#> xmid 0.4056245 2.890440e+01
#> scal 0.7192032 1.820698e+00
#> attr(,"n_boot")
#> [1] 200
plot(profile(mod_uni)) shows the profiles
themselves.
9.2 Predictions: predict()
predict() gives confidence intervals for the model
function and prediction intervals for a new measurement. In the
errors-in-variables model, a new observation has errors in its response
and its predictor, so the prediction variance adds the
effective variance of ODR to the variance of the fitted value: \(se^2 = se_{fit}^2 + s^2(\sigma_y^2 + g^T\Sigma_x
g)\), with \(g\) the slope of
the model at the new predictor value. With known variances (\(s = 1\)) the quantiles are normal,
otherwise \(t\) with
df_resid degrees of freedom. sigma_x = 0 gives
the interval for a new observation at an exactly known predictor:
nd <- data.frame(conc = c(0.1, 1, 10))
predict(mod_3, nd, interval = "confidence")
#> fit lwr upr
#> [1,] 0.08922738 0.03166254 0.1467922
#> [2,] 0.47019992 0.39683873 0.5435611
#> [3,] 1.62064781 1.52005517 1.7212405
predict(mod_3, nd, interval = "prediction")
#> fit lwr upr
#> [1,] 0.08922738 -0.1255113 0.303966
#> [2,] 0.47019992 0.2587028 0.681697
#> [3,] 1.62064781 1.4002620 1.841034
predict(mod_3, nd, interval = "prediction", sigma_x = 0)
#> fit lwr upr
#> [1,] 0.08922738 -0.1150477 0.2935024
#> [2,] 0.47019992 0.2609239 0.6794759
#> [3,] 1.62064781 1.4003447 1.8409509
9.3 Speed: the two algorithms
Both algorithms minimize the same objective and give the same
estimates; they differ in cost. For a sigmoid with \(n = 400\) observations:
set.seed(1)
n <- 400
xt <- seq(0, 10, length.out = n)
DAT_big <- data.frame(x = xt + rnorm(n, 0, 0.2), y = 5 / (1 + exp(-(xt - 5))) + rnorm(n, 0, 0.2))
t_s <- system.time(m_s <- onls(y ~ b1 / (1 + exp(-(x - b2) / b3)), data = DAT_big,
start = list(b1 = 4, b2 = 4, b3 = 1.5)))[["elapsed"]]
t_d <- system.time(m_d <- onls(y ~ b1 / (1 + exp(-(x - b2) / b3)), data = DAT_big,
start = list(b1 = 4, b2 = 4, b3 = 1.5), algorithm = "dense"))[["elapsed"]]
data.frame(algorithm = c("structured", "dense"), seconds = c(t_s, t_d),
rbind(coef(m_s), coef(m_d)))
#> algorithm seconds b1 b2 b3
#> 1 structured 0.012 5.002628 5.024735 1.015022
#> 2 dense 0.507 5.002628 5.024735 1.015022
With the default algorithm = "auto", the structured
algorithm is used for \(n > 100\)
and for fits with parameter bounds; otherwise, the compiled
minpack.lm code of the dense algorithm has less overhead
per iteration. The dense algorithm grows with \(n^3\) (about 15 s at \(n = 1000\) and 150 s at \(n = 2000\)), the structured one linearly
(about 2 s at \(n = 100\,000\)).
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 orthogonality
criterion 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).