Geometrically Designed Splines: GeDS
Introduction
Geometrically designed spline (GeDS) regression is an efficient and versatile free-knot spline regression technique. This is based on a locally-adaptive knot insertion scheme that produces a piecewise linear spline fit, over which smoother higher order spline fits are subsequently built. GeDS was firstly introduced for the univariate Normal case by Kaishev et al. 2016. It was later expanded to the broader context of Generalized Non-linear Models (GNM) —which include Generalized Linear Models (GLM) as a special case—, and to the bivariate case by Dimitrova et al. 2023. More recently, GeDS methodology has been extended towards the benchmark of Generalized Additive Models (GAM) and Functional Gradient Boosting (FGB) by Dimitrova et al., 2026.
GeDS methodology is implemented on the R package GeDS that I introduce herein. For a tutorial on the the GAM-GeDS methodology, see Generalized Additive Models (GAM) with GeD Splines. For a tutorial focused on FGB-GeDS, see Functional Gradient Boosting (FGB) with GeD Splines.
What are (polynomial) splines?
A polynomial spline is a type of mathematical function made up of several polynomial segments that are “smoothly” joined together at specific points called knots. More specifically, a degree-\(d\) polynomial spline consists of a piecewise degree-\(d\) polynomial function, with continuous \(d-1\) first derivatives at each knot. So, for example, a linear spline is obtained by fitting a straight line in each interval that is defined by the knots, requiring continuity at each knot. Similarly, a cubic spline (degree=3) consists of cubic polynomials joined with continuous first and second derivatives at each knot.
Polynomial splines can accommodate a wide range of shapes and are useful for fitting both simple and complex datasets that may not be well-represented by a single polynomial function. Their flexibility arises from the ability to adjust the position and number of knots, allowing for local adjustments to the fitted curve without affecting the entire function.
Polynomial splines are typically represented through a set of basis functions (known as B-splines) spanning the space of the corresponding spline function. This basis representation —also referred to as the B-spline representation—, constitutes a significant advantage since it notably simplifies the mathematical manipulation and computation of these functions. B-splines provide an efficient way to evaluate, differentiate, and integrate polynomial spline functions, facilitating their use in various applications such as computer graphics and computer-aided design, data interpolation and approximation, and regression analysis.
Normal GeD Spline Regression
Let me start introducing the case of univariate Normal GeDS regression. Consider a response variable \(Y\) and a sole independent variable \(X\), with \(X \in [a, b]\), \(a\in\mathbb{R}\), \(b\in\mathbb{R}\), and assume there is a relationship between \(X\) and \(Y\) of the form: \[\begin{equation} Y = f(X) + \epsilon \end{equation}\] where \(f(\cdot)\) is an unknown function and \(\epsilon\) is a random normal error variable with zero mean and variance \(\mathbb{E}[\epsilon^2]=\sigma_\epsilon^2>0\). In other words, assume \(Y|X\) to be normally distributed. A possible solution to the regression problem of estimating \(f(\cdot)\) based on a sample of observations \(\{X_i, Y_i\}^N_{i=1}\), is to approximate \(f\) with an \(n\)-th order (i.e., degree \(n - 1\)) spline function that takes values from an interval \([a,b]\) and maps them into \(\mathbb{R}\).
Let me now illustrate the usage of the NGeDS() function,
which constructs a Geometrically Designed variable-knot spline
regression model for a response having a Normal distribution. The
implementation supports joint spline components with two or more
covariates. But first, let us simulate some data based on the following
test example. Assume the “true” predictor to be \(\eta=f_1(x)\), where, \[\begin{equation}
f_1(x)=40\frac{x}{1+100x^2}+4,\quad x\in[-2,2].
\end{equation}\]
Based on \(f_1(x)\) I generate random samples \(\{X_i,Y_i\}_{i=1}^N\) with correspondingly normally distributed response variable, \(Y\), and uniformly distributed explanatory variable, \(X\). That is, \(Y_i\sim N(\mu_i,\sigma)\) with \(\sigma=0.2\), \(\mu_i=\eta_i=f_1(X_i)\) and \(X_i\sim U[-2,2]\), \(i=1,...,N\), where \(N\) denotes the sample size. For this example, I will consider \(N=500\). The R implementation of this example is as follows:
# Generate a data sample for the response variable
# Y and the single covariate X
set.seed(123)
N <- 500
f_1 <- function(x) (10*x/(1+100*x^2))*4+4
X <- sort(runif(N, min = -2, max = 2))
# Specify a model for the mean of Y to include only a component
# non-linear in X, defined by the function f_1
means <- f_1(X)
# Add (Normal) noise to the mean of Y
Y <- rnorm(N, means, sd = 0.2)Let us plot our simulated data:
A Normal GeDS model is then fit using the NGeDS()
function. This implements the GeDS regression technique assuming a
Normal response variable:
# Fit a Normal GeDS regression using NGeDS
(Gmod <- NGeDS(Y ~ f(X), beta = 0.6, phi = 0.995, Xextr = c(-2,2)))##
## Call:
## NGeDS(formula = Y ~ f(X), beta = 0.6, phi = 0.995, Xextr = c(-2,
## 2))
##
## Number of internal knots of the second order (linear) spline: 12
## Deviances:
## Order 2 Order 3 Order 4
## 19.65 19.82 19.57where beta and phi are tuning parameters
and Xextr stands for the left-most and right-most limits of
the interval embedding the observations of \(X\). GeDS algorithm operates through two
main stages, A and B, that I will now
briefly explain.
Stage A
GeDS regression starts by fitting a least-squares (LS) straight line
model to the data. After this, the residuals are computed and grouped
into clusters based on their sign. Specifically, the first \(d_1\) consecutive residuals with the same
sign form the first cluster, the next \(d_2\) residuals with the same sign form the
second cluster, and so on. Weights for these clusters of residuals are
then calculated. In the calculation of these weights, the parameter
beta determines the weight put on the within-cluster
mean residual value, \(m_j\), and
the weight put on the within-cluster \(X\)-range, \(\eta_j\). The weights, \(j\), for each cluster of residuals are
given by: \[\begin{equation}
w_j = \beta \times m'_j + (1 − \beta)\times\eta'_j,\quad j =
1,..., u
\end{equation}\] where \(u\)
denotes the total number of clusters of residuals of identical signs,
\(m'_j=m'_j/m_{max}\) is the
normalized within-cluster mean residual value and \(\eta'_j=\eta_j/\eta_{max}\) stands for
the normalized within-cluster range of the independent variable \(X\). A new knot is then inserted in the
cluster of residuals that does not already contain a knot and that has
the maximal weight \(w_j\). The
specific location of the new knot is determined by a weighted average of
the \(X\)-coordinates of the residuals
in the selected cluster of residuals. For a detailed explanation see Kaishev et al. 2016.
A LS linear spline fit is then computed including this new knot.
Before continuing to a new GeDS iteration a stopping rule consisting on
a ratio of consecutive residual sum of squares (RSS) is evaluated: \[\begin{equation}
\text{RSS}(\kappa+q)/\text{RSS}(\kappa)
=\sum_{i=1}^N(Y_i-\hat{f}\left(\delta_{\kappa+q,2};X_i)\right)^2
\Bigg/ \sum_{j=1}^N\left(Y_i-\hat{f}(\delta_{\kappa,2};X_i)\right)^2
\geq \phi_{\text{exit}}
\end{equation}\] where \(q \geq
1\) and \(\phi_{\text{exit}}\in(0,1)\) is a certain
pre-specified threshold level close to one. \(\boldsymbol{\delta}_{\kappa,2}=\{\delta_1<
\delta_2<
\delta_3<...<\delta_{\kappa+2}<\delta_{\kappa+3}=\delta_{\kappa+4}\}\)
denotes the knots locations, and \(\kappa\) the number of internal knots of
the linear GeDS fit. In particular, if the ratio of the RSS of the
current iteration and the RSS of the qth previous iteration
is less than phi, we continue iterating, starting by
fitting a LS linear spline model to the data, with one extra knot.
Otherwise, stage A iterations stop.
Stage A fitting process can be visualized as follows:
layout(matrix(c(1,2), nrow=1, byrow=TRUE))
plot(Gmod, n = 2, which = 1:(Gmod$Nintknots+1), legend.pos = NA, pch = 20, col = "darkgrey")
Argument which in plot.GeDS allows us to
determine the iterations of stage A for which the corresponding GeDS
fits should be plotted; n specifies the order of the GeDS
fit that should be plotted; n = 2, means the linear fit is
to be plotted. Indeed, stage A pursues a “stick breaking” procedure
starting from a straight line fit, that is subsequently broken into
smaller linear pieces.
Note that the final stage A GeDS fit is the fit from the \((\kappa + 1)\)th iteration, where \(\kappa\) denotes the number of internal knots of this fit. This is because the first iteration fit is just a straight line LS fit and then at each GeDS iteration one internal knot is added.
Stage B
Stage B, begins computing the averaging knot location1 of the knots vector of the final linear fit obtained via Stage A’s iterations. This is computed for quadratic, and cubic orders; linear averaging knot location coincides with Stage A’s knot location. Secondly, linear, quadratic, and cubic LS spline fits to the data are computed utilizing the corresponding relocated set of knots, thereby completing the GeDS procedure. Note that the final stage A GeDS fit coincides with the final linear GeDS fit.
The final quadratic GeDS fit can be visualized as follows:
where the black line depicts the real underlying function \(f_1(x)\).
Generalized GeD Spline Regression
Normal GeDS regression methodology can be extended to the wider context of GNM/GLM models (Dimitrova et al. 2023).
The GNM (GLM) fitting problem
Take \(Y\) to be a response variable that depends on a single covariate \(X\). In a generalized model, we assume that \(Y |X\) has a distribution belonging to the exponential family. That is, \(Y |X\) has a probability density function (continuous case) or probability mass function (discrete case) of the form: \[\begin{equation} f(y;\vartheta,\phi) = \exp\left(\frac{y\vartheta − b(\vartheta)}{a(\phi)}+ c(y,\phi)\right) \end{equation}\] where \(\vartheta = \vartheta(μ)\) is the so-called natural parameter and \(\phi\) is the so-called scale parameter; \(a\), \(b\) and \(c\) are functions that characterize the distribution. For example, in the Normal case, \(\vartheta = \vartheta(\mu)=\mu\), \(\phi=\sigma^2\), \(a=a(\phi)=\phi\), \(b=b(\vartheta)=\frac{\vartheta^2}{2}=\frac{\mu^2}{2}\) and \(c=c(y,\phi)=-\frac{1}{2}\left[\frac{y^2}{\sigma^2}+\ln(2\pi\sigma^2)\right]\), and thus, the corresponding density function is: \[\begin{equation} f(y;\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y-\mu)^2}{2\sigma^2}\right). \end{equation}\] The objective is to model the expectation of \(Y|X\): \[\begin{equation} E\left[Y|X\right]=\mu=g^{-1}(\eta(\boldsymbol{\theta};X)) \end{equation}\] where \(g\) is some link function that describes how \(\mu\) is connected with the predictor component \(\eta(\boldsymbol{\theta};X)\). A generalized model is therefore defined by three components:
- a random component, that specifies a distribution for \(Y |X\);
- a systematic component, that relates \(\eta\) to the predictor(s), \(X\);
- and a link function, \(g\), that connects the random (the conditional mean \(E\left[Y|X\right])\) and systematic components (the linear predictor \(\eta(\boldsymbol{\theta};X)\)).
Now, we can assume the predictor \(\eta(\boldsymbol{\theta}; X)\) to be a polynomial spline function, whose degree, number and location of knots are unknown parameters that we need to estimate. Let us denote by \(S_{\boldsymbol{t}_{k,n}}\) the linear space of all \(n\)th order spline functions defined on a set of non-decreasing knots \(\boldsymbol{t}_{k,n} = \{t_i\}^{2n+k}_{i=1}\) , where \(t_n = a\), \(t_{n+k+1} = b\). We use splines with simple knots, except for the \(n\) left and right most knots which will be assumed coalescent, i.e. \[\begin{equation} t_{k,n} = \{t_1 = ... = t_n < t_{n+1} < ... < t_{n+k} < t_{n+k+1} = ... = t_{2n+k}\} \end{equation}\] Then, \(\eta(\boldsymbol{\theta}; X)\) can be expressed as: \[\begin{equation} \eta(\boldsymbol{\theta};X)= f(\boldsymbol{t}_{k,n}; X)=\boldsymbol{\theta}^T\boldsymbol{N}_n(X)=\sum_{i=1}^p\theta_iN_{i,n}(X). \end{equation}\] where \(f \in S_{\boldsymbol{t}_{k,n}}\), \(\left(\boldsymbol{\theta} =\theta_1, ..., \theta_p\right)^T\) is a vector of real valued coefficients and \(N_n(x) = \left(N_{1,n}(X), . . ., N_{p,n}(X)\right)^T\), \(p = n + k\), are B-splines of order \(n\), defined on \(t_{k,n}\).
Generalized GeDS estimation
Now, to estimate the unknown linear predictor function \(\eta(\boldsymbol{\theta};X)\) we can use the generalized GeDS method. Indeed, GeDS methodology successfully extends to the more general Exponential Family of distributions in the context of Generalized Non-Linear/Linear Models. Specifically, stage A of Normal GeDS is extended to the more general GNM/GLM context by replacing LS fitting by Iteratively Reweighted Least Squares (IRLS) fitting and the stopping rule by a deviance-based (instead of RSS) stopping criterion.
The usage of function GGeDS() is illustrated with the
same example as before. Code implemenations and plots are presented as
follows.
Example 1. The “true” linear predictor is assumed to be \(\eta = f_1(x)\), where \[\begin{equation} f_1(x) = 40\frac{x}{1 + 100x^2} + 4,\quad x \in [−2, 2]. \end{equation}\]
Random samples, \(\{X_i,Y_i\}^N_{i=1}\), are generated with correspondingly Poisson, Gamma and Binomially distributed response variable, \(Y\), and uniformly distributed explanatory variable, \(X\). That is, \(Y_i ∼ \text{Poisson}(\mu_i)\), with \(\mu_i = \exp(\eta_i)\) and \(\eta_i = f_1(X_i)\); \(Y_i ∼ \text{Gamma}(\mu_i, \phi)\) with \(\phi = 0.1\), \(\mu_i = \exp(\eta_i)\) and \(\eta_i = f_1(X_i)\); \(Y_i ∼ \text{Binomial}(m,\mu_i)\) with \(m = 50\), \(\mu_i = \exp(\eta_i)/(1 + \exp(\eta_i))\), \(\eta_i = f_1(X_i) − 4\); and \(X_i \sim U[−2, 2]\), \(i = 1, ..., N\), where \(N=500\) is the sample size.
# Generate a data sample for the response variable Y and the covariate X
# See section 4.1 in Dimitrova et al. (2023)
set.seed(123)
N <- 500
f_1 <- function(x) (10*x/(1+100*x^2))*4+4
X <- sort(runif(N, min = -2, max = 2))
# Specify a model for the mean of Y to include only a component
# non-linear in X, defined by the function f_1
means <- exp(f_1(X))Poisson:
#########################################
## (A) Y ~ Poisson + log link function ##
#########################################
# Generate Poisson distributed Y according to the mean model
Y <- rpois(N, means)
# Fit a Poisson GeDS regression using GGeDS
Gmod <- GGeDS(Y ~ f(X), beta = 0.2, phi = 0.995, family = "poisson",
Xextr = c(-2,2))
plot(Gmod, f = function(x) exp(f_1(x)), n = 3, pch = 20, col = "darkgrey")
Gamma:
#######################################
## (B) Y ~ Gamma + log link function ##
#######################################
# Generate Gamma distributed Y according to the mean model
Y <- rgamma(N, shape = means, rate = 0.1)
# Fit a Gamma GeDS regression using GGeDS
Gmod <- GGeDS(Y ~ f(X), beta = 0.1, phi = 0.995, family = Gamma(log),
Xextr = c(-2,2))
plot(Gmod, f = function(x) exp(f_1(x))/0.1, n = 3, pch = 20, col = "darkgrey")
(Quasi)Binomial:
While the binomial distribution is used for modeling binary data with fixed success probabilities (e.g., success/failure, yes/no, 0/1) and independent trials, the quasibinomial distribution accommodates scenarios where the probability of success may vary across trials, and does not strictly adhere to a fixed probability of success for each trial. In particular, it is used for binary outcome data that exhibit so-called overdispersion, which occurs when the observed variance is greater than what the binomial distribution predicts. This variance flexibility provided by the quasibinomial distribution allows for a more adaptable approach to modeling binary outcome data.
############################################
## (C) Y ~ Binomial + logit link function ##
############################################
# Generate Binomial distributed Y according to the mean model
eta <- f_1(X) - 4
means <- exp(eta)/(1+exp(eta))
Y <- rbinom(N, size = 50, prob = means) / 50
# Fit a Binomial GeDS regression using GGeDS
Gmod <- GGeDS(Y ~ f(X), beta = 0.1, phi = 0.995, family = "binomial",
Xextr = c(-2,2))
plot(Gmod, f = function(x) exp(f_1(x) - 4)/(1 + exp(f_1(x) - 4)),
n = 3, pch = 20, col = "darkgrey")
Bivariate GeD Spline Regression
Now, take the case where \(Z\) to be a response variable that depends on two covariates \(X_1=X\) and \(X_2=Y\). Dimitrova et al. 2023, introduce a bivariate extension of GeDS by assuming that the predictor component of the GLM is in the form of bivariate spline function: \[\begin{equation} \eta(\boldsymbol{\theta};\boldsymbol{X})=f(\boldsymbol{T}_{k_1\times k_2};\boldsymbol{X}) \end{equation}\] where \(\boldsymbol{T}_{k_1\times k_2}=\boldsymbol{t}_{1;k_1,n_1}\times\boldsymbol{t}_{1;k_1,n_1}\) and \(\boldsymbol{t}*{1;k_1,n_1},* \boldsymbol{t}{2;k_2,n_2}\) are sets of knots with respect to \(X_1=X\) and \(X_2=Y\), where \(k_1\) and \(k_2\) denote respectively the number internal knots.
Similar to the univariate case, bivariate GeDS regression starts applying the IRLS procedure to find a bivariate linear spline fit \(\hat{f}(\Delta_{k_1\times k_2},\hat{\alpha}_p; x)\) with zero internal knots (i.e. a bivariate linear model). Subsequently, a stopping rule based on the deviances at consecutive iterations is checked and, if fulfilled, we move to stage B. Otherwise, a new knot in the \(X_1=X\) or \(X_2=Y\) direction is added. This is done by identifying a potential new knot \(\delta^*\) for each direction, \(X_1=X\) and \(X_2=Y\), based on the univariate GeDS procedure (GeDS). Then, the knot with maximal within-cluster mean residual/within-cluster range-weight is added.
Below, the implementation of Normal bivariate GeDS and Poisson
bivariate GeDS, with NGeDS() and GGeDS()
functions, respectively, is illustrated. In particular the following
example is considered:
Both functions use the smoothed-ratio (SR) stopping rule by default
for joint smoothers. If the selected Normal bivariate linear fit is rank
deficient, NGeDS() uses the latest full-rank Stage A
fit.
Example 2. The “true” linear predictor is assumed to be \(\eta = f(x_1,x_2)\), where \[\begin{equation} f(x_1,x_2) = \sin(2x_1)\sin(2x_2),\quad x \in [0, 3]^2. \end{equation}\]
For the Normal case, \(Z_i = f(X_{i1}, X_{i2}) + \epsilon_i\), where \(\epsilon_i \sim \mathcal{N}(0, 0.2^2)\). For the Poisson case, \(Z_i \sim \text{Poisson}(\mu_i)\), with \(\mu_i = \exp(\sin(2X_{i1}) + \sin(2X_{i2}))\). In both cases, 400 covariate pairs are sampled uniformly from \([0,3]^2\).
Normal Bivariate GeDS:
# bivariate example
# See Dimitrova et al. (2023), section 5
# Generate a data sample for the response variable
# Z and the covariates X and Y assuming Normal noise
set.seed(123)
doublesin <- function(x){
sin(2*x[,1])*sin(2*x[,2])
}
X <- (round(runif(400, min = 0, max = 3),2))
Y <- (round(runif(400, min = 0, max = 3),2))
Z <- doublesin(cbind(X,Y))
Z <- Z + rnorm(400, 0, sd = 0.2)
# Fit a two dimensional GeDS model using NGeDS
(BivGeDS <- NGeDS(Z ~ f(X, Y), beta = 0.3, phi = 0.95,
Xextr = c(0, 3), Yextr = c(0, 3)))##
## Call:
## NGeDS(formula = Z ~ f(X, Y), beta = 0.3, phi = 0.95, Xextr = c(0,
## 3), Yextr = c(0, 3))
##
## No internal knots found## N1 N2 N3 N4 N5 N6
## 0.72279322 0.69203912 -1.30082252 0.22337793 -0.09707304 -0.39127430
## N7 N8 N9
## -0.64017610 -0.65277821 1.57279517## $Xk
## [1] 0 0 0 3 3 3
##
## $Yk
## [1] 0 0 0 3 3 3## [1] 73.62378# RSS w.r.t true function
f_XY <- apply(cbind(X, Y), 1, function(row) doublesin(matrix(row, ncol = 2)))
mean((f_XY - predict(BivGeDS, n = 3))^2)## [1] 0.1437639
Poisson Bivariate GeDS:
# bivariate example
# See Dimitrova et al. (2023), section 5
# Generate a data sample for the response variable
# Z and the covariates X and Y assuming Poisson distributed error
set.seed(123)
doublesin <- function(x) {
# Adjusting the output to ensure it's positive
exp(sin(2*x[,1]) + sin(2*x[,2]))
}
X <- round(runif(400, min = 0, max = 3), 2)
Y <- round(runif(400, min = 0, max = 3), 2)
# Calculate lambda for Poisson distribution
lambda <- doublesin(cbind(X,Y))
# Generate Z from Poisson distribution
Z <- rpois(400, lambda)
data <- data.frame(X, Y, Z)
# Fit a Poisson GeDS regression using GGeDS
(BivGeDS <- GGeDS(Z ~ f(X,Y), beta = 0.2, phi = 0.99, family = "poisson"))##
## Call:
## GGeDS(formula = Z ~ f(X, Y), family = "poisson", beta = 0.2,
## phi = 0.99)
##
## Number of internal knots of the second order (linear) spline in the X direction: 3
## Number of internal knots of the second order (linear) spline in the Y direction: 4
##
##
## Deviances:
## Order 2 Order 3 Order 4
## 374.3 351.3 350.4## N1 N2 N3 N4 N5 N6
## -0.60538119 1.67124804 1.31493910 0.79205177 -3.57274577 0.59411908
## N7 N8 N9 N10 N11 N12
## 2.16932726 1.29156236 2.60737056 1.04359868 0.13737448 -0.31360836
## N13 N14 N15 N16 N17 N18
## -0.55663694 1.55960245 0.84841499 1.51649745 -2.22675456 1.45766755
## N19 N20 N21 N22 N23 N24
## -0.58546764 -4.41447850 0.43636464 -3.17900717 -1.23871511 -6.75072786
## N25 N26 N27 N28 N29 N30
## -1.81599303 1.06415416 1.60856604 -0.08678098 -1.93137103 -1.03636703## $Xk
## [1] 0.000000 0.000000 0.000000 1.076578 2.100068 3.000000 3.000000 3.000000
##
## $Yk
## [1] 0.0000000 0.0000000 0.0000000 0.4752607 1.1527838 1.6695686 2.9900000
## [8] 2.9900000 2.9900000## [1] 351.2995# Poisson deviance w.r.t true function
f_XY <- apply(cbind(X, Y), 1, function(row) doublesin(matrix(row, ncol = 2)))
mean(poisson()$dev.resids(lambda, predict(BivGeDS, n = 3), wt = 1))## [1] 0.09631082
Multivariate Normal GeDS Regression
The Normal GeDS implementation can also fit one joint tensor-product
spline in more than two covariates. For example, a three-dimensional fit
is specified as NGeDS(y ~ f(x0, x1, x2), data = data). This
is structurally different from the additive model
NGeDSgam(y ~ f(x0) + f(x1) + f(x2), data = data): the
former can represent interactions among all three covariates, whereas
the latter estimates three separate one-dimensional effects.
The distinction can be seen clearly with a pure three-way interaction,
\[\begin{equation} f(x_0,x_1,x_2)=\sin(2\pi x_0)\sin(2\pi x_1)\sin(2\pi x_2), \qquad (x_0,x_1,x_2)\in[0,1]^3. \end{equation}\]
For independent Uniform\((0,1)\) covariates, every one-dimensional marginal effect is zero. Consequently, an additive model cannot recover this signal, while a joint three-dimensional spline can.
set.seed(123)
N <- 800L
interaction_function <- function(x0, x1, x2) {
sin(2 * pi * x0) * sin(2 * pi * x1) * sin(2 * pi * x2)
}
x0 <- runif(N)
x1 <- runif(N)
x2 <- runif(N)
truth <- interaction_function(x0, x1, x2)
y <- rnorm(N, mean = truth, sd = 0.1)
interaction_data <- data.frame(y, x0, x1, x2)
additive_fit <- NGeDSgam(
y ~ f(x0) + f(x1) + f(x2),
data = interaction_data
)## f(x0) has less than 2 linear internal knots.
## f(x1) has less than 2 linear internal knots.
## f(x2) has less than 2 linear internal knots. Quadratic averaging knot location is not computed for base learners with less than 2 internal knots.
## f(x0) has less than 3 linear internal knots.
## f(x1) has less than 3 linear internal knots.
## f(x2) has less than 3 linear internal knots. Cubic averaging knot location is not computed for base learners with less than 3 internal knots.The fits can be compared on new covariate values using the known generating function.
set.seed(456)
N_test <- 5000L
newdata <- data.frame(
x0 = runif(N_test, min(x0), max(x0)),
x1 = runif(N_test, min(x1), max(x1)),
x2 = runif(N_test, min(x2), max(x2))
)
test_truth <- with(
newdata,
interaction_function(x0, x1, x2)
)
orders <- 2:4
additive_mse <- vapply(
orders,
function(n) mean((test_truth - predict(additive_fit, newdata, n = n))^2),
numeric(1)
)
joint_mse <- vapply(
orders,
function(n) mean((test_truth - predict(joint_fit, newdata, n = n))^2),
numeric(1)
)
data.frame(
order = orders,
additive_MSE = additive_mse,
joint_MSE = joint_mse
)## order additive_MSE joint_MSE
## 1 2 0.1269506 0.076898070
## 2 3 0.1272846 0.085200313
## 3 4 0.1280459 0.004199663## x0 x1 x2
## 2 1 1## $x0
## NULL
##
## $x1
## NULL
##
## $x2
## NULL## [1] 10.04737## B[x0=1,x1=1,x2=1] B[x0=1,x1=1,x2=2] B[x0=1,x1=1,x2=3] B[x0=1,x1=1,x2=4]
## 0.4517842 -0.8566270 0.9864489 -0.1015185
## B[x0=1,x1=2,x2=1] B[x0=1,x1=2,x2=2]
## -0.2556598 -2.2734181##
## Experimental multivariate Normal GeDS fit
##
## Call:
## NGeDS(formula = y ~ f(x0, x1, x2), data = interaction_data)
##
## Dimensions (3): x0, x1, x2
## Linear internal knots: x0=2, x1=1, x2=1
##
## Residual sum of squares:
## Order 2 (linear) Order 3 (quadratic) Order 4 (cubic)
## 59.33 65.20 10.05
## Least-squares solver: Order 2=dense, Order 3=dense, Order 4=dense
##
## Stage A:
## Stopping rule: SR (phi = 0.99, q = 2)
## Knot-placement beta: 0.5
## Fitted iterations: 7
## Selected iteration: 5
## Maximum coefficients: 100000
## Largest fitted tensor basis: 60 coefficients
## Termination: SR threshold reachedIn this example the cubic joint tensor fit captures the interaction
especially well, while the additive errors remain close to the variance
of the omitted signal. The selected linear internal knots can be
inspected through joint_fit$linear.intknots, and
predictions for orders 2, 3, and 4 are available through
predict(joint_fit, newdata, n = ...).
Because the tensor basis grows multiplicatively across dimensions, an internal safeguard stops Stage A before an excessively large tensor would be fitted. If only a higher-order Stage B tensor is excessive, that order is omitted with a warning while the smaller available fits are retained.
Joint smoothers above two dimensions are also available in
GGeDS(), NGeDSgam(), and
NGeDSboost(); see the higher-dimensional vignette. For
NGeDS() joint smoothers above two dimensions, parametric
terms, offsets, initial knots or fits, and the LR stopping rule are not
yet supported. The dedicated bivariate fitters remain in use, including
for generalized responses.
When the observations form a complete Cartesian mesh and have equal
weights, the multivariate fitter uses the tensor factorization of the
least-squares operator described by Kaishev
1984. The coefficient tensor is computed by successive univariate
QR solves, without constructing the full tensor-product design matrix.
If the data are scattered, the weights are non-uniform, or a univariate
mesh design is rank deficient, the fitter automatically uses the general
dense least-squares calculation instead. The solver selected for each
spline order is reported by print() and
summary().
See Kaishev et al. 2016 for details.↩︎