Regression Splines

Questo metodo per fare regression is not parametersco e it is costruito per risolvere the problema dell’elevata collineartà presente in molti models polinomiali; infatti if ad esempio supponiamo:

\[y = f(x) + \varepsilon\]

Let us consider the following function for \(x\) \[f(x) = sin(2(4x-2)) + 2e^{-(16^2)(x-.5)^2}\]

A high degree (equal to 15) would be needed for polynomial regression to accurately describe the function, but this would lead to high correlation with the following consequences:

• Instability of the estimates
• Sensitivity to outliers

L’idea che sta to the base di this model is that di dividere the spazio da analizzare in more sezioni as divisione it scelgono of the punti detti nodi (knots) attraverso a qualche criterio (solitamente the quantili) e at each sezione it attua a regression.
Using this method as described above, we obtain piecewise functions:

• We define \(K\) knots (internal breakpoints) in the range of values of \(x\): \[\min(x)<\xi_1< \ldots < \xi_K<\max(x)\]

• then we must fit the polynomial model with d degrees in the \(K+1\) intervals

\[(-\infty, \xi_{1}], (\xi_{1}, \xi_{2}], \ldots, (\xi_{K-1}, \xi_{K}], (\xi_{K}, +\infty)\]

Let us therefore consider a general function for \(x\):

\[f(x) = \sum_{j=1}^{p} \beta_j b_j(x)\] where \[b_j(\cdot)\] are functions note chiamate basis functions.

set.seed(123)
x = sort(runif(500))
mu = sin(2*(4*x-2)) + 2*exp(-(16^2)*((x-.5)^2))
y = rnorm(500, mu, .3)
knots = seq(0.1, 0.9, by=.1)
xcut = cut(x, c(-Inf,knots,Inf) )

Step function d=0 \[b_1(X) = I\{x < \xi_1\}, b_2(x) = I\{\xi_1 \leq x < \xi_2\} ,\ldots, b_{K}(x) = I\{\xi_{K-1} \leq x < \xi_K\}, b_{K+1}(x) = I\{ x \geq \xi_K\}\]

plot(x,y ,col="lightgray")
lines(x, mu, col=2, lwd=2)
abline(v=knots, lty=2)
for (i in 1:length(levels(xcut)) ){
sub = (xcut==levels(xcut)[i])
lines(x[sub], fitted(lm(y[sub] ~ 1)) )
}

Piecewise linear regression d=1

\[b_1(x) =I\{x < \xi_1\}, b_2(x) = x \cdot I\{x < \xi_1\}, b_3(x) = I\{\xi_1 \leq x < \xi_2\}, b_4(x) = x \cdot I\{\xi_1 \leq x < \xi_2\},\ldots, b_{2(K+1)}(x) = x \cdot I\{ x \geq \xi_K\}\]

Piecewise polynomial (cubic) regression d=3 \[b_1(x) = 1, b_2(x)= x, b_{3}(x) = x^2, b_{4}(x) = x^3\]

plot(x,y ,col="lightgray")
lines(x, mu, col=2, lwd=2)
abline(v=knots, lty=2)
for (i in 1:length(levels(xcut)) ){
sub = (xcut==levels(xcut)[i])
lines(x[sub], fitted(lm(y[sub] ~ poly(x[sub],3))) )
}

In this first solution, as we can observe, continuity is not guaranteed. To solve this problem, it was devised to impose that at each knot the function has continuous derivatives from degree 0 to degree d-1:
VINCOLI DI CONTINUITA’
\[f(\xi_i) = Y_i \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=1,...,K \] \[f(\xi_i^{+}) = f(\xi_i^{-}) \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=2,...,Z-1\] \[f'(\xi_i^{+}) = f'(\xi_i^{-}) \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=2,...,K-1\] \[f''(\xi_i^{+}) = f''(\xi_i^{-}) \ \ \ \ \ \ \ \ con \ \ \ \ \ \ \ i=2,...,K-1\] \[.\]\[.\]\[.\]

Trade-off between degrees and knots

Below there is an application showing how the two parameters influence the algorithm's behaviour (this application was made with the free version of shinyapp, so I apologise for any inconvenience; I also provide the .Rmd code to facilitate visualisation: codice.Rmd)

This reasoning, while correct, poses two problems:

• Beyond the second derivative, it is difficult to evaluate the continuity of a function.
• The parameters to estimate are \((d+1)(K-1)\) while the constraints are only \(K+(d(K-d+1))\)

Regarding the first problem, a maximum degree of \(d=3\) e it parla di Natural cubic splines; to solve the fact that a unique function is not identified, it is usually imposed that the second derivatives in the extreme intervals are zero.

In definitiva the natural cubic spline with \(K\) knots \(\xi_1,\ldots,\xi_K\) is a piecewise function \(f\) such that:
\(f\) is a polynomial of degree 3 \[[\xi_1,\xi_2],[\xi_2,\xi_3],\ldots,[\xi_{K-1},\xi_{K}]\]
\(f\) is linear in the intervals \((-\infty, \xi_1]\) and \([\xi_K,\infty)\)
\(f\), \(f'\) and \(f''\) are continuous at each knot \(\xi_1,\ldots,\xi_K\): \(f(\xi_k^-)=f(\xi_k^+)\), \(f'(\xi_k^-)=f'(\xi_k^+)\) and \(f''(\xi_k^-)=f''(\xi_k^+)\), \(k=1,\ldots,K\), where \(\xi_k^+\) and \(\xi_k^-\) denote the right and left limit of the function \(f(\cdot)\) at the knot \(\xi_k\)

Number of knots

We can entirely avoid the knot selection problem by formulating a penalised minimisation problem:

\[\underset{f \in \mathcal{F}''} \min \,\, \left\{ \sum_{i=1}^{n}(y_i - f(x_i))^2 + \lambda \int \{f''(t)\}^2 dt \right\}\] where \(f\) ranges over \(\mathcal{F}''\) of twice-differentiable functions

• The first term is RSS, and compares \(f(x_i)\) con \(y_i\) at each \(x_i\)

• Il according to termine is a roughness penalty (rugosità dela function) e controllo as fluttua (wiggly) as it muove \(f(x)\)

• The penalty depends on the tuning parameter \(\lambda\geq 0\):
• \(\lambda=0\) imposes no condition on \(f\) and is computed by interpolating the data

• \(\lambda = \infty\) makes curvatures impossible, and therefore the function \(f\) will be linear

Theorem
To solve this penalised minimisation problem \(f\) among all possible functions of \(\mathcal{F}''\) deve essere a natural cubic spline

• Theorem (Green and Silverman, 1994)

Out of all twice-differentiable functions \(f\), the one that minimizes \[\sum_{i=1}^{n}(y_i - f(x_i))^2 + \lambda \int \{f''(t)\}^2 dt\] is a natural cubic spline with knots at every unique value of \(x_i\)

import numpy as np
from sklearn.preprocessesng import SplineTransformer
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.metrics import mean_squared_error

np.random.seed(42)
n = 200
X = np.sort(np.random.uniform(0, 10, n)).reshape(-1, 1)
y_true = np.sin(X.ravel()) + 0.5 * np.cos(2*X.ravel()) + 0.1 * X.ravel()
y = y_true + np.random.normal(0, 0.3, n)

for n_knots in [3, 6, 10, 15]:
    model = make_pipeline(SplineTransformer(n_knots=n_knots, degree=3), LinearRegression())
    model.fit(X, y)
    y_pred = model.predict(X)
    mse = mean_squared_error(y, y_pred)
    print(f"  n_knots={n_knots:2d} -> MSE: {mse:.4f}, n_params: {model[-1].coef_.shape[0]}")

# Output:
#   n_knots= 3 -> MSE: 0.1892, n_params: 5
#   n_knots= 6 -> MSE: 0.0923, n_params: 8
#   n_knots=10 -> MSE: 0.0854, n_params: 12
#   n_knots=15 -> MSE: 0.0841, n_params: 17

for degree in [1, 2, 3]:
    model = make_pipeline(SplineTransformer(n_knots=8, degree=degree), LinearRegression())
    model.fit(X, y)
    mse = mean_squared_error(y, model.predict(X))
    name = {1: 'Linear', 2: 'Quadratic', 3: 'Cubic'}[degree]
    print(f"  {name} spline (degree={degree}): MSE = {mse:.4f}")

# Output:
#   Linear spline (degree=1): MSE = 0.1054
#   Quadratic spline (degree=2): MSE = 0.0878
#   Cubic spline (degree=3): MSE = 0.0865

st = SplineTransformer(n_knots=5, degree=3)
X_grid = np.linspace(0, 10, 100).reshape(-1, 1)
basis = st.fit_transform(X_grid)
print(f"Number of basis functions: {basis.shape[1]}")
print(f"Each observation is represented by {basis.shape[1]} spline features")
print(f"Basis function values at X=5: {basis[50].round(3)}")

# Output:
# Number of basis functions: 7
# Each observation is represented by 7 spline features
# Basis function values at X=5: [0.    0.    0.167 0.667 0.167 0.    0.   ]