Quick start

This is a short tour from raw data to a fitted GAM, a model summary, and predictions. Each step is only a few lines. The rest of the guide covers every piece in depth.

Create some data

We start with a simple nonlinear relationship: y = \sin(x) + \varepsilon, where \varepsilon is Gaussian noise. This is a classic test case for smooth function estimation.

import numpy as np
import whittaker as wk

# Generate 200 observations from a noisy sine curve
rng = np.random.default_rng(23)
n = 200
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + rng.normal(0, 0.3, n)

# Whittaker accepts data as a dictionary of arrays
data = {"x": x, "y": y}

The data is a dictionary mapping column names to 1-D NumPy arrays. Every column referenced in the formula must be present.

Specify and fit the model

A GAM formula looks like an R formula. Wrapping a predictor in s() tells Whittaker to model it as a smooth function. The default basis is a thin plate regression spline (TPRS) with 10 basis functions.

# Create a GAM with a smooth term for x
model = wk.GAM("y ~ s(x)")

# Fit the model (smoothing parameters are selected automatically via REML)
model.fit(data, method="REML")
GAM(y ~ s(x), family=Gaussian(link='identity'), fitted)

Whittaker automatically selects the smoothing parameter \lambda by REML (restricted maximum likelihood). No manual tuning is required. The fit() method returns self, so you can chain calls: wk.GAM("y ~ s(x)").fit(data).

Inspect the summary

# Print a summary of the fitted model
print(model.summary())
GAM fit summary
============================================================
Formula:    y ~ s(x)
Family:     Gaussian(link='identity')
Observations: 200
Coefficients: 10

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                 -0.0160     0.0226     -0.709     0.4789

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(x)                       6.99      7    896.560    < 1e-16

Total EDF:  7.99
Deviance:   19.6215
Null dev:   111.9322
Dev. expl:  82.5%
GCV score:  0.106445
Scale est:  0.102191
AIC:        119.39
BIC:        145.74

The summary reports:

  • Smooth terms: the effective degrees of freedom (EDF) for each smooth. An EDF near 1 means the smooth is approximately linear (higher values indicate more complex curvature).
  • Model fit statistics: deviance explained, GCV score, and scale estimate \hat\phi.

Predict on new data

# Create a fine grid for smooth predictions
x_new = np.linspace(0, 2 * np.pi, 100)
new_data = {"x": x_new}

# Predict on the response scale
preds = model.predict(new_data)

# The result contains fitted values
print(f"Prediction shape: {preds.values.shape}")
print(f"First 5 predictions: {preds.values[:5].round(3)}")
Prediction shape: (100,)
First 5 predictions: [-0.001  0.058  0.118  0.177  0.236]

The predict() method returns a PredictionResult with a .values attribute containing predictions on the response scale (\hat\mu = g^{-1}(\hat\eta). For Gaussian with identity link, this is simply X\hat\beta).

Predictions with standard errors

# Predict with standard errors
preds_se = model.predict(new_data, se=True)

# Standard errors are on the linear predictor scale
print(f"SE shape: {preds_se.se.shape}")
print(f"First 5 SEs: {preds_se.se[:5].round(4)}")
SE shape: (100,)
First 5 SEs: [0.1161 0.1025 0.09   0.0792 0.0706]

Setting se=True additionally computes standard errors from the Bayesian posterior covariance matrix V_\beta = \hat\phi (X^\top W X + \sum_j \lambda_j S_j)^{-1}, where S_j are the penalty matrices and \lambda_j the estimated smoothing parameters.

Visualize the fit

import altair as alt

# Build a DataFrame for plotting
x_plot = np.linspace(0, 2 * np.pi, 200)
preds_plot = model.predict({"x": x_plot}, se=True)

# Compute 95% confidence band on the response scale
z = 1.96
lower = preds_plot.values - z * preds_plot.se
upper = preds_plot.values + z * preds_plot.se

# Observed data points
points = alt.Chart({"values": [{"x": float(x[i]), "y": float(y[i])} for i in range(n)]}).mark_circle(
    size=15, opacity=0.3, color="steelblue"
).encode(
    x=alt.X("x:Q", title="x"),
    y=alt.Y("y:Q", title="y"),
)

# Fitted curve
fit_data = [
    {"x": float(x_plot[i]), "fit": float(preds_plot.values[i]),
     "lower": float(lower[i]), "upper": float(upper[i])}
    for i in range(len(x_plot))
]
line = alt.Chart({"values": fit_data}).mark_line(color="firebrick", strokeWidth=2).encode(
    x="x:Q", y="fit:Q"
)

# Confidence band
band = alt.Chart({"values": fit_data}).mark_area(opacity=0.2, color="firebrick").encode(
    x="x:Q", y="lower:Q", y2="upper:Q"
)

# True function
true_data = [{"x": float(x_plot[i]), "true": float(np.sin(x_plot[i]))} for i in range(len(x_plot))]
true_line = alt.Chart({"values": true_data}).mark_line(
    color="gray", strokeDash=[4, 4], strokeWidth=1.5
).encode(x="x:Q", y="true:Q")

# Combine all layers
(band + points + line + true_line).properties(
    width="container", height=300,
    title="GAM fit: y ~ s(x)"
)

The red curve is the estimated smooth \hat{f}(x), the shaded band is the 95% pointwise confidence interval, the gray dashed line is the true \sin(x), and the blue points are the observed data. The GAM recovers the true shape closely, with the confidence band covering the truth everywhere.

A model with multiple smooths

GAMs shine when you have multiple predictors, each with a potentially nonlinear effect. Here we generate data with two smooth effects and a linear term:

# Generate data with two smooth effects
rng = np.random.default_rng(23)
n = 300
x1 = np.linspace(0, 2 * np.pi, n)
x2 = rng.uniform(0, 1, n)
x3 = rng.normal(0, 1, n)

# True relationship: sin(x1) + x2^2 + 0.5*x3 + noise
y = np.sin(x1) + x2**2 + 0.5 * x3 + rng.normal(0, 0.3, n)

data_multi = {"x1": x1, "x2": x2, "x3": x3, "y": y}
# Fit a GAM with two smooth terms and one linear term
model_multi = wk.GAM("y ~ s(x1) + s(x2) + x3")
model_multi.fit(data_multi, method="REML")

# Print the summary
print(model_multi.summary())
GAM fit summary
============================================================
Formula:    y ~ s(x1) + s(x2) + x3
Family:     Gaussian(link='identity')
Observations: 300
Coefficients: 20

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                  0.3329     0.0166     20.046    < 1e-16
  x3                           0.5087     0.0170     29.865    < 1e-16

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(x1)                      7.57      8   1910.864    < 1e-16
  s(x2)                      2.51      3    251.012    < 1e-16

Total EDF:  12.08
Deviance:   23.8113
Null dev:   265.7658
Dev. expl:  91.0%
GCV score:  0.086171
Scale est:  0.082701
AIC:        115.69
BIC:        160.43
# Predicted vs. observed for the multi-predictor model
preds_multi = model_multi.predict(data_multi)

obs_vs_pred = [
    {"observed": float(y[i]), "predicted": float(preds_multi.values[i])}
    for i in range(len(y))
]

scatter_multi = alt.Chart({"values": obs_vs_pred}).mark_circle(
    size=20, opacity=0.4, color="steelblue"
).encode(
    x=alt.X("observed:Q", title="Observed y"),
    y=alt.Y("predicted:Q", title="Predicted y"),
)

# 1:1 reference line
y_range = [float(min(y.min(), preds_multi.values.min())),
           float(max(y.max(), preds_multi.values.max()))]
ref_line = alt.Chart(
    {"values": [{"v": y_range[0]}, {"v": y_range[1]}]}
).mark_line(
    color="gray", strokeDash=[4, 4], strokeWidth=1.5
).encode(
    x=alt.X("v:Q", title="Observed y"),
    y=alt.Y("v:Q", title="Predicted y"),
)

(ref_line + scatter_multi).properties(
    width="container", height=350,
    title="Multi-predictor GAM: predicted vs. observed"
)

The formula "y ~ s(x1) + s(x2) + x3" specifies:

  • s(x1): a smooth function of x1 (captures the sine wave)
  • s(x2): a smooth function of x2 (captures the quadratic)
  • x3: a plain linear term (enters the model as \beta \cdot x_3)

The summary shows that s(x1) uses more effective degrees of freedom (capturing the sine wave’s curvature) while s(x2) uses fewer (a quadratic is a simpler shape).

Non-Gaussian responses

For count data, binary outcomes, or other non-Gaussian responses, specify a family:

# Poisson count data
rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)
mu = np.exp(0.5 + 0.8 * np.sin(x))  # true log-linear rate
y_counts = rng.poisson(mu)

# Fit a Poisson GAM
model_pois = wk.GAM("y ~ s(x)", family=wk.Poisson())
model_pois.fit({"x": x, "y": y_counts.astype(float)}, method="REML")

# Predictions are on the response scale (counts)
preds_pois = model_pois.predict({"x": x})
print(f"Mean predicted count: {preds_pois.values.mean():.2f}")
print(f"Mean observed count: {y_counts.mean():.2f}")
Mean predicted count: 1.93
Mean observed count: 1.93
# Scatter of observed counts + fitted rate curve + true rate
x_fine = np.linspace(0, 2 * np.pi, 200)
preds_fine = model_pois.predict({"x": x_fine})
true_rate = np.exp(0.5 + 0.8 * np.sin(x_fine))

obs_counts = [
    {"x": float(x[i]), "y": float(y_counts[i])} for i in range(len(x))
]
points_pois = alt.Chart({"values": obs_counts}).mark_circle(
    size=15, opacity=0.3, color="steelblue"
).encode(
    x=alt.X("x:Q", title="x"),
    y=alt.Y("y:Q", title="Count"),
)

fit_rate = [
    {"x": float(x_fine[i]), "rate": float(preds_fine.values[i])}
    for i in range(len(x_fine))
]
fitted_line = alt.Chart({"values": fit_rate}).mark_line(
    color="firebrick", strokeWidth=2
).encode(x="x:Q", y="rate:Q")

true_rate_data = [
    {"x": float(x_fine[i]), "rate": float(true_rate[i])}
    for i in range(len(x_fine))
]
truth_line = alt.Chart({"values": true_rate_data}).mark_line(
    color="gray", strokeDash=[4, 4], strokeWidth=1.5
).encode(x="x:Q", y="rate:Q")

(points_pois + fitted_line + truth_line).properties(
    width="container", height=300,
    title="Poisson GAM: fitted rate vs. true rate"
)

The Poisson family uses a log link, so the model is \log(\mu) = \beta_0 + f(x) and predictions are on the count scale after applying \exp.

Where to go next