Fit Order-2 Logistic Trend Filtering

Fits a piecewise-quadratic logit curve with an L1 penalty on its third divided differences. An R interface calls native C++ dense or banded solvers. No Python or external optimizer is required.

Usage

fit.logistic.trend.filtering(
  x,
  y,
  lambda,
  weights = NULL,
  covariates = NULL,
  offset = NULL,
  nuisance.ridge = 0,
  curve.ridge = 0,
  lambda.selection = c("fixed", "cv"),
  foldid = NULL,
  control = list(),
  warn = TRUE,
  backend = c("auto", "dense", "banded"),
  cv.penalty.scale = c("training_mass", "none")
)

Arguments

x

Finite numeric predictor, with at least four distinct values.

y

Numeric response in [0, 1]. Binary observations use unit weights; binomial proportions use trial counts as weights.

lambda

Nonnegative penalties, or Inf for the quadratic limit. Required explicitly: no data-dependent default penalty is assumed.

weights

Positive observation weights, default one.

covariates

Optional finite numeric matrix of additional predictors, without an intercept. Columns are centered and scaled on training data.

offset

Optional finite vector added to the linear predictor.

nuisance.ridge, curve.ridge

Nonnegative explicit ridge penalties, default zero. See Details for their precise scaling.

lambda.selection

"fixed" returns the supplied path without selecting a penalty; "cv" selects minimum held-out weighted binomial log loss.

foldid

Integer fold labels, required for lambda.selection="cv". Supply the same label for observations that must be held out together.

control

Named list overriding tolerance=1e-9, gap.tolerance=1e-10, max.iterations=1000L, or max.workspace.bytes=1024^3. Tolerances apply to mean weighted loss.

warn

Logical; warn when a path fit fails convergence checks.

backend

Native solver: "auto" uses "banded" above 256 distinct locations or when the dense workspace estimate exceeds its limit, and "dense" otherwise. The banded backend uses compact working sets and checks every knot; it does not thin the predictor locations.

cv.penalty.scale

"training_mass" (default) multiplies each fold's lambda and ridge penalties by its training/full weight-sum ratio. This preserves their coefficients relative to mean loss. "none" retains the earlier same-numerical-penalty convention. Fixed-path fits are unchanged.

Value

An object of class logistic.trend.filtering.fit, with fits, lambda, diagnostics, training transformations, optional cv.table and cv.fold.lambda, cv.fold.penalty.factors, cv.fold.ridge, the selected cv.penalty.scale, and selected.index. Failed fits are preserved for diagnosis but standard fitted/prediction methods return NA for them unless explicitly permitted.

Details

Predictor locations are normalized to [0,1] using the training range. Lambda refers to third divided differences on that normalized scale, not on the original units of x. The minimized objective is $$\sum_i w_i[\log(1+\exp(\eta_i))-y_i\eta_i] +\lambda\|D^{(3)}f\|_1 +\frac{r_z}{2}\|\gamma\|^2 +\frac{r_f}{2\sum_i w_i}\sum_i w_i(f_i-\bar f_w)^2.$$ Here eta is the curve plus standardized covariates and offset, gamma contains standardized covariate coefficients, and the curve includes an unpenalized intercept. Third differences annihilate quadratic polynomials. Both backends preserve the falling-factorial penalty exactly. The dense solver uses a scaled global basis and a logarithmic barrier. The banded solver uses locally supported linear combinations of the same functions, a fully corrective working set, banded curve Hessians, and a small Schur complement for covariates/centering. Every stopping check scans all knots; omitted working-set coefficients remain exactly zero and violating knots are admitted until the global KKT checks pass. It does not approximate the penalty or restrict the permitted knot locations. Exact ties share a curve value but retain their own likelihood contributions, covariates, and offsets. No jittering or approximate thinning is performed. Extra covariates can be confounded with the curve; without identifying constraints, their individual coefficients are not effects identifiable separately from the curve.

Convergence requires finite arithmetic, small stationarity residual, dual feasibility and small L1 complementarity. The reported dual vector lies in [-1,1]; complementarity checks whether it is also an approximate subgradient at the returned coefficients. This is an approximate KKT check, not a separately computed primal-dual objective gap. Iteration limits and numerical failures remain explicit. Separation can prevent a finite minimizer, including in the unpenalized quadratic component. The classical Gaussian trend-filtering minimax result is not a guarantee for this model.

The dense backend is intended for moderate knot counts. The banded backend avoids observation-by-full-knot and full-knot-square arrays during fitting and prediction. Its storage is linear in observations and knots for a fixed number of covariates and penalties; computation also depends on working-set changes and Newton convergence. The conservative workspace estimate is checked before allocation. It is a solver workspace guard, not a process- wide memory limit. The quadratic initialization uses the dense solver only on the three polynomial columns and optional covariates.

CV recomputes training normalization, centering, and scaling in each fold. With the default cv.penalty.scale="training_mass", the penalty coefficients relative to weighted mean loss are held constant. The predictor range and standardized covariates still belong to each training split; scaling does not hold a common raw-coordinate operator fixed. A candidate with any failed training fit is ineligible for selection. All fits are refitted on the complete training data. Prediction outside its range uses the nearest boundary curve value; covariates and offsets continue to contribute normally. This boundary policy is not extrapolation of a quadratic polynomial.

See also

assess.logistic.trend.filtering for a global constancy test and nested-CV predictive improvement for independent binary observations.

Examples

x <- seq(0, 1, length.out = 40)
y <- rep(c(0, 1, 0, 1, 1), 8)
model <- fit.logistic.trend.filtering(x, y, lambda = c(1, 0.1))
predict(model, seq(0, 1, length.out = 10), lambda.index = 2)
#>  [1] 0.5415826 0.5544658 0.5674219 0.5804312 0.5934730 0.6065260 0.6195687
#>  [8] 0.6325792 0.6455353 0.6584149
#> attr(,"outside.range")
#>  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> attr(,"lambda")
#> [1] 0.1