check_o: Check the fit for orthogonality of all points

In onls: Orthogonal Nonlinear Least-Squares Regression

Check the fit for orthogonality of all points

Description

Checks, for every observation, whether the foot point (x_{0i}, y_{0i}) returned by onls is a genuine stationary point of the orthogonal-distance objective that onls minimizes, i.e. whether the first-order condition \partial S/\partial \delta_i = 0 of the joint ODRPACK-type problem (see onls) is met. For unweighted (isotropic) fits, orthogonality is checked via the angle between the tangent to the model at (x_{0i}, y_{0i}) and the Euclidean vector from the foot point to (x_i, y_i), which should be 90^{\circ}. For weighted/heteroscedastic fits (see sigma_x, sigma_y, weights in onls), a plain right angle is no longer the correct criterion, so check_o instead checks the underlying first-order KKT stationarity condition directly. Multivariate models (more than one predictor) are checked one predictor axis at a time. See 'Details'.

Usage

  check_o(object, plot = TRUE, tol_deg = 0.05, tol_kkt = 0.001) 

Arguments

object	an object returned from onls.
plot	logical. If TRUE, the orthogonality diagnostic (\alpha-values for unweighted fits, or relative KKT residuals for weighted fits) is plotted for a quick overview of all points. For multivariate models, one panel per predictor is drawn.
tol_deg	tolerance in degrees for the unweighted criterion: a point is orthogonal if |\alpha_i - 90^{\circ}| < tol_deg. Default 0.05, i.e. 89.95^{\circ} < \alpha_i < 90.05^{\circ}.
tol_kkt	tolerance for the weighted criterion: a point is orthogonal if its relative KKT residual is smaller than tol_kkt. Default 0.001.

Details

Which criterion is used. check_o automatically chooses between two criteria, based on whether the fitted onls model is effectively weighted, i.e. if any of the following holds:

• object$known_sigma is TRUE,

• object$sigma_x was supplied to onls,

• object$weights has more than one distinct value.

Otherwise (default onls call, unit sigma_y, no sigma_x, constant or absent weights) the fit is treated as unweighted/isotropic. Note that explicitly passing sigma_y to onls (even sigma_y = 1) sets known_sigma = TRUE and therefore selects the weighted (KKT) criterion.

Background. onls minimizes S = \sum_i \left[ Qyy_i (y_i - f(x_i + \delta_i, \theta))^2 + \delta_i^T Qx_i \delta_i \right] jointly over the parameters and the foot-point corrections \delta_i = x_{0i} - x_i. Setting \partial S/\partial \delta_i = 0 gives the stationarity condition

Qx_i \, (x_{0i} - x_i) = Qyy_i \, (y_i - y_{0i}) \, \nabla_x f(x_{0i}, \hat\theta),

which both criteria below test. For unit precisions it is exactly the statement that the vector from the foot point to the observation is orthogonal to the model surface.

Unweighted case: tangent-angle criterion.
Let dx_i = x_i - x_{0i}, dy_i = y_i - y_{0i} be the residual vector from the foot point to the observation, and m_i = df(x, \hat\theta)/dx evaluated at x = x_{0i} the slope of the tangent, whose direction vector is (1, m_i). The function calculates the angle \alpha_i between the residual vector and the tangent,

\alpha_i[^{\circ}] = \mathrm{atan2}\left(\left|m_i \, dx_i - dy_i\right|, \; \left|dx_i + m_i \, dy_i\right|\right) \cdot \frac{180}{\pi},

which lies in [0^{\circ}, 90^{\circ}] and equals 90^{\circ} exactly when dx_i + m_i \, dy_i = 0, the stationarity condition above with unit precisions. This is algebraically identical to the classical expression \tan(\alpha_i) = |(m_i - n_i)/(1 + m_i n_i)| with the slope n_i = dy_i/dx_i of the residual vector, but avoids the division by dx_i, so that points whose foot point coincides with the observation (dx_i = 0, for example where the model slope is zero, or observations at the boundary of a flat model) are evaluated correctly instead of returning NaN. A point is flagged orthogonal when |\alpha_i - 90^{\circ}| < tol_deg. If the residual vector has (numerically) zero length, \sqrt{dx_i^2 + dy_i^2} \le \sqrt{\epsilon_{mach}}\,(|x_i| + |y_i|), the observation lies on the fitted curve, the angle is undefined, and the point is reported with alpha = NA and Ortho = TRUE. This criterion is the appropriate one only when the underlying objective is the plain (unweighted) Euclidean distance, i.e. Qyy_i = Qx_i = 1 in the notation of onls.

Weighted case: KKT-residual criterion.
Once Qyy_i and/or Qx_i rescale the response/predictor axes anisotropically (see onls), a plain right angle is no longer the correct orthogonality condition. Instead check_o checks directly the first-order stationarity condition that the foot point x_{0i} satisfies at the onls optimum,

Qyy_i\,(y_i - y_{0i})\,m_i = Qx_i\,(x_{0i} - x_i),

by computing the relative residual of this equality,

\text{rel\_resid}_i = \frac{\left|\,Qyy_i(y_i-y_{0i})m_i \,-\, Qx_i(x_{0i}-x_i)\,\right|}{\max\left(\left|Qyy_i(y_i-y_{0i})m_i\right|,\ \left|Qx_i(x_{0i}-x_i)\right|,\ \epsilon_{mach}\right)},

where \epsilon_{mach} is .Machine$double.eps (used only to avoid division by zero when both terms vanish). A point is flagged orthogonal when \text{rel\_resid}_i < tol_kkt. Qyy_i is taken from object$Q_yy; Qx_i is taken from object$Q_x (its diagonal, or observation-specific row, as applicable, see onls). For single-predictor models, Qyy_i/Qx_i (stored in the original observation order) are internally realigned to match the sorted-predictor order used elsewhere in the returned data frame (see 'Value'), so they are correctly paired with the corresponding observation.

Multivariate models (p>1).
One check (angle-based or KKT-residual-based, per the rule above) is performed per predictor axis k, using the partial derivative \partial f/\partial x_k at the foot point in place of m_i, and (for the weighted case) the k-th diagonal precision element in place of Qx_i. For a global, correlated (non-diagonal) sigma_x covariance matrix, only the diagonal entries of the precision matrix are used for this per-axis check; the true joint stationarity condition couples all predictor axes simultaneously, so this diagnostic is an approximation in that case. The overall Ortho flag for an observation is TRUE only if every predictor axis passes its individual check.

Slopes. The slopes m_i (or partial derivatives) are computed inside check_o by central differences of the model formula at the foot points, evaluated in the environment of the formula, so that user-defined functions and constants used in the model are found.

Interpretation of failures. For observations with very small residuals, the angle \alpha_i is highly sensitive to the exact position of the foot point (an error in the foot point of relative size \varepsilon produces an angle error of roughly \varepsilon divided by the length of the residual vector), so that a loose convergence tolerance in onls can lead to a few flagged points although the fit is otherwise converged. Tightening control = list(ftol = 1e-12, ptol = 1e-12) in onls sharpens the angles. Foot points held at a bound set through extend/window in onls are, by construction, not stationary.

Value

A data frame, with columns depending on whether the model has a single predictor or several, and on whether the fit is unweighted or weighted (see 'Details' for when each regime applies): For single-predictor models (p=1): the observed predictor x_i, the foot point x_{0i}, the observed response y_i, the fitted foot-point response y_{0i}, either alpha (unweighted; \alpha_i in degrees, NA for an observation that lies on the fitted curve) or rel_resid (weighted; the relative KKT residual), the model slope df/dx at the foot point, and a logical Ortho that is TRUE when |\alpha_i - 90^{\circ}| < tol_deg (unweighted; also TRUE if alpha is NA) or when the relative KKT residual is < tol_kkt (weighted). For multivariate models (p>1): the observed predictors x_{i1},\dots,x_{ip}, the foot-point predictors x0_<name> for each predictor, the observed response, the fitted foot-point response y0, either alpha_<name> or rel_resid_<name> for each predictor axis, the model partial derivative df/dx_<name> at the foot point for each predictor axis, and an overall logical Ortho that is TRUE only if the per-axis criterion passes on every predictor axis. If plot = TRUE, the diagnostic is additionally plotted, in black where Ortho = TRUE and dark red otherwise.

For single-predictor models, rows are returned in the internally-used sorted-predictor order (matching object$pred/object$resp/x0/y0), NOT the original row order of the input data; for multivariate models, sorting is a no-op and the original observation order is used.

Author(s)

Andrej-Nikolai Spiess

Examples

## Univariate fit
set.seed(123)
x <- 1:20
y <- 10 + 3*x^2
y <- sapply(y, function(a) rnorm(1, a, 0.1 * a))
DAT <- data.frame(x, y)
mod1 <- onls(y ~ a + b * x^2, data = DAT, start = list(a = 1, b = 1))
check_o(mod1)

## Multivariate fit => one check per predictor axis
set.seed(123)
n  <- 30
x1 <- runif(n, 1, 5); x2 <- runif(n, 1, 5)
z  <- 5 + 2 * x1 + 1.5 * x2^2 + rnorm(n, 0, 2)
DAT2 <- data.frame(x1 = x1, x2 = x2, z = z)
mod2 <- onls(z ~ b1 + b2 * x1 + b3 * x2^2, data = DAT2,
             start = list(b1 = 1, b2 = 1, b3 = 1))
check_o(mod2)

## Foot point equal to the observation (model slope zero at x = 0, dx = 0):
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)
DAT3 <- data.frame(x, y)
mod3 <- onls(y ~ b1 + b2 * (exp(b3 * x) - 1)^2, data = DAT3, 
             start = list(b1 = 1500, b2 = -50, b3 = -0.1))
check_o(mod3, plot = FALSE)

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