integrate.lgspline: Definite Integral of a Fitted lgspline

In lgspline: Lagrangian Multiplier Smoothing Splines for Smooth Function Estimation

Definite Integral of a Fitted lgspline

Description

Given a fitted lgspline object, computes the definite integral of the fitted surface over a rectangular domain using Gauss–Legendre quadrature on predict().

Usage

## S3 method for class 'lgspline'
integrate(
  f,
  lower,
  upper,
  vars = NULL,
  initial_values = NULL,
  B_predict = NULL,
  link_scale = FALSE,
  n_quad = 50L,
  ...
)

Arguments

f	A fitted lgspline object.
lower	Numeric vector of lower bounds, one per integration variable. Scalar values are recycled.
upper	Numeric vector of upper bounds, one per integration variable. Scalar values are recycled.
vars	Default: NULL. Character or integer vector identifying which predictor(s) to integrate over. When NULL all numeric predictors are integrated simultaneously.
initial_values	Default: NULL. Numeric vector of length q supplying fixed values for predictors not among the integration variables. When NULL, non-integration predictors are held at the midpoint of their training range.
B_predict	Default: NULL. Optional list of coefficient vectors, one per partition. For additive fits, this may be the nested post_draw_coefficients list returned by generate_posterior(). When NULL the fitted coefficients are used.
link_scale	Default: FALSE. Logical; when TRUE the integral is computed on the link (linear predictor) scale \eta rather than the response scale \mu.
n_quad	Default: 50. Number of Gauss–Legendre nodes per integration dimension.
...	Additional arguments (currently unused; present for S3 method compatibility).

Value

A numeric scalar: the estimated definite integral.

Method

The integration domain is discretised into a tensor-product grid of Gauss–Legendre quadrature nodes. Predicted values at each node come from the model's predict() method, which correctly handles partition assignment and piecewise polynomial evaluation. The integral is the weighted sum

\int_{a_1}^{b_1} \cdots \int_{a_d}^{b_d} \hat{f}(\mathbf{t})\,\mathrm{d}t_1 \cdots \mathrm{d}t_d \;\approx\; \sum_{i=1}^{M} w_i\,\hat{f}(\mathbf{t}_i)

where M = n_{\mathrm{quad}}^d and each weight incorporates the Jacobian (b_j - a_j)/2 for the affine map from [-1, 1] to [a_j, b_j]. Nodes and weights on [-1, 1] are computed via the Golub–Welsch algorithm (eigenvalues of the symmetric tridiagonal Jacobi matrix).

For smooth polynomials, 30–50 nodes per dimension is typically sufficient; highly partitioned models (large K) may benefit from more. Total evaluation points scale as n_{\mathrm{quad}}^d, so problems with d \ge 4 may require reducing n_quad.

Integration scale

By default (link_scale = FALSE), integration is on the response scale \mu = g^{-1}(\eta). Setting link_scale = TRUE integrates the linear predictor \eta = \mathbf{x}^{\top}\boldsymbol{\beta} directly, which is useful when the quantity of interest is the area under the link-transformed surface rather than the response. For the identity link the two coincide.

Examples

## 1-D: integral of fitted sin(t) over [-pi, pi] should be near 0
set.seed(1234)
t <- seq(-pi, pi, length.out = 1000)
y <- sin(t) + rnorm(length(t), 0, 0.01)
fit <- lgspline(t, y, K = 4, opt = FALSE)
integrate(fit, lower = -pi, upper = pi)

## Base R integrate still works as expected
integrate(sin, lower = -pi, upper = pi)

## 2-D: volume under fitted volcano surface
data(volcano)
vlong <- cbind(
  rep(seq_len(nrow(volcano)), ncol(volcano)),
  rep(seq_len(ncol(volcano)), each = nrow(volcano)),
  as.vector(volcano)
)
colnames(vlong) <- c("Length", "Width", "Height")
fit_v <- lgspline(vlong[, 1:2], vlong[, 3], K = 18,
                  include_quadratic_interactions = TRUE, opt = FALSE)
integrate(fit_v, lower = c(1, 1), upper = c(87, 61))

lgspline documentation built on Aug. 5, 2026, 1:10 a.m.