Fit a Local Polynomial Smoother

Fits a local polynomial smoother (LPS) and selects its support size, polynomial degree, and kernel by cross-validation. By default, the smoother works in the observed ambient coordinates: each prediction point uses its nearest training points in Euclidean distance, centers the support at the prediction point, fits a weighted local polynomial, and uses the fitted intercept as the prediction.

Usage

fit.lps(
  X,
  y,
  foldid = NULL,
  support.grid = c(10L, 15L, 20L),
  degree.grid = 0:2,
  kernel.grid = c("gaussian", "tricube"),
  cv.folds = 5L,
  cv.seed = 1L,
  X.eval = NULL,
  coordinate.method = c("coordinates", "local.pca"),
  chart.dim = NULL,
  chart.dim.grid = NULL,
  local.chart.method = c("pca", "second.order.svd"),
  auto.chart.support.metric = c("coordinates", "operator", "both"),
  auto.chart.selection.metric = c("coordinates", "operator"),
  backend = c("auto", "R", "cpp", "cpp.local.pca"),
  design.basis = c("orthogonal.polynomial.drop", "monomial", "weighted.qr",
    "weighted.qr.drop"),
  design.drop.tol = 1e-08,
  ridge.multiplier.grid = c(0, 1e-10, 1e-08),
  ridge.condition.max = 1e+12,
  unstable.action = c("na", "mean"),
  outcome.family = c("gaussian", "bernoulli", "binomial"),
  bandwidth.multiplier.grid = 1,
  keep.cv.predictions = FALSE,
  ridge.shrinkage.target = c("zero", "local.mean"),
  selection.strategy = c("grid", "sparse_kd", "plateau_kd"),
  chart.activation = c("none", "subject.od"),
  chart.activation.response = NULL,
  chart.activation.control = list(),
  chart.dim.max = NULL,
  design.margin = 2L
)

Arguments

X

Numeric design/coordinate matrix with one observation per row.

y

Numeric response vector with length nrow(X).

foldid

Optional positive integer vector assigning rows to CV folds.

support.grid

Integer candidate neighborhood sizes.

degree.grid

Integer polynomial degrees. Currently degrees 0, 1, and 2 are supported.

kernel.grid

Candidate kernels. Supported kernels are "gaussian", "tricube", "epanechnikov", and "triangular".

cv.folds

Number of folds used when foldid is not supplied.

cv.seed

Random seed used to generate folds when foldid is not supplied.

X.eval

Optional matrix of prediction locations. Defaults to X.

coordinate.method

Local coordinate system. "coordinates" uses centered ambient coordinates. "local.pca" uses a local PCA chart.

chart.dim

Chart dimension for coordinate.method = "local.pca". If NULL, defaults to ncol(X). The special value "auto" estimates one global chart dimension from observed X only. The experimental special value "local.auto" estimates a local chart dimension separately for each prediction anchor in the ordinary local-PCA R backend.

chart.dim.grid

Optional candidate chart dimensions for experimental coupled support-size by chart-dimension selection. When supplied, coordinate.method must be "local.pca" and the evaluated candidates use scalar numeric chart dimensions after feasibility filtering. The default NULL preserves the historical single-chart.dim behavior.

local.chart.method

Local chart constructor used when coordinate.method = "local.pca". "pca" preserves the ordinary local-PCA chart path. "second.order.svd" uses an experimental curvature-corrected second-order local SVD chart and records compact chart fallback diagnostics. This option is opt-in and does not affect ambient coordinate fits.

auto.chart.support.metric

Support system used when chart.dim = "auto" or "local.auto". Included for consistency with LPL-TF and S-LPL-TF; because this smoother uses coordinate supports, "operator" is equivalent to "coordinates".

auto.chart.selection.metric

Which auto chart-dimension diagnostic to select when both diagnostics are requested.

backend

Computation backend. "auto" uses the C++ backend for coordinate.method = "coordinates" and the R reference backend for coordinate.method = "local.pca". "R" always uses the reference implementation. "cpp" requires ambient coordinates. "cpp.local.pca" is an opt-in prototype backend for coordinate.method = "local.pca" with local.chart.method = "pca".

design.basis

Local polynomial design backend. "monomial" uses the ordinary raw monomial design. "weighted.qr" solves the same design using an explicit weighted-QR numerical path. "weighted.qr.drop" first drops numerically dependent local design columns by weighted QR before solving. "orthogonal.polynomial.drop" replaces the local monomial design by a weighted-orthogonal basis spanning its estimable polynomial directions, dropping numerically rank-deficient directions first.

design.drop.tol

Relative QR tolerance used by design.basis = "weighted.qr.drop" or design.basis = "orthogonal.polynomial.drop".

ridge.multiplier.grid

Nonnegative scale-relative ridge multipliers tried by the R local-solve backend. The smallest multiplier whose penalized local normal equations pass ridge.condition.max is used.

ridge.condition.max

Maximum allowed condition number for the penalized local normal equations. Use Inf to disable this guard.

unstable.action

Action when no local solve passes the rank and condition guards. "mean" preserves the historical weighted-mean fallback. "na" returns NA, causing CV candidates with unstable predictions to be avoided.

outcome.family

Response family. "gaussian" preserves the ordinary numeric-response LPS behavior. "bernoulli" treats 0/1 responses as numeric conditional-expectation targets for \(\Pr(Y=1\mid X)\), keeps the same local least-squares fitting core, clips reported response-scale probabilities to [0,1], selects candidates by the observed CV Brier score of the clipped predictions (E2.12: the selection score is the deployed metric, which requires per-point CV predictions, so "bernoulli" – like "binomial" – always uses the R backend: backend = "auto" resolves to "R" and an explicit C++ backend is an error), and records Brier/log-loss probability diagnostics. "binomial" uses local weighted logistic polynomial fits, selects candidates by observed CV log loss with the log-loss probability truncation pinned at 1e-6 (E2.12) and with any candidate having a non-finite CV prediction scored Inf – unselectable, the same rule as the gaussian/bernoulli selection scores (E2.15) – and reports probability diagnostics on the fitted probabilities. The local logistic IRLS uses deviance step-halving (E2.14): a Newton update is accepted only if the weighted binomial deviance does not increase by more than 1e-8; otherwise the step is halved toward the current iterate (at most 30 times, after which the solve is declared non-convergent). The deviance is evaluated on the same [-35, 35] clamped linear predictor the IRLS update uses, so it is finite for every finite iterate. A solve that does not converge within the iteration cap (including under exact separation, where the unpenalized logistic MLE does not exist) falls back deterministically to the local weighted event rate under unstable.action = "mean" or to NA under unstable.action = "na"; every fallback is counted in logistic.diagnostics (fallback.path, event.rate.fallback, na.failure).

bandwidth.multiplier.grid

Nonnegative bandwidth multipliers added to the CV candidate grid. For a candidate with multiplier \(b\), the local kernel bandwidth becomes \(h = b \cdot d_{(K)}\) where \(d_{(K)}\) is the distance to the \(K\)-th nearest support neighbor, decoupling the kernel scale from the support size. The default 1 reproduces the historical behavior exactly (the bandwidth equals the support radius). Any grid other than exactly 1 requires the R reference backend: backend = "auto" then resolves to "R", and explicit backend = "cpp" or "cpp.local.pca" is an error. The selected multiplier is returned in $selected$bandwidth.multiplier and as a bandwidth.multiplier column of cv.table.

keep.cv.predictions

Logical; if TRUE, the returned object additionally carries cv.predictions, the matrix of out-of-fold CV predictions with one column per cv.table row, so selection scores can be recomputed from the actual CV predictions (E2.12). NULL on the C++ CV paths (reachable only for outcome.family = "gaussian"), which do not materialize per-point predictions. Default FALSE leaves the returned object exactly as before.

ridge.shrinkage.target

Shrinkage target of the local ridge penalty in the least-squares solve (E2.13). The default "zero" preserves the historical behavior bit-for-bit: in the orthogonal design bases the penalty acts on every transformed direction, including the constant, so a large ridge multiplier shrinks the local prediction toward 0. "local.mean" – the statistically recommended setting – leaves the constant direction unpenalized via a weighted-centering reparametrization, so a large ridge shrinks the prediction toward the local weighted mean instead, and rho = 0 remains the unpenalized weighted least-squares solve. The two settings coincide for the non-orthogonal design bases ("monomial", "weighted.qr", "weighted.qr.drop"), whose constant column is already unpenalized, and at rho = 0. The setting applies to the least-squares solve (outcome.family "gaussian" / "bernoulli"); it has no effect on the "binomial" local logistic solver, whose ridge keeps the historical structure (a warning is issued if combined).

selection.strategy

Candidate-selection strategy. "grid" preserves the full Cartesian candidate grid. The experimental "sparse_kd" strategy evaluates a sparse coupled support-size by chart-dimension skeleton when chart.dim.grid is supplied. "plateau_kd" is a geometry-only selector: from observed X only, it estimates the local PCA variance dimension over the supplied support grid, finds the initial support-size plateau where that dimension is stable from the smallest support size, aggregates plateau endpoints across representative anchors, and evaluates the resulting single (support.size, chart.dim) candidate.

chart.activation

Optional sparse-response chart activation rule. "none" preserves the ordinary LPS behavior. "subject.od" is intended for subject-occupation density workflows: a local chart whose support contains too little subject occupation mass, too few positive subject-visited points, or only fringe occupation receives fitted value zero without constructing the local polynomial fit.

chart.activation.response

Optional nonnegative response used only by chart.activation = "subject.od" to decide whether a chart is active. When omitted, y is used.

chart.activation.control

List controlling sparse chart activation. Supported fields are mass.min, n.positive.min, positive.tol, core.weight.rule, core.weight.quantile, and core.weight.min. The OD default is two positive support points and a chart-specific 0.25 weight quantile.

chart.dim.max

Optional explicit maximum chart dimension for the experimental coupled selector.

design.margin

Nonnegative prefit margin used to mark coupled (support.size, chart.dim) candidates infeasible when the full local polynomial design would be underdetermined.

Value

A list of class "lps" with response-scale fitted.values, unmodified local least-squares fitted.values.raw, selected parameters, a candidate CV table, the requested local.chart.method, and the effective chart method used for reporting. In outcome.family = "bernoulli" or "binomial" mode, fitted.values are response-scale probabilities in [0,1], fitted.values.raw are the un-clipped conditional-expectation estimates for "bernoulli" and the fitted probabilities for "binomial", cv.table$cv.brier.observed is the observed CV Brier score of the response-scale (clipped) predictions with Inf for candidates having any non-finite prediction, and probability.diagnostics records raw/clipped probability ranges, out-of-range fractions, and Brier/log-loss diagnostics. The Brier and log-loss diagnostics are defined only when the fitted predictions have the same length as the training response, which is the default X.eval = X path. In "binomial" mode, logistic.diagnostics records local logistic solve attempts, convergence statuses, fallback-path counts, event-rate fallback counts, and NA failure counts separately for CV and final fitting.

Details

The optional coordinate.method = "local.pca" mode keeps the same support and kernel weighting rule, but builds the local polynomial in a local PCA chart centered at each prediction point. With chart.dim = "auto", the chart dimension is estimated as one global scalar from observed X only, using the same shared local-PCA dimension helper used by LPL-TF and S-LPL-TF. The experimental chart.dim = "local.auto" mode estimates an input-only local chart dimension separately for each prediction anchor.

Choosing controls

First fit: supply X, y, foldid, and small support.grid, degree.grid, and kernel.grid candidates. See the installed function guide.

Selection: cv.folds, cv.seed, and selection.strategy control the search. Candidate combinations and folds multiply the fitting work; use lps.nested.cv() for an outer assessment of the tuning procedure.

Geometry: coordinate.method, chart.dim, and local chart controls describe neighborhoods and coordinate systems, not a different response model.

Numerics: backend, design.basis, rank thresholds and ridge controls govern stable local solves. Inspect lps.backend.diagnostics() before changing them. A numerical ridge is distinct from neighborhood selection.

Diagnostics: selected, cv.table, and diagnostics retain selection and local-fit information. A dense smoother.matrix needs one entry per evaluation/training pair; it is restricted to supported configurations.

References

Fan, J. and Gijbels, I. (1996). Local Polynomial Modelling and Its Applications. Chapman and Hall/CRC. ISBN 9780412983214.

Examples

X <- matrix(seq(0, 1, length.out = 20), ncol = 1)
y <- sin(2 * pi * X[, 1])
fit <- fit.lps(
  X, y, support.grid = 8L, degree.grid = 1L,
  kernel.grid = "tricube", cv.folds = 2L, backend = "R"
)
head(fit$fitted.values)
#> [1] 0.09208568 0.31586792 0.52932117 0.73582381 0.85205011 0.87594357