Compositional data that are normalized so that the largest coordinate
equals one lie on a surface made of flat faces that meet at right
angles. The neighborhood of a point near an edge of this surface can
extend over two faces. Local coordinates computed by PCA directly from
the observed coordinates then describe a surface with a right-angle
fold. Rotating the neighboring face into the plane of the point’s own
face removes the fold. This vignette shows how
fit.local.likelihood() does that for a binary response.
Unfolding corrects the geometry only. It does not ensure that the response varies smoothly, that a projection to a few dimensions keeps the information about the response, or that a neighborhood contains enough cases and controls.
fit.local.likelihood() can compute the local coordinates
by weighted PCA or by weighted partial least squares (PLS1). PLS uses
the response to choose its directions, so in cross-validation the
directions must be computed from the training folds only. Both methods
return estimated probabilities of the binary outcome at the evaluation
points you supply. If the training sample was selected according to the
outcome and no correction is made, these are not estimates of risk in
the population.
Rows are max-normalized to the nonnegative surface where at least one
coordinate is one. chart.geometry="linf.simplex" uses the
shortest Euclidean-length path within that surface, computed by
dgraphs::linf.simplex.distances(). This is neither direct
Euclidean distance nor Chebyshev distance. Intermediate faces are
allowed in the path calculation. Every training point is considered
during neighbor selection; ties are ordered by training row number.
The same intrinsic neighborhoods and kernel weights are used before and after unfolding. Thus changes in the fitted chart can be attributed to its coordinates, not to a change in which observations were selected. No three-dimensional visualization coordinates enter the calculation.
X <- rbind(c(1, .8, .2), c(.8, 1, .2), c(1, .9, .2))
query <- c(1, .95, .2)
chart <- unfold.simplex.chart(X, query)
chart[c("certified", "status", "distance.error")]
#> $certified
#> [1] TRUE
#>
#> $status
#> [1] "two.faces"
#>
#> $distance.error
#> [1] 0
chart$coordinates
#> [,1] [,2] [,3]
#> 1 0.8 0.2
#> 1 1.2 0.2
#> 1 0.9 0.2For query face i and adjacent face j, unfolding sets the reference
coordinate to one and the adjacent coordinate to 2-u[i],
leaving other coordinates unchanged. Tied maxima and neighborhoods
spanning more than two faces are explicitly uncertified. A two-face
chart must also stay away from other faces, and its pairwise unfolded
distances must agree with intrinsic distances to 1e-8
absolute error. The certification is conservative: a neighborhood that
it rejects might still be unfolded in some other way, but the function
does not attempt this.
A one-dimensional seam provides a direct PCA check. Its ambient image bends through a right angle, while the certified unfolded points lie on a line. The table reports the proportion of centered variance retained by one principal component; it is a geometry check, not a response-model performance measure.
t.seam <- c(seq(-.2, -.01, length.out=20), seq(.01, .2, length.out=20))
seam <- cbind(ifelse(t.seam < 0, 1, 1-t.seam),
ifelse(t.seam < 0, 1+t.seam, 1), .2)
flat <- unfold.simplex.chart(seam, c(1, .95, .2))
variance.one <- function(x) {
s <- svd(scale(x, center=TRUE, scale=FALSE), nu=0, nv=0)$d^2
s[1]/sum(s)
}
data.frame(chart=c("Ambient PCA", "Unfolded PCA"),
first.PC.variance=c(variance.one(seam), variance.one(flat$coordinates)))
#> chart first.PC.variance
#> 1 Ambient PCA 0.8118812
#> 2 Unfolded PCA 1.0000000The following small synthetic example creates two adjacent faces and an outcome that varies smoothly across the seam. Its purpose is to illustrate the API; one realization does not compare the methods’ general performance.
set.seed(2710)
n <- 100L
t <- runif(n, -.2, .2)
U <- cbind(ifelse(t < 0, 1, 1-t), ifelse(t < 0, 1+t, 1), runif(n, .1, .3))
y <- rbinom(n, 1, plogis(-1 + 8*t))
fit <- fit.local.likelihood(U, y, X.eval=U[81:100, ],
likelihood.family="bernoulli", chart.geometry="linf.simplex",
coordinate.method="local.pca", unfold.faces=TRUE,
support.size=25, chart.dim=2, lambda.ridge=1)
fit$diagnostics$status.counts
#>
#> converged
#> 20
head(fit$diagnostics$per.eval)
#> status fallback chart.status chart.certified radius clearance
#> 1 converged FALSE one.face TRUE 0.07454397 0.1380609
#> 2 converged FALSE one.face TRUE 0.10509837 0.1337739
#> 3 converged FALSE two.faces TRUE 0.09954185 0.8908660
#> 4 converged FALSE one.face TRUE 0.12837765 0.1466193
#> 5 converged FALSE one.face TRUE 0.10142618 0.1366675
#> 6 converged FALSE one.face TRUE 0.15171540 0.1557173
#> distance.error dimension pca.variance.fraction effective.support gradient
#> 1 0.000000e+00 2 1 24.45218 3.569367e-13
#> 2 0.000000e+00 2 1 24.47495 1.110223e-15
#> 3 1.387779e-16 2 1 24.38710 2.331468e-15
#> 4 0.000000e+00 2 1 24.42415 9.628409e-12
#> 5 0.000000e+00 2 1 24.36282 5.523387e-12
#> 6 0.000000e+00 2 1 24.42656 9.458989e-12The example above is a fit evaluated on a subset of its training
points, not a prediction assessment. For a held-out assessment, exclude
those rows from X and y; the next example does
so. PLS uses the same API with
coordinate.method="local.pls".
pls <- fit.local.likelihood(U[1:80, ], y[1:80], X.eval=U[81:100, ],
likelihood.family="bernoulli", chart.geometry="linf.simplex",
coordinate.method="local.pls", unfold.faces=TRUE,
support.size=25, chart.dim=2, lambda.ridge=1)
range(pls$fitted.values)
#> [1] 0.05613232 0.85911118Coordinates are centered at the local weighted training mean. PCA
uses weighted variance; PLS1 uses weighted covariance with the training
response and deflation. Scores are scaled to unit weighted standard
deviation. Kernel times quadrature weights are normalized to the number
of neighbors. The local logistic objective penalizes slopes by
lambda.ridge * sum(slopes^2)/2, and adds
2*intercept.mass*log(1+exp(intercept))-intercept.mass*intercept
at the weighted chart center. The default half-success and half-failure
keep homogeneous neighborhoods finite. This penalty is a modeling
choice, not a confidence-interval or calibration guarantee. The query is
projected using the training center, directions and scales; its
probability is evaluated there, not at the intercept.
Convergence requires the maximum absolute gradient to be less than
tol. At machine-precision objective differences, a step
must reduce the gradient; a small step alone is not accepted as
convergence. The diagnostics retain solver failures, effective weighted
support, realized dimension, PCA variance fraction, chart certificates
and fallback counts. fallback="degree0" uses the smoothed
local proportion for an unsupported or failed chart;
fallback="na" returns NA. Compare fallback rates as well as
predictive scores.
foldid can keep all observations from a subject or
matched set together. The caller must construct these group-respecting
folds. Grid CV rebuilds neighbors, charts, PCA/PLS directions and score
scaling from each training fold. Candidate selection minimizes mean
held-out Brier loss; the table also reports log loss and fallback
counts. All rows must have finite predictions for a candidate to be
eligible. Candidate neighborhoods larger than a training fold are
excluded. A fresh outer test split is still needed to evaluate the
selected procedure. In-sample PLS scores and the best tuning score are
not such an evaluation.
fold <- rep(1:4, length.out=80)
selected <- fit.local.likelihood(U[1:80, ], y[1:80], X.eval=U[81:100, ],
likelihood.family="bernoulli", chart.geometry="linf.simplex",
coordinate.method="local.pls", unfold.faces=TRUE, support.size=20,
chart.dim.grid=c(1,2), lambda.ridge.grid=c(.1,1), foldid=fold)
selected$cv.table
#> support.size degree kernel bandwidth.multiplier lambda.ridge chart.dim
#> 1 20 1 gaussian 1 0.1 1
#> 2 20 1 gaussian 1 1.0 1
#> 3 20 1 gaussian 1 0.1 2
#> 4 20 1 gaussian 1 1.0 2
#> logloss brier fallback.count status
#> 1 0.4960247 0.1556948 0 evaluated
#> 2 0.4858100 0.1529047 0 evaluated
#> 3 0.4952880 0.1594561 0 evaluated
#> 4 0.4801041 0.1537859 0 evaluatedThe independent Python reference supplies 810 prediction values: five
synthetic scenarios (constant probability, seam signal, high-variance
signal, low-variance signal and nonlinear signal), six held-out queries
each, three methods, three chart dimensions and three ridge penalties.
Each uses 320 training observations and 60 neighbors. R predictions
agree within 2e-7, with no local fallback in these
fixtures. These deterministic comparisons test the implementation across
methods and settings; six queries per scenario are not a statistical
performance study. Geometry tests additionally cover intermediate-face
shortcuts, ties, corners, and one-dimensional supports. Fold tests
compare CV predictions to independent training-fold fits and verify that
changing held-out responses does not change their predictions when
tuning parameters are fixed.
The earlier research simulation used 150 independent datasets, with
separate training, tuning and test samples. That study tested the Python
prototype, not this R release. Package parity connects these
implementations on specified examples without claiming a repeat of the
entire simulation. Source fixtures and the generator are recorded in
dev/simplex-charts in the source repository.
The default, chart.geometry="legacy", is the smoother of
earlier versions. The simplex charts must be requested explicitly. They
support Bernoulli fits of degree zero or one, a fixed chart dimension,
and cross-validation over a grid. They do not yet support density
likelihoods, quadratic local models, or automatic choice of the chart
dimension. chart.geometry= "euclidean" provides the
weighted PCA/PLS path on ordinary coordinates without unfolding. Neither
PCA nor PLS is uniformly preferable: low-variance directions may matter
for the response, while supervised directions can overfit noisy small
supports. Local balance and held-out calibration need separate
evaluation.
Neighbor search computes one query-to-training distance vector at a time. Unfolding certification stores a support-by-support matrix. Work therefore grows with all query/training pairs and with the square of local support size; CV repeats these calculations. No approximate neighbor pruning is used. This implementation favors explicit geometry and reproducibility over speed at very large sample sizes.
See the function guide for the remaining smoothers.
The companion PCA and PLS comparison executes five synthetic examples with held-out prediction scores and figures.