lgspline: Fit Lagrangian Multiplier Smoothing Splines

In lgspline: Lagrangian Multiplier Smoothing Splines for Smooth Function Estimation

Fit Lagrangian Multiplier Smoothing Splines

Description

lgspline fits penalized piecewise-polynomial regression splines in a monomial basis, with smoothness imposed directly as linear constraints at partition boundaries rather than absorbed into a special spline basis.

The estimator is obtained with Lagrangian multipliers and enforces the usual smoothing restrictions:

• Equivalent fitted values at knots

• Equivalent first derivatives at knots, with respect to predictors

• Equivalent second derivatives at knots, with respect to predictors

The coefficients are penalized by the closed-form cubic smoothing-spline penalty, with optional predictor- and partition-specific modifications. The same framework extends to GLMs, shape constraints, and correlated responses, while keeping the fitted partition-wise polynomials explicit and interpretable.

Usage

lgspline(
  predictors = NULL,
  y = NULL,
  formula = NULL,
  response = NULL,
  standardize_response = TRUE,
  standardize_predictors_for_knots = TRUE,
  standardize_expansions_for_fitting = TRUE,
  family = gaussian(),
  glm_weight_function = default_glm_weight_function,
  schur_correction_function = function(X, y, B, dispersion, order_list, K, family,
    observation_weights, ...) {
     lapply(1:(K + 1), function(k) 0)
 },
  need_dispersion_for_estimation = FALSE,
  dispersion_function = function(mu, y, order_indices, family, observation_weights,
    VhalfInv, ...) {
     if (!is.null(VhalfInv)) {
         VhalfInv <-
    VhalfInv[order_indices, order_indices]
         c(mean((tcrossprod(VhalfInv, t(y -
    mu)))^2/family$variance(mu)))
     }
     else {
         c(mean((y -
    mu)^2/family$variance(mu)))
     }
 },
  K = NULL,
  custom_knots = NULL,
  cluster_on_indicators = FALSE,
  make_partition_list = NULL,
  previously_tuned_penalties = NULL,
  smoothing_spline_penalty = NULL,
  opt = TRUE,
  use_custom_bfgs = TRUE,
  delta = NULL,
  tol = 10 * sqrt(.Machine$double.eps),
  tuning_criterion = "loo",
  gcv_gamma = 1.4,
  initial_wiggle = c(2e-12, 2e-07, 2e-04, 0.2),
  initial_flat = c(0.5, 5),
  wiggle_penalty = 2e-07,
  flat_ridge_penalty = 0.5,
  unique_penalty_per_partition = TRUE,
  unique_penalty_per_predictor = TRUE,
  meta_penalty = 1e-08,
  predictor_penalties = NULL,
  partition_penalties = NULL,
  include_quadratic_terms = TRUE,
  include_cubic_terms = TRUE,
  include_quartic_terms = NULL,
  include_2way_interactions = TRUE,
  include_3way_interactions = TRUE,
  include_quadratic_interactions = FALSE,
  offset = c(),
  just_linear_with_interactions = NULL,
  just_linear_without_interactions = NULL,
  exclude_interactions_for = NULL,
  exclude_these_expansions = NULL,
  custom_basis_fxn = NULL,
  include_constrain_fitted = TRUE,
  include_constrain_first_deriv = TRUE,
  include_constrain_second_deriv = TRUE,
  include_constrain_interactions = TRUE,
  add_first_and_second_derivative_constraints = NULL,
  qr_pivot_smoothing_constraints = TRUE,
  cl = NULL,
  chunk_size = NULL,
  parallel_eigen = TRUE,
  parallel_trace = FALSE,
  parallel_aga = FALSE,
  parallel_matmult = FALSE,
  parallel_qr = FALSE,
  parallel_bfgs = FALSE,
  parallel_grideval = TRUE,
  parallel_qr_qp = FALSE,
  parallel_unconstrained = TRUE,
  parallel_find_neighbors = TRUE,
  parallel_penalty = FALSE,
  parallel_make_constraint = TRUE,
  unconstrained_fit_fxn = unconstrained_fit_default,
  keep_weighted_Lambda = FALSE,
  iterate_tune = TRUE,
  iterate_final_fit = TRUE,
  blockfit = TRUE,
  qp_score_function = function(X, y, mu, order_list, dispersion, VhalfInv,
    observation_weights, ...) {
     default_qp_score_function(X, y, mu, order_list,
    dispersion, VhalfInv, observation_weights, family, ...)
 },
  qp_observations = NULL,
  qp_Amat = NULL,
  qp_bvec = NULL,
  qp_meq = 0,
  qp_positive_derivative = FALSE,
  qp_negative_derivative = FALSE,
  qp_positive_2ndderivative = FALSE,
  qp_negative_2ndderivative = FALSE,
  qp_monotonic_increase = FALSE,
  qp_monotonic_decrease = FALSE,
  qp_range_upper = NULL,
  qp_range_lower = NULL,
  qr_pivot_inequality_constraints = FALSE,
  qp_Amat_fxn = NULL,
  qp_bvec_fxn = NULL,
  qp_meq_fxn = NULL,
  constraint_values = cbind(),
  constraint_vectors = cbind(),
  return_G = TRUE,
  return_Ghalf = TRUE,
  return_U = TRUE,
  estimate_dispersion = TRUE,
  unbias_dispersion = NULL,
  return_varcovmat = TRUE,
  exact_varcovmat = FALSE,
  return_lagrange_multipliers = FALSE,
  custom_penalty_mat = NULL,
  cluster_args = c(custom_centers = NA, nstart = 10),
  dummy_dividor = 1.2345672152894e-22,
  dummy_adder = 2.234567210529e-18,
  verbose = FALSE,
  verbose_tune = FALSE,
  dummy_fit = FALSE,
  auto_encode_factors = TRUE,
  observation_weights = NULL,
  do_not_cluster_on_these = c(),
  neighbor_tolerance = 1 + 1e-08,
  null_constraint = NULL,
  critical_value = qnorm(1 - 0.05/2),
  data = NULL,
  weights = NULL,
  no_intercept = FALSE,
  correlation_id = NULL,
  spacetime = NULL,
  correlation_structure = NULL,
  VhalfInv = NULL,
  Vhalf = NULL,
  VhalfInv_fxn = NULL,
  Vhalf_fxn = NULL,
  VhalfInv_par_init = c(),
  REML_grad = NULL,
  custom_VhalfInv_loss = NULL,
  VhalfInv_logdet = NULL,
  include_warnings = TRUE,
  penalty_args = NULL,
  tuning_args = NULL,
  expansion_args = NULL,
  constraint_args = NULL,
  qp_args = NULL,
  parallel_args = NULL,
  covariance_args = NULL,
  return_args = NULL,
  glm_args = NULL,
  ...
)

Arguments

predictors	Default: NULL. Numeric matrix or data frame of predictor variables, or a formula when using the formula interface.
y	Default: NULL. Numeric response variable vector.
formula	Default: NULL. Optional statistical formula for model specification, supporting spl() (and the alias s()) for spline terms.
response	Default: NULL. Alternative name for response variable.
standardize_response	Default: TRUE. Logical indicator controlling whether the response variable should be centered and scaled before model fitting. Only offered for identity link functions.
standardize_predictors_for_knots	Default: TRUE. Logical flag controlling whether predictors are internally standardized for partitioning / knot placement. The exact transformation is handled inside make_partitions and depends on the effective clustering dimension.
standardize_expansions_for_fitting	Default: TRUE. Logical switch to standardize polynomial basis expansions during model fitting. Design matrices, variance-covariance matrices, and coefficients are backtransformed after fitting. \mathbf{U} and \mathbf{G} remain on the transformed scale; B_raw corresponds to coefficients on the expansion-standardized scale.
family	Default: gaussian(). GLM family specifying the error distribution and link function. Minimally requires: family name, link name, linkfun, linkinv, variance.
glm_weight_function	Default: GLM working weight family$mu.eta(eta)^2 / family$variance(mu), optionally multiplied by observation_weights.
schur_correction_function	Default: function returning list of zeros. Computes Schur complements \mathbf{S} added to \mathbf{G}: \mathbf{G}^{*} = (\mathbf{G}^{-1} + \mathbf{S})^{-1}.
need_dispersion_for_estimation	Default: FALSE. Logical indicator specifying whether a dispersion parameter is required for coefficient estimation (e.g. Weibull AFT).
dispersion_function	Default: function returning mean squared residuals. Custom function for estimating the exponential dispersion parameter.
K	Default: NULL. Integer specifying the number of knot locations. Intuitively, total partitions minus 1.
custom_knots	Default: NULL. Optional matrix providing user-specified knot locations in 1-D.
cluster_on_indicators	Default: FALSE. Logical flag for whether indicator variables should be used for clustering knot locations.
make_partition_list	Default: NULL. Optional list allowing direct specification of custom partition assignments. The make_partition_list returned by one model can be supplied here to reuse knot locations.
previously_tuned_penalties	Default: NULL. Optional list of pre-computed penalty components from a previous model fit.
smoothing_spline_penalty	Default: NULL. Optional custom smoothing spline penalty matrix.
opt	Default: TRUE. Logical switch controlling automatic penalty optimization.
use_custom_bfgs	Default: TRUE. Selects between a native damped-BFGS implementation with closed-form gradients or base R's BFGS with finite-difference gradients. The native path is usually faster, while the finite-difference fallback can be preferable when tuning under exact LOO and the leverage derivative is numerically noisy.
delta	Default: NULL. Numeric pseudocount for stabilizing optimization in non-identity link function scenarios.
tol	Default: 10*sqrt(.Machine$double.eps). Numeric convergence tolerance.
tuning_criterion	Default: "loo". Character scalar selecting the tuning criterion. Use "loo" for exact leave-one-out on the transformed tuning problem, or "gcv" for the generalized cross-validation criterion. The LOO path computes the needed hat-matrix diagonal exactly from blockwise constrained-\mathbf{G} quantities, without explicitly forming the full projection matrix or full hat matrix. In empirical diagnostics, the observation-wise derivative of the LOO leverage term can be numerically delicate even when the overall tuning criterion and fitted penalties remain well behaved; users who prefer a more conservative optimization path can set use_custom_bfgs = FALSE. Penalty gradients reuse the partition penalty-matrix derivatives \partial\ell/\partial\boldsymbol{\Lambda}_k, so additional tuned penalty directions require only trace products. For very large samples, generalized cross-validation is often the more practical choice; as a rough guideline, "gcv" is recommended once the sample size is above about 250,000.
gcv_gamma	Default: 1.4. Numeric scalar, at least 1, used only when tuning_criterion = "gcv". Multiplies the effective degrees of freedom in the GCV denominator during automatic penalty tuning. It is accepted but ignored when tuning_criterion = "loo".
initial_wiggle	Default: c(2e-12, 2e-7, 2e-4, 0.2). Numeric vector of initial grid points for wiggle penalty optimization, on the raw (non-negative) scale.
initial_flat	Default: c(0.5, 5). Numeric vector of initial grid points for ridge penalty optimization, on the raw scale (ratio of ridge to wiggle).
wiggle_penalty	Default: 2e-7. Numeric penalty on the integrated squared second derivative, governing function smoothness.
flat_ridge_penalty	Default: 0.5. Numeric flat ridge penalty for intercepts and linear terms only. Multiplied by wiggle_penalty to obtain total ridge penalty.
unique_penalty_per_partition	Default: TRUE. Logical flag allowing penalty magnitude to differ across partitions.
unique_penalty_per_predictor	Default: TRUE. Logical flag allowing penalty magnitude to differ between predictors.
meta_penalty	Default: 1e-8. Numeric regularization coefficient for predictor- and partition-specific penalties during tuning. On the raw scale, the implemented meta-penalty shrinks these penalty multipliers toward 1; the wiggle penalty receives only a tiny stabilizing penalty by default.
predictor_penalties	Default: NULL. Optional vector of custom penalties per predictor, on the raw (positive) scale.
partition_penalties	Default: NULL. Optional vector of custom penalties per partition, on the raw (positive) scale.
include_quadratic_terms	Default: TRUE. Logical switch to include squared predictor terms.
include_cubic_terms	Default: TRUE. Logical switch to include cubic predictor terms.
include_quartic_terms	Default: NULL. Includes quartic terms; when NULL, set to FALSE for single predictor and TRUE otherwise. Highly recommended for multi-predictor models to avoid over-specified constraints.
include_2way_interactions	Default: TRUE. Logical switch for linear two-way interactions.
include_3way_interactions	Default: TRUE. Logical switch for three-way interactions.
include_quadratic_interactions	Default: FALSE. Logical switch for linear-quadratic interaction terms.
offset	Default: Empty vector. Column indices/names to include as offsets. Coefficients for offset terms are automatically constrained to 1.
just_linear_with_interactions	Default: NULL. Integer or character vector specifying predictors to retain as linear terms while still allowing interactions.
just_linear_without_interactions	Default: NULL. Integer or character vector specifying predictors to retain only as linear terms without interactions. Eligible for blockfitting.
exclude_interactions_for	Default: NULL. Integer or character vector of predictors to exclude from all interaction terms.
exclude_these_expansions	Default: NULL. Character vector of basis expansions to exclude. Named columns of data, or in the form "_1_", "_2_", "_1_x_2_", "_2_^2" etc.
custom_basis_fxn	Default: NULL. Optional user-defined function for custom basis expansions. See get_polynomial_expansions.
include_constrain_fitted	Default: TRUE. Logical switch to constrain fitted values at knot points.
include_constrain_first_deriv	Default: TRUE. Logical switch to constrain first derivatives at knot points.
include_constrain_second_deriv	Default: TRUE. Logical switch to constrain second derivatives at knot points.
include_constrain_interactions	Default: TRUE. Logical switch to constrain interaction terms at knot points.
add_first_and_second_derivative_constraints	Default: NULL. Logical switch controlling how first- and second-derivative smoothness constraints are assembled. If TRUE, the corresponding first- and second-derivative rows are added before entering the equality constraint matrix. If FALSE, they are kept as separate equality constraints. If NULL, they are combined only when more than one predictor is spline-expanded, so a model with one spline effect and any remaining non-spline effects keeps the derivative constraints separate.
qr_pivot_smoothing_constraints	Default: TRUE. Logical switch to reduce the smoothness/equality constraint matrix to a linearly independent set before fitting. Disabling this keeps the original equality columns.
cl	Default: NULL. Parallel processing cluster object (use parallel::makeCluster()).
chunk_size	Default: NULL. Integer specifying custom chunk size for parallel processing.
parallel_eigen	Default: TRUE. Logical flag for parallel eigenvalue decomposition. Ignored inside tuning fits when parallel_grideval or parallel_bfgs is using the cluster.
parallel_trace	Default: FALSE. Logical flag for parallel trace computation.
parallel_aga	Default: FALSE. Logical flag for parallel \mathbf{G} and \mathbf{A} matrix operations.
parallel_matmult	Default: FALSE. Logical flag for parallel block-diagonal matrix multiplication.
parallel_qr	Default: FALSE. Logical flag for the tall-skinny least-squares and rank-reduction steps that arise in transformed constraint solves. When active and a cluster is supplied, these steps use row-chunked cross-products with small dense fallback solves instead of relying entirely on base QR; unstable cases fall back automatically to .lm.fit() or qr().
parallel_bfgs	Default: False. Logical flag for parallel evaluation of damped BFGS step candidates during penalty tuning. When active and a cluster is supplied, multiple damping factors are evaluated across workers and inner parallel flags are ignored for those fits.
parallel_grideval	Default: TRUE. Logical flag for parallel evaluation of the initial tuning grid. When active and a cluster is supplied, grid points are distributed across workers and inner parallel flags are ignored for those fits.
parallel_qr_qp	Default: FALSE. Logical flag for parallel QR pivot reduction of partition-local inequality constraint columns. When active and a cluster is supplied, each partition's reducible QP block is handled independently across workers before the final QP matrix is assembled.
parallel_unconstrained	Default: TRUE. Logical flag for parallel unconstrained MLE for non-identity-link-Gaussian models.
parallel_find_neighbors	Default: TRUE. Logical flag for parallel neighbor identification.
parallel_penalty	Default: FALSE. Logical flag for parallel penalty matrix construction.
parallel_make_constraint	Default: TRUE. Logical flag for parallel constraint matrix generation.
unconstrained_fit_fxn	Default: unconstrained_fit_default. Custom function for fitting unconstrained models per partition.
keep_weighted_Lambda	Default: FALSE. Logical flag to retain GLM weights in penalty constraints using Tikhonov parameterization. Advised for non-canonical GLMs.
iterate_tune	Default: TRUE. Logical switch for iterative optimization during penalty tuning.
iterate_final_fit	Default: TRUE. Logical switch for iterative optimization in final model fitting.
blockfit	Default: TRUE. Logical switch for backfitting with mixed spline and non-interactive linear terms. When the blockfit conditions are met, both tuning and the final fit use blockfit_solve; otherwise the code uses get_B. Any failure falls back to get_B.
qp_score_function	Default: GLM score using family$mu.eta(eta) / family$variance(mu). Used for quadratic programming, blockfit, and GEE formulations. With dense VhalfInv, the same score weight is applied after whitening. Accepts arguments "X, y, mu, order_list, dispersion, VhalfInv, observation_weights, ...".
qp_observations	Default: NULL. Either a numeric vector of observation indices at which every active built-in QP constraint is evaluated, or a named list keyed by "var:qp_<type>" (or bare "qp_<type>") giving different built-in constraints different observation subsets. The known types are qp_range_lower, qp_range_upper, qp_positive_derivative, qp_negative_derivative, qp_positive_2ndderivative, qp_negative_2ndderivative, qp_monotonic_increase, and qp_monotonic_decrease. For range and monotonicity the canonical keys are the bare forms such as "qp_range_lower" and "qp_monotonic_increase", because those constraints are not tied to a specific variable; prefixed entries are still accepted and unioned. Derivative entries dispatch per variable, so different variables and constraint types may use different subsets. Unknown keys are ignored with a warning when include_warnings = TRUE.
qp_Amat	Default: NULL. Optional pre-built QP constraint matrix. In the current pipeline its presence marks QP handling as active, but the built-in constructor does not merge it into the assembled constraint set; use qp_Amat_fxn for custom assembled constraints.
qp_bvec	Default: NULL. Optional pre-built QP right-hand side paired with qp_Amat. Like qp_Amat, it is currently treated as an advanced placeholder rather than merged into the built-in constructor.
qp_meq	Default: 0. Optional number of equality constraints paired with qp_Amat. Like qp_Amat, it is currently treated as an advanced placeholder rather than merged into the built-in constructor.
qp_positive_derivative	Default: FALSE. Require nonnegative first derivatives. Accepts FALSE (inactive), TRUE (all predictors), or a character / integer vector selecting the predictor variables to constrain.
qp_negative_derivative	Default: FALSE. Require nonpositive first derivatives. Same input types as qp_positive_derivative; may be used simultaneously on different predictors.
qp_positive_2ndderivative	Default: FALSE. Require nonnegative second derivatives (convexity). Same input types as qp_positive_derivative.
qp_negative_2ndderivative	Default: FALSE. Require nonpositive second derivatives (concavity). Same input types as qp_positive_derivative.
qp_monotonic_increase	Default: FALSE. Logical only. Require fitted values to be nondecreasing in observation order.
qp_monotonic_decrease	Default: FALSE. Logical only. Require fitted values to be nonincreasing in observation order.
qp_range_upper	Default: NULL. Optional upper bound on constrained fitted values.
qp_range_lower	Default: NULL. Optional lower bound on constrained fitted values.
qr_pivot_inequality_constraints	Default: FALSE. Logical switch to reduce partition-local inequality constraint columns to QR pivot columns before solving. Built-in range and monotonicity constraints are left unchanged, and more generally any inequality columns spanning multiple partitions are also left unchanged.
qp_Amat_fxn	Default: NULL. Custom function generating Amat.
qp_bvec_fxn	Default: NULL. Custom function generating bvec.
qp_meq_fxn	Default: NULL. Custom function generating meq.
constraint_values	Default: cbind(). Optional matrix encoding nonzero equality targets paired with constraint_vectors. When left empty, added equality constraints are treated as homogeneous.
constraint_vectors	Default: cbind(). Optional matrix of user-supplied equality-constraint vectors, appended to the internally generated smoothness constraints.
return_G	Default: TRUE. Logical switch to return the unscaled unconstrained variance-covariance matrix \mathbf{G}.
return_Ghalf	Default: TRUE. Logical switch to return \mathbf{G}^{1/2}.
return_U	Default: TRUE. Logical switch to return the constraint projection matrix \mathbf{U}.
estimate_dispersion	Default: TRUE. Logical flag to estimate dispersion after fitting.
unbias_dispersion	Default: NULL. Logical switch to multiply dispersion by N/(N - \mathrm{trace}(\mathbf{H})). When NULL, set to TRUE for Gaussian identity link and FALSE otherwise.
return_varcovmat	Default: TRUE. Logical switch to return the variance-covariance matrix of estimated coefficients. Needed for Wald inference.
exact_varcovmat	Default: FALSE. Logical switch to replace the default asymptotic (Bayesian posterior) variance-covariance matrix with the exact frequentist variance-covariance matrix of the constrained estimator. The asymptotic version uses the Hessian of the penalized log-likelihood: \tilde{\sigma}^{2}\mathbf{U}\mathbf{G}\mathbf{U}^{\top}. The exact version additionally corrects for the penalty's contribution as a shrinkage prior, giving: \tilde{\sigma}^{2}\mathbf{U}\mathbf{G}\mathbf{U}^{\top} - \tilde{\sigma}^{2}\mathbf{U}\mathbf{G}\boldsymbol{\Lambda}\mathbf{G}\mathbf{U}^{\top} When a correlation structure is present (VhalfInv non-NULL), \mathbf{G}_{\mathrm{correct}} replaces the block-diagonal \mathbf{G}. For Gaussian identity link (with or without correlation structure), the result is the exact variance-covariance matrix of the constrained estimate. The returned object still stores the result in varcovmat. Requires return_varcovmat = TRUE.
return_lagrange_multipliers	Default: FALSE. Logical switch to return the Lagrangian multiplier vector.
custom_penalty_mat	Default: NULL. Optional p \times p custom penalty matrix for individual partitions, replacing the default ridge on linear/intercept terms. Run with dummy_fit = TRUE first to inspect expansion structure.
cluster_args	Default: c(custom_centers = NA, nstart = 10). Named vector of arguments controlling clustering. If the first argument is not NA, it is treated as custom cluster centers (typically an (K+1) \times q matrix). Otherwise, default k-means is used.
dummy_dividor	Default: 0.00000000000000000000012345672152894. Small numeric constant to prevent division by zero.
dummy_adder	Default: 0.000000000000000002234567210529. Small numeric constant to prevent division by zero.
verbose	Default: FALSE. Logical flag to print general progress messages.
verbose_tune	Default: FALSE. Logical flag to print detailed progress during penalty tuning.
dummy_fit	Default: FALSE. Runs the full pipeline but sets coefficients to zero, allowing inspection of design matrix structure, penalty matrices, and partitioning. Replaces the deprecated expansions_only argument.
auto_encode_factors	Default: TRUE. Logical switch to automatically one-hot encode factor or character variables when using the formula interface.
observation_weights	Default: NULL. Numeric vector of observation-specific weights for generalized least squares estimation.
do_not_cluster_on_these	Default: c(). Predictor columns to exclude from clustering. Accepts numeric column indices or character column names.
neighbor_tolerance	Default: 1 + 1e-8. Numeric tolerance for determining neighboring partitions using k-means clustering. Intended for internal use.
null_constraint	Default: NULL. Alternative parameterization for a nonzero equality target when constraint_vectors is supplied and constraint_values is left empty.
critical_value	Default: qnorm(1-0.05/2). Numeric critical value for Wald confidence intervals.
data	Default: NULL. Optional data frame for formula-based model specification.
weights	Default: NULL. Alias for observation_weights.
no_intercept	Default: FALSE. Logical flag to constrain intercept to 0. Formulas with "0+" set this to TRUE automatically.
correlation_id, spacetime	Default: NULL. N-length vector and N-row matrix of cluster ids and longitudinal/spatial variables, respectively.
correlation_structure	Default: NULL. Native implementations: "exchangeable", "spatial-exponential", "squared-exponential", "ar(1)", "spherical", "gaussian-cosine", "gamma-cosine", "matern", and aliases. Estimated via REML.
VhalfInv	Default: NULL. Fixed custom N \times N square-root-inverse covariance matrix \mathbf{V}^{-1/2}. Triggers GLS with known covariance. Post-fit inference recomputed from whitened Gram matrices.
Vhalf	Default: NULL. Fixed custom N \times N square-root covariance \mathbf{V}^{1/2}. Computed as inverse of VhalfInv if not supplied.
VhalfInv_fxn	Default: NULL. Parametric function for \mathbf{V}^{-1/2}; takes single numeric vector "par", returns N \times N matrix. If VhalfInv is not supplied, correlation parameters are optimized from this function; when VhalfInv_par_init is omitted, the optimizer starts at 1e-2. Helper functions that also use correlation_id, spacetime, or matching extra arguments are supported.
Vhalf_fxn	Default: NULL. Optional function for efficient computation of \mathbf{V}^{1/2} from the same parameter vector used by VhalfInv_fxn. When omitted, Vhalf is obtained by explicit matrix inversion of VhalfInv.
VhalfInv_par_init	Default: c(). Initial parameter values for VhalfInv_fxn optimization, on unbounded transformed scale. If left empty while VhalfInv_fxn is supplied and VhalfInv is NULL, it is set internally to 1e-2.
REML_grad	Default: NULL. Function for the gradient of the negative REML (or custom loss) with respect to the parameters of VhalfInv_fxn. Takes "par", "model_fit", and "...".
custom_VhalfInv_loss	Default: NULL. Alternative to negative REML for the correlation parameter objective function. Takes "par", "model_fit", and "...".
VhalfInv_logdet	Default: NULL. Function for efficient \log|\mathbf{V}^{-1/2}| computation. Takes same "par" as VhalfInv_fxn.
include_warnings	Default: TRUE. Logical switch to control display of warnings.
penalty_args	Default: NULL. Optional named list grouping penalty-related arguments. See section "Grouped Argument Lists".
tuning_args	Default: NULL. Optional named list grouping tuning-related arguments.
expansion_args	Default: NULL. Optional named list grouping basis expansion arguments.
constraint_args	Default: NULL. Optional named list grouping constraint arguments.
qp_args	Default: NULL. Optional named list grouping quadratic programming arguments.
parallel_args	Default: NULL. Optional named list grouping parallel processing arguments.
covariance_args	Default: NULL. Optional named list grouping correlation structure arguments.
return_args	Default: NULL. Optional named list grouping return-control arguments.
glm_args	Default: NULL. Optional named list grouping GLM customization arguments.
...	Additional arguments passed to the unconstrained model fitting function.

Details

Main features include:

• Multiple predictors and interaction terms

• Various GLM families and link functions

• Correlation structures for longitudinal/clustered data

• Shape constraints via quadratic programming

• Parallel computation for large datasets

• Comprehensive inference tools

Value

A list of class "lgspline" containing model components:

y

Original response vector.

ytilde

Fitted/predicted values on the scale of the response.

X

List of design matrices \mathbf{X}_{k} for each partition k, containing basis expansions including intercept, linear, quadratic, cubic, and interaction terms as specified. Returned on the unstandardized scale.

A

Constraint matrix \mathbf{A} encoding smoothness constraints at knot points and any user-specified linear constraints. When qr_pivot_smoothing_constraints = TRUE, only a linearly independent subset of columns is retained via pivoted QR decomposition; otherwise the original equality columns are kept.

B

List of fitted coefficients \boldsymbol{\beta}_{k} for each partition k on the original, unstandardized scale of the predictors and response.

B_raw

List of fitted coefficients for each partition on the predictor-and-response standardized scale.

K

Number of interior knots with one predictor (number of partitions minus 1 with > 1 predictor).

p

Number of basis expansions of predictors per partition.

q

Number of predictor variables.

P

Total number of coefficients (p \times (K+1)).

N

Number of observations.

penalties

List containing optimized penalty matrices and components:

• Lambda: Baseline per-partition penalty matrix corresponding to \lambda_w(\boldsymbol{\Lambda}_s + \lambda_r\boldsymbol{\Lambda}_r + \sum_l \lambda_{l,k}\boldsymbol{\Lambda}_{l,k}) before any partition-specific block is added.

• L1: Implementation label for the curvature penalty matrix \boldsymbol{\Lambda}_s.

• L2: Implementation label for the ridge penalty matrix \boldsymbol{\Lambda}_r.

• L_predictor_list: Implementation lists for additional penalty matrices \boldsymbol{\Lambda}_{l,k} associated with predictor-specific tuning directions.

• L_partition_list: Implementation lists for additional penalty matrices \boldsymbol{\Lambda}_{l,k} associated with partition-specific tuning directions.

knot_scale_transf

Function for transforming predictors to standardized scale used for knot placement.

knot_scale_inv_transf

Function for transforming standardized predictors back to original scale.

knots

Matrix of knot locations on original unstandardized predictor scale for one predictor.

partition_codes

Vector assigning observations to partitions.

partition_bounds

Vector or matrix specifying the boundaries between partitions.

knot_expand_function

Internal function for expanding data according to partition structure.

predict

Function for generating predictions on new data. For multi-predictor models, take_first_derivatives = TRUE, take_second_derivatives returns derivatives as a named list of components per predictor variable, rather than a concatenated vector. When new_predictors contains columns not present in the data, extraneous columns are silently dropped before prediction.

assign_partition

Function for assigning new observations to partitions.

family

GLM family object specifying the error distribution and link function.

estimate_dispersion

Logical indicating whether dispersion parameter was estimated.

unbias_dispersion

Logical indicating whether dispersion estimates should be unbiased.

backtransform_coefficients

Function for converting standardized coefficients to original scale.

forwtransform_coefficients

Function for converting coefficients to standardized scale.

mean_y, sd_y

Mean and standard deviation of response if standardized.

og_order

Original ordering of observations before partitioning.

order_list

List containing observation indices for each partition.

constraint_values, constraint_vectors

Matrices specifying linear equality constraints if provided.

make_partition_list

List containing partition information for > 1-D cases.

expansion_scales

Vector of scaling factors used for standardizing basis expansions.

take_derivative, take_interaction_2ndderivative

Functions for computing derivatives of basis expansions.

get_all_derivatives_insample

Function for computing all derivatives on training data.

numerics

Indices of numeric predictors used in basis expansions.

power1_cols, power2_cols, power3_cols, power4_cols

Column indices for linear through quartic terms.

quad_cols

Column indices for all quadratic terms (including interactions).

interaction_single_cols, interaction_quad_cols

Column indices for linear-linear and linear-quadratic interactions.

triplet_cols

Column indices for three-way interactions.

nonspline_cols

Column indices for terms excluded from spline expansion.

return_varcovmat

Logical indicating whether variance-covariance matrix was computed.

raw_expansion_names

Names of basis expansion terms.

std_X, unstd_X

Functions for standardizing/unstandardizing design matrices.

parallel_cluster_supplied

Logical indicating whether a parallel cluster was supplied.

weights

Original observation weights on the data scale. When no weights were supplied, this is a vector of ones.

G

List of unscaled partition-wise information inverses \mathbf{G}_{k} if return_G = TRUE. These are the blockwise quantities stored on the fitting scale; correlation-aware trace, posterior, and variance calculations additionally use dense GLS analogues internally when needed.

Ghalf

List of \mathbf{G}_{k}^{1/2} matrices if return_Ghalf = TRUE. As with G, dense GLS square-root factors may also be constructed internally for correlation-aware post-fit calculations.

U

Constraint projection matrix \mathbf{U} if return_U = TRUE. For K=0 and no constraints, returns identity. Otherwise, returns \mathbf{U} = \mathbf{I} - \mathbf{G}\mathbf{A}(\mathbf{A}^{\top}\mathbf{G}\mathbf{A})^{-1}\mathbf{A}^{\top}. Used for computing the variance-covariance matrix \sigma^{2}\mathbf{U}\mathbf{G}.

sigmasq_tilde

Estimated (or fixed) dispersion parameter \tilde{\sigma}^{2}. For Gaussian identity fits without correlation, this is the weighted mean squared residual with optional bias correction. When VhalfInv is non-NULL, Gaussian-identity residuals are whitened before this calculation.

trace_XUGX

Effective degrees of freedom (\mathrm{trace}(\mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^{\top})), where \mathbf{X}\mathbf{U}\mathbf{G}\mathbf{X}^{\top} serves as the "hat" matrix. When VhalfInv is non-NULL, computed as \|\mathbf{V}^{-1/2}\mathbf{X}\mathbf{U}\mathbf{G}_{\mathrm{correct}}^{1/2}\|_{F}^{2} using the full penalized GLS information.

varcovmat

Variance-covariance matrix of coefficient estimates if return_varcovmat = TRUE. Computed as \sigma^{2}(\mathbf{U}\mathbf{G}^{1/2})(\mathbf{U}\mathbf{G}^{1/2})^{\top} for numerical stability. When VhalfInv is non-NULL, uses the full \mathbf{G}_{\mathrm{correct}}^{1/2} in place of the block-diagonal \mathbf{G}^{1/2}.

lagrange_multipliers

Vector of Lagrangian multipliers if return_lagrange_multipliers = TRUE. For equality-only fits these correspond to the active columns of \mathbf{A}; when quadratic-programming constraints are active they are taken directly from solve.QP and therefore refer to the combined equality/inequality constraint system. NULL if no constraints are active (\mathbf{A} is NULL or K == 0).

VhalfInv

The \mathbf{V}^{-1/2} matrix used for implementing correlation structures, if specified.

VhalfInv_fxn, Vhalf_fxn, VhalfInv_logdet, REML_grad

Functions for generating \mathbf{V}^{-1/2}, \mathbf{V}^{1/2}, \log|\mathbf{V}^{-1/2}|, and gradient of REML if provided.

VhalfInv_params_estimates

Vector of estimated correlation parameters when using VhalfInv_fxn.

VhalfInv_params_vcov

Approximate variance-covariance matrix of estimated correlation parameters from BFGS optimization.

wald_univariate

Function for computing univariate Wald statistics and confidence intervals. Returns an S3 object of class "wald_lgspline" with dedicated print, summary, plot, coef, and confint methods. The print method uses printCoefmat() for standard R coefficient table formatting with significance stars.

critical_value

Critical value used for confidence interval construction.

generate_posterior

Function for drawing from the posterior distribution of coefficients. When VhalfInv is non-NULL, draws are from the correct joint posterior \mathbf{U}\mathbf{G}_{\mathrm{correct}}^{1/2}\mathbf{z} using the full penalized GLS information, reflecting cross-partition posterior covariance induced by off-diagonal blocks of \mathbf{V}^{-1/2}.

find_extremum

Function for optimizing the fitted function. Accepts both numeric column indices and character column names for vars. When select_vars_fl = TRUE, L-BFGS-B bounds are correctly subsetted to the optimized variables.

plot

Function for visualizing fitted curves.

qp_info

Metadata for inequality-constrained solves, or NULL. Includes active constraints, multipliers, and a method string describing the solver path: "active_set" = partition-wise active-set on the standard block-diagonal path; "active_set_full" = active-set with full-system equality re-solves; "active_set_woodbury" = active-set with Woodbury equality re-solves on the correlated low-rank path; "dense_qp_gee_gaussian" = dense QP fallback for correlated Gaussian fits; "dense_qp_gee_glm" = dense SQP / dense QP-subproblem fallback for correlated non-Gaussian fits.

quadprog_list

List containing the assembled QP objects qp_Amat, qp_bvec, qp_meq, and when applicable a copy of qp_info.

.fit_call_args

List containing the arguments passed to lgspline.

The returned object has class "lgspline" and provides comprehensive tools for model interpretation, inference, prediction, and visualization. All coefficients and predictions can be transformed between standardized and original scales using the provided transformation functions. The object includes both frequentist and Bayesian inference capabilities through Wald statistics and posterior sampling. S3 methods logLik.lgspline and confint.lgspline are available for standard log-likelihood extraction and confidence interval computation, respectively. Advanced customization options are available for analyzing arbitrarily complex study designs.

Response and Predictor Setup

These arguments control the primary data inputs and the initial standardization steps applied before knot placement and fitting.

GLM Customization

These options let you override the default GLM working-weight, dispersion, and partition-wise unconstrained fitting behavior.

Knots and Partitioning

These arguments determine how the predictor space is partitioned and how knot locations are chosen or reused.

Penalty

These arguments control the per-partition penalty \boldsymbol{\Lambda}_k = \lambda_w(\boldsymbol{\Lambda}_s + \lambda_r\boldsymbol{\Lambda}_r + \sum_{l=1}^{L}\lambda_{l,k}\boldsymbol{\Lambda}_{l,k}), and hence the full block-diagonal penalty \boldsymbol{\Lambda} = \mathrm{blockdiag}(\boldsymbol{\Lambda}_0, \ldots, \boldsymbol{\Lambda}_K). By default, tuning uses exact leave-one-out on the transformed problem; modified GCV is also available.

Basis Expansions

These arguments control which polynomial and interaction terms are included in the partition-specific design matrices.

Constraints

These arguments govern the smoothness equalities at partition boundaries and any additional user-supplied linear equality constraints. Shape restrictions such as monotonicity, convexity, or range bounds are handled separately in the quadratic-programming stage below.

Quadratic Programming

These arguments activate built-in or custom linear inequality constraints on fitted values or derivatives. Internally they are assembled in the form \mathbf{C}^{\top}\boldsymbol{\beta} \ge \mathbf{c} and passed to the constrained solver only when requested.

Parallel Processing

These arguments control which computational subroutines may run in parallel and how work is chunked across cluster workers.

Tuning Control

These options control iterative updates during penalty tuning and the final constrained fit.

Return Control

These arguments determine which intermediate matrices and inferential quantities are retained in the returned fit object.

Correlation Structures

These arguments enable built-in or custom working-correlation structures for longitudinal, clustered, or spatially indexed responses. In the notation of Details, the correlated penalized information is written as \mathbf{G}^{-1} = \mathbf{G}_{\mathrm{on}}^{-1} + \mathbf{G}_{\mathrm{off}}^{-1}. When the cross-partition part is low rank, the Woodbury path factors \mathbf{G}_{\mathrm{off}}^{-1} = \mathbf{E}\mathbf{J}\mathbf{E}^{\top}, define \mathbf{N} = \mathbf{G}_{\mathrm{on}}^{1/2}\mathbf{E}, and works through \mathbf{G}^{1/2} = \mathbf{G}_{\mathrm{on}}^{1/2}\mathbf{F}^{1/2} without changing the final estimator. If that low-rank path is unavailable or fails its numerical checks, the code falls back to the dense correlated solve.

Grouped Argument Lists

For convenience, related arguments can be bundled into named lists. When a grouped argument is non-NULL, its entries overwrite the corresponding individual arguments. Individual arguments remain available for backward compatibility.

penalty_args

Groups: wiggle_penalty, flat_ridge_penalty, unique_penalty_per_partition, unique_penalty_per_predictor, meta_penalty, predictor_penalties, partition_penalties, custom_penalty_mat, previously_tuned_penalties, smoothing_spline_penalty.

tuning_args

Groups: opt, use_custom_bfgs, delta, tol, tuning_criterion, gcv_gamma, initial_wiggle, initial_flat, iterate_tune, iterate_final_fit, blockfit.

expansion_args

Groups: include_quadratic_terms, include_cubic_terms, include_quartic_terms, include_2way_interactions, include_3way_interactions, include_quadratic_interactions, just_linear_with_interactions, just_linear_without_interactions, exclude_interactions_for, exclude_these_expansions, custom_basis_fxn, offset.

constraint_args

Groups: include_constrain_fitted, include_constrain_first_deriv, include_constrain_second_deriv, include_constrain_interactions, add_first_and_second_derivative_constraints, qr_pivot_smoothing_constraints, constraint_values, constraint_vectors, null_constraint, no_intercept.

qp_args

Groups: qp_score_function, qp_observations, qp_Amat, qp_bvec, qp_meq, qp_positive_derivative, qp_negative_derivative, qp_positive_2ndderivative, qp_negative_2ndderivative, qp_monotonic_increase, qp_monotonic_decrease, qp_range_upper, qp_range_lower, qr_pivot_inequality_constraints, qp_Amat_fxn, qp_bvec_fxn, qp_meq_fxn.

parallel_args

Groups: cl, chunk_size, parallel_eigen, parallel_trace, parallel_aga, parallel_matmult, parallel_qr, parallel_bfgs, parallel_grideval, parallel_qr_qp, parallel_unconstrained, parallel_find_neighbors, parallel_penalty, parallel_make_constraint.

covariance_args

Groups: correlation_id, spacetime, correlation_structure, VhalfInv, Vhalf, VhalfInv_fxn, Vhalf_fxn, VhalfInv_par_init, REML_grad, custom_VhalfInv_loss, VhalfInv_logdet.

return_args

Groups: return_G, return_Ghalf, return_U, estimate_dispersion, unbias_dispersion, return_varcovmat, exact_varcovmat, return_lagrange_multipliers.

glm_args

Groups: glm_weight_function, schur_correction_function, need_dispersion_for_estimation, dispersion_function, unconstrained_fit_fxn, keep_weighted_Lambda.

Miscellaneous

These remaining arguments affect inference defaults, numerical safeguards, verbosity, and developer-oriented diagnostics.

See Also

• lgspline.fit for the low-level fitting interface

• logLik.lgspline for log-likelihood extraction

• confint.lgspline for confidence interval extraction

• leave_one_out for leave-one-out cross-validated predictions

• blockfit_solve for the standalone backfitting solver

• solve.QP for quadratic programming optimization

• plot_ly for interactive plotting

• kmeans for k-means clustering

• optim for general purpose optimization routines

Examples

## ## ## ## Simple Examples ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## Simulate some data, fit using default settings without tuning, and plot
set.seed(1234)
t <- runif(2500, -10, 10)
y <- 2*sin(t) + -0.06*t^2 + rnorm(length(t))
model_fit <- lgspline(t, y, opt = FALSE)
plot(t, y, main = 'Observed Data vs. Fitted Function, Colored by Partition',
     ylim = c(-10, 10))
plot(model_fit, add = TRUE)


## Repeat using logistic regression, with univariate inference shown
# and alternative function call
y <- rbinom(length(y), 1, 1/(1+exp(-std(y))))
df <- data.frame(t = t, y = y)
model_fit <- lgspline(y ~ spl(t),
                      df,
                      family = binomial())
plot(t, y, main = 'Observed Data vs Fitted Function with Formulas and Derivatives',
  ylim = c(-0.5, 1.05), cex.main = 0.8)
plot(model_fit,
     show_formulas = TRUE,
     text_size_formula = 0.65,
     legend_pos = 'bottomleft',
     legend_args = list(y.intersp = 1.1),
     add = TRUE)
## Notice how the coefficients match the formula, and expansions are
# homogenous across partitions without reparameterization
print(summary(model_fit))

## Overlay first and second derivatives of fitted function respectively
derivs <- predict(model_fit,
                  new_predictors = sort(t),
                  take_first_derivatives = TRUE,
                  take_second_derivatives = TRUE)
points(sort(t), derivs$first_deriv, col = 'gold', type = 'l')
points(sort(t), derivs$second_deriv, col = 'goldenrod', type = 'l')
legend('bottomright',
       col = c('gold','goldenrod'),
       lty = 1,
       legend = c('First Derivative', 'Second Derivative'))

## Simple 2D example - including a non-spline effect
z <- seq(-2, 2, length.out = length(y))
df <- data.frame(Predictor1 = t,
                 Predictor2 = z,
                 Response = sin(y)+0.1*z)
model_fit <- lgspline(Response ~ spl(Predictor1) + Predictor1*Predictor2,
                      df)

## Notice, while spline effects change over partitions,
# interactions and non-spline effects are constrained to remain the same
coefficients <- Reduce('cbind', coef(model_fit))
colnames(coefficients) <- paste0('Partition ', 1:(model_fit$K+1))
print(coefficients)

## One or two variables can be selected for plotting at a time
# even when >= 3 predictors are present
plot(model_fit,
      custom_title = 'Marginal Relationship of Predictor 1 and Response',
      vars = 'Predictor1',
      custom_response_lab = 'Response',
      show_formulas = TRUE,
      legend_pos = 'bottomright',
      digits = 4,
      text_size_formula = 0.5)

## 3D plots are implemented as well, retaining closed-formulas
my_plot <- plot(model_fit,
                show_formulas = TRUE,
                custom_response_lab = 'Response')
my_plot


## ## ## ## More Detailed 1D Example ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## 1D data generating functions
t <- seq(-9, 9, length.out = 1000)
slinky <- function(x) {
  (50 * cos(x * 2) -2 * x^2 + (0.25 * x)^4 + 80)
}
coil <- function(x) {
  (100 * cos(x * 2) +-1.5 * x^2 + (0.1 * x)^4 +
  (0.05 * x^3) + (-0.01 * x^5) +
     (0.00002 * x^6) -(0.000001 * x^7) + 100)
}
exponential_log <- function(x) {
  unlist(c(sapply(x, function(xx) {
    if (xx <= 1) {
      100 * (exp(xx) - exp(1))
    } else {
      100 * (log(xx))
    }
  })))
}
scaled_abs_gamma <- function(x) {
  2*sqrt(gamma(abs(x)))
}

## Composite function
fxn <- function(x)(slinky(t) +
                   coil(t) +
                   exponential_log(t) +
                   scaled_abs_gamma(t))

## Bind together with random noise
dat <- cbind(t, fxn(t) + rnorm(length(t), 0, 50))
colnames(dat) <- c('t', 'y')
t <- dat[,'t']
y <- dat[,'y']

## Fit Model, 4 equivalent ways are shown below
model_fit <- lgspline(t, y, opt = FALSE)
model_fit <- lgspline(y ~ spl(t), as.data.frame(dat), opt = FALSE)
model_fit <- lgspline(response = y, predictors = t, opt = FALSE)
model_fit <- lgspline(data = as.data.frame(dat), formula = y ~ ., opt = FALSE)

# This is not valid: lgspline(y ~ ., t)
# This is not valid: lgspline(y, data = as.data.frame(dat))
# Do not put operations in formulas, not valid: lgspline(y ~ log(t) + spl(t))

## Basic Functionality
predict(model_fit, new_predictors = rnorm(1)) # make prediction on new data
loo_vals <- suppressWarnings(head(leave_one_out(model_fit)))
loo_vals # may contain NA when leverage is too high
coef(model_fit) # extract coefficients
summary(model_fit) # model information and Wald inference
generate_posterior(model_fit) # generate draws of parameters from posterior distribution
find_extremum(model_fit, minimize = TRUE) # find the minimum of the fitted function

## Incorporate range constraints, custom knots, keep penalization identical
# across partitions and predictors
model_fit <- lgspline(y ~ spl(t),
                      unique_penalty_per_partition = FALSE,
                      unique_penalty_per_predictor = FALSE,
                      custom_knots = cbind(c(-2, -1, 0, 1, 2)),
                      data = data.frame(t = t, y = y),
                      qp_range_lower = -150,
                      qp_range_upper = 150,
                      qp_observations = sample(1:length(t), 50),
                      opt = FALSE)

## Plotting the constraints and knots
plot(model_fit,
     custom_title = 'Fitted Function Constrained to Lie Between (-150, 150)',
     cex.main = 0.75)
# knot locations
abline(v = model_fit$knots)
# lower bound from quadratic program
abline(h = -150, lty = 2)
# upper bound from quadratic program
abline(h = 150, lty = 2)
# observed data
points(t, y, cex = 0.24)

## Enforce monotonic increasing constraints on fitted values
# K = 4 => 5 partitions
t <- seq(-10, 10, length.out = 100)
y <- 5*sin(t) + t + 2*rnorm(length(t))
model_fit <- lgspline(t,
                      y,
                      K = 4,
                      qp_monotonic_increase = TRUE)
plot(t, y, main = 'Monotonic Increasing Function with Respect to Fitted Values')
plot(model_fit,
     add = TRUE,
     show_formulas = TRUE,
     legend_pos = 'bottomright',
     custom_predictor_lab = 't',
     custom_response_lab = 'y')

## Posterior draws under constraint
draw <- generate_posterior(model_fit, enforce_qp_constraints = TRUE)
pr <- predict(model_fit, B_predict = draw$post_draw_coefficients)
points(t, pr, col = 'grey')

## ## ## ## 2D Example using Volcano Dataset ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## Prep
data('volcano')
volcano_long <-
  Reduce('rbind', lapply(1:nrow(volcano), function(i){
    t(sapply(1:ncol(volcano), function(j){
      c(i, j, volcano[i,j])
    }))
  }))
colnames(volcano_long) <- c('Length', 'Width', 'Height')

## Fit, with 50 partitions
# When fitting with > 1 predictor and large K, including quartic terms
# is highly recommended, and/or dropping the second-derivative constraint.
# Otherwise, the constraints can impose all partitions to be equal, with one
# cubic function fit for all (there is not enough degrees of freedom to fit
# unique cubic functions due to the massive amount of constraints).
# Below, quartic terms are included and the constraint of second-derivative
# smoothness at knots is ignored.
model_fit <- lgspline(volcano_long[,c(1, 2)],
                      volcano_long[,3],
                      include_quadratic_interactions = TRUE,
                      K = 49,
                      opt = FALSE,
                      return_U = FALSE,
                      return_varcovmat = FALSE,
                      estimate_dispersion = TRUE,
                      return_Ghalf = FALSE,
                      return_G = FALSE,
                      include_constrain_second_deriv = FALSE,
                      unique_penalty_per_predictor = FALSE,
                      unique_penalty_per_partition = FALSE,
                      wiggle_penalty = 1e-10, # the fixed wiggle penalty
                      flat_ridge_penalty = 1e-2) # the ridge penalty / wiggle penalty

## Plotting on new data with interactive visual + formulas
new_input <- expand.grid(seq(min(volcano_long[,1]),
                             max(volcano_long[,1]),
                             length.out = 250),
                         seq(min(volcano_long[,2]),
                             max(volcano_long[,2]),
                             length.out = 250))
plot(model_fit,
     new_predictors = new_input,
     show_formulas = TRUE,
     custom_response_lab = "Height",
     custom_title = 'Volcano 3-D Map',
     digits = 2)

## Get AUC
area_under_volcano <- integrate(model_fit,
                               lower = apply(volcano_long, 2, min)[1:2],
                               upper = apply(volcano_long, 2, max)[1:2])

## ## ## ## Advanced Techniques using Trees Dataset ## ## ## ## ## ## ## ## ## ## ## ## ##
## Goal here is to introduce how lgspline works with non-canonical GLMs and
# demonstrate some custom features
data('trees')

## L1-regularization constraint function on standardized coefficients
# Bound all coefficients to be less than a certain value (l1_bound) in absolute
# magnitude such that | B^{(j)}_k | < lambda for all j = 1....p coefficients,
# and k = 1...K+1 partitions.
l1_constraint_matrix <- function(p, K) {
  ## Total number of coefficients
  P <- p * (K + 1)

  ## Create diagonal matrices for L1 constraint
  # First matrix: lamdba > -bound
  # Second matrix: -lambda > -bound
  first_diag <- diag(P)
  second_diag <- -diag(P)

  ## Combine matrices
  l1_Amat <- cbind(first_diag, second_diag)

  return(l1_Amat)
}

## Bounds absolute value of coefficients to be < l1_bound
l1_bound_vector <- function(qp_Amat,
                            scales,
                            l1_bound) {

  ## Combine matrices
  l1_bvec <- rep(-l1_bound, ncol(qp_Amat)) * c(1, scales)

  return(l1_bvec)
}

## Fit model, using predictor-response formulation, assuming
# Gamma-distributed response and custom quadratic-programming constraints
# as well as quartic terms, keeping the effect of height constant across
# partitions, and 3 partitions in total. The default GLM score and weights
# handle the non-canonical log link.
# You can modify this code for performing l1-regularization in general.
model_fit <- lgspline(
  Volume ~ spl(Girth) + Height*Girth,
  data = with(trees, cbind(Girth, Height, Volume)),
  family = Gamma(link = 'log'),
  keep_weighted_Lambda = TRUE,
  unbias_dispersion = TRUE, # multiply dispersion by N/(N-trace(XUGX^{T}))
  K = 2, # 3 partitions
  opt = FALSE, # keep penalties fixed
  unique_penalty_per_partition = FALSE,
  unique_penalty_per_predictor = FALSE,
  flat_ridge_penalty = 1e-64,
  wiggle_penalty = 1e-64,
  qp_Amat_fxn = function(N, p, K, X, colnm, scales, deriv_fxn, ...) {
    l1_constraint_matrix(p, K)
  },
  qp_bvec_fxn = function(qp_Amat, N, p, K, X, colnm, scales, deriv_fxn, ...) {
    l1_bound_vector(qp_Amat, scales, 25)
  },
  qp_meq_fxn = function(qp_Amat, N, p, K, X, colnm, scales, deriv_fxn, ...) 0
)

## Notice, interaction effect is constant across partitions as is the effect
# of Height alone
Reduce('cbind', coef(model_fit))

## Many constraints, many coefficients, and small sample size makes inference
#  using asymptotic variance-covariance matrix untrustworthy.
#  Confidence intervals are often too wide or narrow, even for "good" fit.
#  Consider bootstrapping or alternative.
print(summary(model_fit))

## Plot results
plot(model_fit, custom_predictor_lab1 = 'Girth',
     custom_predictor_lab2 = 'Height',
     custom_response_lab = 'Volume',
     custom_title = 'Girth and Height Predicting Volume of Trees',
     show_formulas = TRUE)

## Verify magnitude of unstandardized coefficients does not exceed bound (25)
print(max(abs(unlist(model_fit$B))))

## Find height and girth where tree volume is closest to 42
# Uses custom objective that minimizes MSE discrepancy between predicted
# value and 42.
# The vanilla find_extremum function can be thought of as
# using "function(mu)mu" aka the identity function as the
# objective, where mu = "f(t)", our estimated function. The derivative is then
# d_mu = "df/dt" with respect to predictors t.
# But with more creative objectives, and since we have machinery for
# df/dt already available, we can compute gradients for (and optimize)
# arbitrary differentiable functions of our predictors too.
# For any objective, differentiate w.r.t. to mu, then multiply by d_mu to
# satisfy chain rule.
# Here, we have objective function: 0.5*(mu-42)^2
# and gradient                    : (mu-42)*d_mu
# and L-BFGS-B will be used to find the height and girth that most closely
# yields a prediction of 42 within the bounds of the observed data.
# The d_mu also takes into account link function transforms automatically
# for most common link functions, and will return warning + instructions
# on how to program the link-function derivatives otherwise.

## Custom acquisition functions for Bayesian optimization could be coded here.
find_extremum(
  model_fit,
  minimize = TRUE,
  custom_objective_function = function(mu, sigma, ybest, ...){
    0.5*(mu - 42)^2
  },
  custom_objective_derivative = function(mu, sigma, ybest, d_mu, ...){
    (mu - 42) * d_mu
  }
)

## ## ## ## How to Use Formulas in lgspline ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## Demonstrates splines with multiple mixed predictors and interactions

## Generate data
n <- 2500
t <- rnorm(n)
u <- rnorm(n)
z <- sin(t)*mean(abs(u))/2

## Categorical predictors
cat1 <- rbinom(n, 1, 0.5)
cat2 <- rbinom(n, 1, 0.5)
cat3 <- rbinom(n, 1, 0.5)

## Response with mix of effects
response <- u + z + 0.1*(2*cat1 - 1)

## Continuous predictors re-named
continuous1 <- t
continuous2 <- z

## Combine data
dat <- data.frame(
  response = response,
  continuous1 = continuous1,
  continuous2 = continuous2,
  cat1 = cat1,
  cat2 = cat2,
  cat3 = cat3
)

## Example 1: Basic Model with Default Terms, No Intercept
# standardize_response = FALSE often needed when constraining intercepts to 0
fit1 <- lgspline(
  formula = response ~ 0 + spl(continuous1, continuous2) +
    cat1*cat2*continuous1 + cat3,
  K = 2,
  standardize_response = FALSE,
  data = dat
)
## Examine coefficients included
rownames(fit1$B$partition1)
## Verify intercept term is near 0 up to some numeric tolerance
abs(fit1$B[[1]][1]) < 1e-8

## Example 2: Similar Model with Intercept, Other Terms Excluded
fit2 <- lgspline(
  formula = response ~ spl(continuous1, continuous2) +
    cat1*cat2*continuous1 + cat3,
  K = 1,
  standardize_response = FALSE,
  include_cubic_terms = FALSE,
  exclude_these_expansions = c( # Not all need to actually be present
    '_batman_x_robin_',
    '_3_x_4_', # no cat1 x cat2 interaction, coded using column indices
    'continuous1xcontinuous2', # no continuous1 x continuous2 interaction
    'thejoker'
  ),
  data = dat
)
## Examine coefficients included
rownames(Reduce('cbind',coef(fit2)))
# Intercept will probably be present and non-0 now
abs(fit2$B[[1]][1]) < 1e-8

## ## ## ## Compare Inference to survreg for Weibull AFT Model Validation ##
# Only linear predictors, no knots, no penalties, using Weibull AFT Model
# The goal here is to ensure that for the special case of no spline effects
# and no knots, this implementation will be consistent with other model
# implementations.
# Also note, that when using models (like Weibull AFT) where dispersion is
# being estimated and is required for estimating beta coefficients,
# we use a schur complement correction function to adjust (or "correct") our
# variance-covariance matrix for both estimation and inference to account for
# uncertainty in estimating the dispersion.
# Typically the schur_correction_function would return a negative-definite
# matrix, as its output is elementwise added to the information matrix prior
# to inversion.
if (requireNamespace("survival", quietly = TRUE)) {
  data("pbc", package = "survival")
  df <- data.frame(na.omit(
    pbc[, c("time", "trt", "stage", "hepato", "bili", "age", "status")]
  ))

  ## Weibull AFT using lgspline, showing how some custom options can be used to
  # fit more complicated models
  model_fit <- lgspline(time ~ trt + stage + hepato + bili + age,
                        df,
                        family = weibull_family(),
                        need_dispersion_for_estimation = TRUE,
                        dispersion_function = weibull_dispersion_function,
                        glm_weight_function = weibull_glm_weight_function,
                        schur_correction_function = weibull_schur_correction,
                        unconstrained_fit_fxn = unconstrained_fit_weibull,
                        opt = FALSE,
                        wiggle_penalty = 0,
                        flat_ridge_penalty = 0,
                        K = 0,
                        status = df$status != 0)
  print(summary(model_fit))

  ## Survreg results match closely on estimates and inference for coefficients
  survreg_fit <- survival::survreg(
    survival::Surv(time, status != 0) ~ trt + stage + hepato + bili + age,
    df
  )
  print(summary(survreg_fit))

  ## sigmasq_tilde = scale^2 of survreg
  print(c(sqrt(model_fit$sigmasq_tilde), survreg_fit$scale))
}

## ## ## ## Modelling Correlation Structures ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
## Setup
n_blocks <- 200 # Number of correlation_ids (subjects)
block_size <- 5 # Size of each correlation_ids (number of repeated measures per subj.)
N <- n_blocks * block_size # total sample size (balanced here)
rho_true <- 0.25  # True correlation

## Generate predictors and mean structure
t <- seq(-9, 9, length.out = N)
true_mean <- sin(t)

## Create block compound symmetric errors = I(1-p) + Jp
errors <- Reduce('rbind',
                 lapply(1:n_blocks,
                        function(i){
                          sigma <- diag(block_size) + rho_true *
                            (matrix(1, block_size, block_size) -
                               diag(block_size))
                          matsqrt(sigma) %*% rnorm(block_size)
                        }))

## Generate response with correlated errors
y <- true_mean + errors * 0.5

## Fit model with correlation structure
# include_warnings = FALSE is a good idea here, since many proposed
# correlations will not work
model_fit <- lgspline(t,
                      y,
                      K = 4,
                      correlation_id = rep(1:n_blocks, each = block_size),
                      correlation_structure = 'exchangeable',
                      include_warnings = FALSE
)

## Assess overall fit
plot(t, y, main = 'Sinosudial Fit Under Correlation Structure')
plot(model_fit, add = TRUE, show_formulas = TRUE, custom_predictor_lab = 't')

## Compare estimated vs true correlation
# Built-in exchangeable uses rho = exp(-exp(par)), so par in (-Inf, Inf)
# maps to rho in (0, 1). Only positive correlation is supported.
rho_est <- exp(-exp(model_fit$VhalfInv_params_estimates))
print(c("True correlation:" = rho_true,
        "Estimated correlation:" = rho_est))

## Also check SD (should be close to 0.5)
print(sqrt(model_fit$sigmasq_tilde))

## Toeplitz Simulation Setup, with demonstration of custom functions
# and boilerplate. Toep is not implemented by default, because it makes
# strong assumptions on the study design and missingness that are rarely met,
# with non-obvious workarounds.
# If a GLM was to-be-fit, you would also submit a function "Vhalf_fxn" analogous
# to VhalfInv_fxn with same argument (par) and an output of an N x N matrix
# that yields the inverse of VhalfInv_fxn output.
n_blocks <- 250   # Number of correlation_ids
block_size <- 8   # Observations per correlation_id
N <- n_blocks * block_size # total sample size
sigma_true <- 2   # Marginal standard deviation

## True Toeplitz components
# This example uses a convex combination of two geometric lag kernels:
# corr(h) = mix * rho_fast^h + (1 - mix) * rho_slow^h
# which is Toeplitz and positive definite for mix in (0, 1) and
# rho_fast, rho_slow in (0, 1).
rho_fast_true <- 0.25
rho_slow_true <- 0.75
mix_true <- 0.40

## Create time and correlation_id variables
time_var <- rep(1:block_size, n_blocks)
correlation_id_var <- rep(1:n_blocks, each = block_size)

## Create nonlinear predictor-response relationship
# Not sinusoidal and not polynomial.
t_base <- seq(-2, 2, length.out = block_size)
t <- rep(t_base, n_blocks) + rnorm(N, sd = 0.10)
f_true <- function(t) {
  1.4 + 0.9 * atan(1.8 * t) + 0.8 * exp(-1.2 * (t - 0.4)^2)
}

## Generate mean structure
mu_true <- f_true(t)

## Toeplitz correlation helper
corr_from_components <- function(rho_fast, rho_slow, mix) {
  corr <- matrix(0, block_size, block_size)
  for(i in 1:block_size) {
    for(j in 1:block_size) {
      lag <- abs(i - j)
      if(lag == 0) {
        corr[i, j] <- 1
      } else {
        corr[i, j] <- mix * rho_fast^lag + (1 - mix) * rho_slow^lag
      }
    }
  }
  corr
}

## Toeplitz correlation function
# Custom functions can use any parameterization. Here we map:
#   par[1] -> rho_fast = exp(-exp(par[1]))
#   par[2] -> rho_slow = exp(-exp(par[2]))
#   par[3] -> mix      = plogis(par[3])
# so the parameter space is unconstrained, while the resulting Toeplitz
# correlation matrix remains valid.
corr_from_par <- function(par) {
  rho_fast <- exp(-exp(par[1]))
  rho_slow <- exp(-exp(par[2]))
  mix <- plogis(par[3])
  corr_from_components(rho_fast, rho_slow, mix)
}

## Create block Toeplitz errors from the same family we will fit
corr_true <- corr_from_components(rho_fast_true, rho_slow_true, mix_true)
errors <- Reduce('c',
                 lapply(1:n_blocks, function(i) {
                   c(matsqrt(corr_true) %*% rnorm(block_size))
                 }))

## Generate response with correlated errors and nonlinear covariate effect
y <- mu_true + sigma_true * errors

VhalfInv_fxn <- function(par) {
  corr <- corr_from_par(par)
  kronecker(diag(n_blocks), matinvsqrt(corr))
}

Vhalf_fxn <- function(par) {
  corr <- corr_from_par(par)
  kronecker(diag(n_blocks), matsqrt(corr))
}

## Determinant function (for efficiency)
# This avoids taking determinant of N by N matrix
VhalfInv_logdet <- function(par) {
  corr <- corr_from_par(par)
  log_det_invsqrt_corr <- -0.5 * determinant(corr, logarithm = TRUE)$modulus[1]
  n_blocks * log_det_invsqrt_corr
}

## GLM weights for REML gradient helper
# For Gaussian identity, these are all 1.
glm_weight_function <- function(mu, y, order_indices, family,
                                dispersion, observation_weights, ...) {
  rep(1, length(mu))
}

## REML gradient function
# The helper reml_grad_from_dV computes the three REML terms once dV / dpar
# is supplied. For this parameterization, dV / dpar has closed form.
REML_grad <- function(par, model_fit, ...) {
  rho_fast <- exp(-exp(par[1]))
  rho_slow <- exp(-exp(par[2]))
  mix <- plogis(par[3])

  dV1_block <- matrix(0, block_size, block_size)
  dV2_block <- matrix(0, block_size, block_size)
  dV3_block <- matrix(0, block_size, block_size)

  for(i in 1:block_size) {
    for(j in 1:block_size) {
      lag <- abs(i - j)
      if(lag > 0) {
        ## d/dpar[1] through rho_fast = exp(-exp(par[1]))
        dV1_block[i, j] <- -mix * lag * exp(par[1]) * rho_fast^lag
        ## d/dpar[2] through rho_slow = exp(-exp(par[2]))
        dV2_block[i, j] <- -(1 - mix) * lag * exp(par[2]) * rho_slow^lag
        ## d/dpar[3] through mix = plogis(par[3])
        dV3_block[i, j] <- mix * (1 - mix) * (rho_fast^lag - rho_slow^lag)
      }
    }
  }

  dV1 <- kronecker(diag(n_blocks), dV1_block)
  dV2 <- kronecker(diag(n_blocks), dV2_block)
  dV3 <- kronecker(diag(n_blocks), dV3_block)

  gradient <- numeric(3)
  gradient[1] <- lgspline::reml_grad_from_dV(dV1, model_fit,
                                   glm_weight_function, ...)
  gradient[2] <- reml_grad_from_dV(dV2, model_fit,
                                   glm_weight_function, ...)
  gradient[3] <- reml_grad_from_dV(dV3, model_fit,
                                   glm_weight_function, ...)
  gradient
}

## Visualization
plot(t, y, col = correlation_id_var,
     main = 'Simulated Data with Toeplitz Correlation')

## Fit model with custom Toeplitz
model_fit <- lgspline(
  response = y,
  predictors = t,
  K = 4,
  standardize_response = FALSE,
  VhalfInv_fxn = VhalfInv_fxn,
  Vhalf_fxn = Vhalf_fxn,
  VhalfInv_logdet = VhalfInv_logdet,
  REML_grad = REML_grad,
  VhalfInv_par_init = c(0, -1, 0),
  include_warnings = FALSE
)

## Print comparison of true and estimated correlations
lag_values <- 1:(block_size - 1)
corr_true_by_lag <- sapply(lag_values, function(h) {
  mix_true * rho_fast_true^h + (1 - mix_true) * rho_slow_true^h
})
rho_fast_est <- exp(-exp(model_fit$VhalfInv_params_estimates[1]))
rho_slow_est <- exp(-exp(model_fit$VhalfInv_params_estimates[2]))
mix_est <- plogis(model_fit$VhalfInv_params_estimates[3])
corr_est_by_lag <- sapply(lag_values, function(h) {
  mix_est * rho_fast_est^h + (1 - mix_est) * rho_slow_est^h
})
cat('Toeplitz Correlation Estimates by Lag:\n')
print(data.frame(
  Lag = lag_values,
  True.Correlation = round(corr_true_by_lag, 4),
  Estimated.Correlation = round(corr_est_by_lag, 4)
))

## Quantify uncertainty in Toeplitz component estimates
#  CI is constructed on the working scale and back-transformed
ci_transformed <- confint(model_fit)
rho_fast_ci <- sort(exp(-exp(ci_transformed['Correlation parameter 1', ])))
rho_slow_ci <- sort(exp(-exp(ci_transformed['Correlation parameter 2', ])))
mix_ci <- sort(plogis(ci_transformed['Correlation parameter 3', ]))
print("95% CI for Toeplitz components:")
print(rbind(rho_fast = rho_fast_ci,
            rho_slow = rho_slow_ci,
            mix = mix_ci))

## Should be ~ 2
print(sqrt(model_fit$sigmasq_tilde))

## ## ## ## Nested Fit ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
t <- sort(runif(120, -3, 3))
s <- runif(120, -2, 2)
y <- sin(t) + 0.2*s + 0.4*t*s + rnorm(120, sd = 0.25)
df <- data.frame(t = t, s = s, y = y)

## Fit full model
model_fit <- lgspline(y ~ spl(t) + t*s,
                      data = df,
                      K = 3,
                      opt = FALSE)

## Re-fit constraining the interaction term to 0
linear_constraint <- matrix(0, nrow = model_fit$P, ncol = model_fit$K + 1)
rownames(linear_constraint) <- names(unlist(model_fit$B))
for(k in 1:(model_fit$K + 1)){
  linear_constraint[paste0('partition', k, '.txs'), k] <- 1
}
null_value <- cbind(rep(0, ncol(linear_constraint)))

nested_fit <- lgspline(y ~ spl(t) + t*s,
                       data = df,
                       K = model_fit$K,
                       opt = FALSE,
                       make_partition_list = model_fit$make_partition_list,
                       previously_tuned_penalties = model_fit$penalties,
                       constraint_vectors = linear_constraint,
                       null_constraint = null_value)

## Interaction coefficient is 0 up to numerical tolerance
nested_coefficients <- Reduce('cbind', coef(nested_fit))
print(nested_coefficients['txs', ])

## ## ## ## Parallelism ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
if (requireNamespace("parallel", quietly = TRUE)) {
  ## Data generating function
  a <- runif(500000, -9, 9)
  b <- runif(500000, -9, 9)
  c <- rnorm(500000)
  d <- rpois(500000, 1)
  y <- sin(a) + cos(b) - 0.2*sqrt(a^2 + b^2) +
    abs(a) + b +
    0.5*(a^2 + b^2) +
    (1/6)*(a^3 + b^3) +
    a*b*c -
    c +
    d +
    rnorm(500000, 0, 5)

  ## Set up cores
  cl <- parallel::makeCluster(1)
  on.exit(try(parallel::stopCluster(cl), silent = TRUE), add = TRUE)

  ## This example shows some options for what operations can be parallelized
  # By default, only parallel_eigen and parallel_unconstrained are TRUE
  # parallel_unconstrained is only for GLMs, for identity link Gaussian
  # response, use parallel_matmult=TRUE to ensure parallel fitting across
  # partitions.
  # G, G^{-1/2}, and G^{1/2} are computed in parallel across each of the
  # K+1 partitions.
  # However, parallel_unconstrained only affects GLMs without corr. components
  # - it does not affect fitting here
  system.time({
    parfit <- lgspline(y ~ spl(a, b) + a*b*c + d,
                       data = data.frame(y = y,
                                         a = a,
                                         b = b,
                                         c = c,
                                         d = d),
                       cl = cl,
                       K = 1,
                       parallel_eigen = TRUE,
                       parallel_unconstrained = TRUE,
                       parallel_aga = FALSE,
                       parallel_find_neighbors = TRUE,
                       parallel_trace = FALSE,
                       parallel_matmult = TRUE,
                       parallel_qr = FALSE,
                       parallel_make_constraint = TRUE,
                       parallel_penalty = FALSE)
  })
  print(summary(parfit))
}

lgspline documentation built on May 8, 2026, 5:07 p.m.