View source: R/fit_functions.R
| check_o | R Documentation |
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.
The check is the angle between the residual vector (from the foot point to the observation) and the model curve or surface at the foot point, which is 90^{\circ} at a stationary foot point. For weighted fits the angle is measured after each axis has been rescaled by its precision, so the same criterion covers unweighted, weighted, multivariate and correlated-error fits. See 'Details'.
check_o(object, plot = TRUE, tol_deg = 0.05, tol_kkt = NULL)
object |
an object returned from |
plot |
logical. If |
tol_deg |
tolerance in degrees: a point is orthogonal if |
tol_kkt |
no longer used (the relative KKT residual of version 0.2 has been replaced by the angle in whitened coordinates); supplying it gives a warning. |
Stationarity condition. 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
Qx_i \, (x_{0i} - x_i) = Qyy_i \, (y_i - y_{0i}) \, g_i \quad \text{with} \quad g_i = \nabla_x f(x_{0i}, \hat\theta)
Whitened coordinates. Let L_i be a factor with L_i^T L_i = Qx_i. In the coordinates (L_i x, \; Qyy_i^{1/2} y) the weighted distance d_i of onls is an ordinary Euclidean distance, and the condition above states that the residual vector
v_i = \left(L_i (x_i - x_{0i}), \; Qyy_i^{1/2} (y_i - y_{0i})\right)
is orthogonal to the tangent space of the model surface at the foot point, spanned by the columns of T_i = (L_i^T, \; Qyy_i^{1/2} g_i)^T. check_o computes the angle between v_i and this tangent space,
\cos \alpha_i = \frac{\sqrt{k_i^T (T_i^T T_i)^{-1} k_i}}{\|v_i\|}; k_i = T_i^T v_i = Qyy_i (y_i - y_{0i}) g_i - Qx_i (x_{0i} - x_i); T_i^T T_i = Qx_i + Qyy_i g_i g_i^T,
with \|v_i\|^2 = d_i^2. Hence \alpha_i = 90^{\circ} exactly when the stationarity condition holds. For a single predictor with unit precisions this is the classical angle between the residual vector and the tangent (1, m_i), m_i = df/dx, as in earlier versions. The criterion is scale-free and uses the full precision matrix, so it is valid for correlated predictor errors (a global p \times p sigma_x) and for predictors with very small sigma_x, for which the per-axis criteria of version 0.2 reported false failures.
A point is flagged orthogonal when |\alpha_i - 90^{\circ}| < tol_deg. If the residual vector has (numerically) zero length, d_i^2 \le \epsilon_{mach} (x_i^T Qx_i x_i + Qyy_i y_i^2), the observation lies on the fitted surface, the angle is undefined, and the point is reported with alpha = NA and Ortho = TRUE; this includes observations with zero weight.
Slopes. The derivatives g_i are computed 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 is sensitive to the exact position of the foot point, 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.
A data frame with one row per observation, in the original observation order (after subset and na.action):
for single-predictor models (p=1) the observed predictor x_i, the foot point x0, the observed response y_i, the foot-point response y0, the angle alpha in degrees (NA for an observation on the fitted curve), the model slope df/dx at the foot point, and a logical Ortho;
for multivariate models (p>1) the observed predictors, the foot-point predictors x0_<name>, the observed response, y0, alpha, the partial derivatives df/dx_<name> at the foot point, and Ortho.
If plot = TRUE, the angles are additionally plotted.
Andrej-Nikolai Spiess
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.