Spline Regression › Implementation in Scikit-learn
Machine Learning from ScratchAll chapters Spline Regression › Implementation in Scikit-learnOne section · Sections in this chapter5 / 8 ↓
Reading tools Definitions & inline links
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
Scikit-learn provides excellent tools for implementing spline regression through the SplineTransformer and LinearRegression components. This approach separates the basis function construction from the regression fitting, making the process modular and interpretable. Let's walk through a complete implementation step by step.
Step 1: Data Preparation and Setup
We'll start by importing the necessary libraries and generating synthetic data with a non-linear relationship.
import numpy as np
from sklearn.model_selection import train_test_split
# Set random seed for reproducibility
np.random.seed(42)
# Generate synthetic data with non-linear relationship
n_samples = 200
x = np.linspace(0, 10, n_samples)
y_true = np.sin(x) + 0.1 * x + 0.05 * x**2
noise = np.random.normal(0, 0.3, n_samples)
y = y_true + noise
# Split into training and test sets
X_train, X_test, y_train, y_test = train_test_split(
x.reshape(-1, 1), y, test_size=0.3, random_state=42
)This synthetic dataset combines a sinusoidal pattern with linear and quadratic terms, then adds noise to create a complex non-linear relationship that demonstrates how spline regression behaves.
Step 2: Building the Spline Regression Pipeline
We'll create a pipeline that first transforms the input using B-spline basis functions, then fits a linear regression model on these transformed features.
# Create spline regression pipeline
spline_reg = Pipeline(
[
("spline", SplineTransformer(n_knots=5, degree=3, include_bias=True)),
("linear", LinearRegression()),
]
)
# Train the model
spline_reg.fit(X_train, y_train)
# Make predictions
y_pred_train = spline_reg.predict(X_train)
y_pred_test = spline_reg.predict(X_test)
# Calculate performance metrics
train_mse = mean_squared_error(y_train, y_pred_train)
test_mse = mean_squared_error(y_test, y_pred_test)
train_r2 = r2_score(y_train, y_pred_train)
test_r2 = r2_score(y_test, y_pred_test)
# Perform cross-validation
cv_scores = cross_val_score(
spline_reg, X_train, y_train, cv=5, scoring="neg_mean_squared_error"
)
cv_mse = -cv_scores.mean()
cv_std = cv_scores.std()The pipeline approach is elegant because it encapsulates both the transformation and regression steps, making the model easy to use and maintain.
Step 3: Evaluating Model Performance
Let's examine the model's performance metrics to understand how well it captures the underlying relationship.
=== Spline Regression Performance === Training MSE: 0.0795 Test MSE: 0.0706 Training R²: 0.9792 Test R²: 0.9810 Cross-Validation MSE: 0.0841 (+/- 0.0451)
The test MSE and R² scores indicate strong predictive performance. The R² value close to 1.0 suggests the model explains most of the variance in the target variable. The small difference between training and test metrics indicates good generalization without overfitting. The cross-validation results confirm the model's stability across different data splits.
Step 4: Inspecting Model Components
Let's examine the internal structure of our fitted spline model to understand how it represents the data.
# Extract model components
spline_transformer = spline_reg.named_steps["spline"]
linear_model = spline_reg.named_steps["linear"]
# Get basis function information
n_basis_functions = spline_transformer.transform(X_train).shape[1]
model_coefficients = linear_model.coef_
model_intercept = linear_model.intercept_=== Model Structure === Number of basis functions: 7 Model intercept: 1.5774 First 5 coefficients: [-6.41834866 -0.67433626 0.01269634 -2.25044819 3.96674687] Coefficient range: [-6.4183, 4.0750]
The model uses multiple B-spline basis functions to represent the non-linear relationship. Each coefficient represents the weight assigned to its corresponding basis function. The range of coefficients indicates how different basis functions contribute to the final prediction across different regions of the input space.
Step 5: Comparing Different Spline Configurations
To find the optimal configuration, let's compare different combinations of knot counts and polynomial degrees.
# Define configurations to test
configurations = [
{"n_knots": 3, "degree": 3, "name": "Cubic, 3 knots"},
{"n_knots": 5, "degree": 3, "name": "Cubic, 5 knots"},
{"n_knots": 7, "degree": 3, "name": "Cubic, 7 knots"},
{"n_knots": 5, "degree": 2, "name": "Quadratic, 5 knots"},
{"n_knots": 5, "degree": 4, "name": "Quartic, 5 knots"},
]
results = []
for config in configurations:
# Create and fit model
model = Pipeline(
[
(
"spline",
SplineTransformer(
n_knots=config["n_knots"],
degree=config["degree"],
include_bias=True,
),
),
("linear", LinearRegression()),
]
)
# Cross-validation
cv_scores = cross_val_score(
model, X_train, y_train, cv=5, scoring="neg_mean_squared_error"
)
cv_mse_config = -cv_scores.mean()
# Fit on full training set and test
model.fit(X_train, y_train)
y_pred_config = model.predict(X_test)
test_mse_config = mean_squared_error(y_test, y_pred_config)
test_r2_config = r2_score(y_test, y_pred_config)
n_features = model.named_steps["spline"].transform(X_train).shape[1]
results.append(
{
"config": config["name"],
"cv_mse": cv_mse_config,
"test_mse": test_mse_config,
"test_r2": test_r2_config,
"n_features": n_features,
}
)
# Find best configuration
best_result = min(results, key=lambda x: x["cv_mse"])=== Configuration Comparison === Configuration | CV MSE | Test MSE | Test R² | Features --------------------------------------------------------------------------- Cubic, 3 knots | 0.0950 | 0.0951 | 0.9745 | 5 Cubic, 5 knots | 0.0841 | 0.0706 | 0.9810 | 7 Cubic, 7 knots | 0.0825 | 0.0733 | 0.9803 | 9 Quadratic, 5 knots | 0.1055 | 0.1214 | 0.9674 | 6 Quartic, 5 knots | 0.0819 | 0.0759 | 0.9796 | 8 Best configuration: Quartic, 5 knots (CV MSE = 0.0819)
The comparison reveals important trade-offs between model complexity and performance. More knots generally provide better fit but risk overfitting, while fewer knots may underfit complex patterns. The cubic splines (degree 3) typically offer the best balance between flexibility and smoothness. The best configuration minimizes cross-validation error while maintaining good test performance.
Step 6: Natural Cubic Splines Implementation
For applications requiring stable extrapolation beyond the training data range, we can implement natural cubic splines with custom boundary conditions.
from sklearn.base import BaseEstimator, TransformerMixin
class NaturalCubicSplineTransformer(BaseEstimator, TransformerMixin):
"""Custom transformer for natural cubic splines with zero second derivative at boundaries."""
def __init__(self, n_knots=5):
self.n_knots = n_knots
def fit(self, X, y=None):
# Place knots at evenly spaced quantiles
self.knots_ = np.linspace(X.min(), X.max(), self.n_knots + 2)[1:-1]
return self
def transform(self, X):
n_samples = X.shape[0]
n_basis = self.n_knots + 2
basis_matrix = np.zeros((n_samples, n_basis))
# Linear terms
basis_matrix[:, 0] = 1 # intercept
basis_matrix[:, 1] = X.flatten() # linear term
# Cubic terms with natural boundary conditions
for i, knot in enumerate(self.knots_):
basis_matrix[:, i + 2] = np.maximum(0, X.flatten() - knot) ** 3
return basis_matrix
# Create and train both models
standard_spline = Pipeline(
[
("spline", SplineTransformer(n_knots=5, degree=3, include_bias=True)),
("linear", LinearRegression()),
]
)
natural_spline = Pipeline(
[
("spline", NaturalCubicSplineTransformer(n_knots=5)),
("linear", LinearRegression()),
]
)
standard_spline.fit(X_train, y_train)
natural_spline.fit(X_train, y_train)
# Make predictions
y_pred_standard = standard_spline.predict(X_test)
y_pred_natural = natural_spline.predict(X_test)
# Calculate metrics
mse_standard = mean_squared_error(y_test, y_pred_standard)
mse_natural = mean_squared_error(y_test, y_pred_natural)
r2_standard = r2_score(y_test, y_pred_standard)
r2_natural = r2_score(y_test, y_pred_natural)=== Standard vs Natural Cubic Splines === Standard B-spline - Test MSE: 0.0706, R²: 0.9810 Natural cubic spline - Test MSE: 0.0746, R²: 0.9800 Performance difference: 0.0041 MSE
The comparison between standard and natural cubic splines reveals subtle differences in boundary behavior. Natural cubic splines enforce zero second derivative at the boundaries, making them linear at the edges. This constraint improves extrapolation stability but may sacrifice some fitting flexibility. The small performance difference suggests both approaches work well for this dataset, with the choice depending on whether boundary behavior or interior fit is more important for your application.
Scikit-learn focuses on the more general B-spline basis, which is why we implement natural cubic splines manually. For applications requiring extrapolation or where boundary conditions are critical (such as economic time series or physical measurements), natural cubic splines offer theoretical advantages despite similar in-sample performance.
Key Parameters
Below are the main parameters that control spline regression behavior and performance.
-
n_knots: Number of interior knots (default: 5). More knots increase model flexibility but risk overfitting. Start with 3-5 knots and increase based on cross-validation performance. For datasets with 100-200 observations, 5-7 knots typically work well. -
degree: Degree of the spline polynomials (default: 3 for cubic). Cubic splines (degree=3) provide C² continuity and are most commonly used. Linear splines (degree=1) are simpler but less smooth. Higher degrees (4+) rarely improve performance and increase computational cost. -
include_bias: Whether to include an intercept column in the basis matrix (default: True). Set to True when using with LinearRegression to avoid multicollinearity issues. Set to False if your regression model already includes an intercept term. -
knots: Explicit knot positions (default: None, places knots at quantiles). When None, knots are automatically placed at evenly spaced quantiles of the training data. Specify custom positions when domain knowledge suggests specific break points. -
extrapolation: How to handle predictions beyond the training data range (default: 'constant'). Options include 'constant' (use boundary values), 'linear' (linear extrapolation), 'continue' (extend polynomial), or 'periodic' (wrap around). Choose based on expected behavior outside the training range.
Key Methods
The following methods are used to work with spline regression models in scikit-learn.
-
fit(X, y): Trains the spline regression model by first transforming X into spline basis functions, then fitting the linear regression on these transformed features. Call this once on your training data. -
predict(X): Returns predictions for new data X by transforming it through the spline basis and applying the learned linear coefficients. Use this for making predictions on test data or new observations. -
transform(X): (SplineTransformer only) Converts input features X into the spline basis representation without fitting a regression model. Useful when you want to use the spline features with a different estimator or inspect the basis functions directly. -
score(X, y): Returns the R² score (coefficient of determination) for the model on the given data. Higher values (closer to 1.0) indicate better fit. Use this for quick model evaluation during hyperparameter tuning.