Use a By-Variable Smooth

Fit a separate smooth per group using by= without assuming a shared functional form.

A by-variable smooth fits a distinct smooth for each level of a grouping factor. Unlike a shared smooth, it places no constraint on the shape: one group can have a rising curve while another falls. Use s(x, by=group_indicator) with dummy-coded columns (one binary (0/1) column per group level). This is the right approach when you expect group-specific non-linear shapes rather than a common curve shifted up or down.

Generate data

Simulate 200 observations split between two groups. Group A follows a hump shape that peaks near the center of x. Group B follows an inverted hump that troughs there. Dummy-code the groups into separate binary columns.

import whittaker as wk
import numpy as np

# Simulate x values for two equal-size groups
rng = np.random.default_rng(23)
n = 100
x = rng.uniform(0, 1, size=2 * n)

# Build opposing group-specific response shapes
y = np.concatenate([
    np.sin(np.pi * x[:n]) + rng.normal(0, 0.15, n),       # group A: hump
    -np.sin(np.pi * x[n:]) + rng.normal(0, 0.15, n),      # group B: inverted hump
])

# Dummy-code groups into binary indicator columns
group_a = np.concatenate([np.ones(n), np.zeros(n)])
group_b = np.concatenate([np.zeros(n), np.ones(n)])

data = {"x": x, "y": y, "group_a": group_a, "group_b": group_b}

Group A occupies the first 100 rows; Group B the second. The two shapes are mirror images of each other. A shared smooth would average them into a flat line.

Fit

Fit a model with two by-variable smooths: one for each group. The by=group_a term fits a smooth that is active only where group_a == 1, and similarly for group_b.

# Fit separate smooth per group
model = wk.GAM(
    "y ~ s(x, by=group_a) + s(x, by=group_b)"
).fit(data)

model.summary()

/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/scipy/optimize/_optimize.py:2358: RuntimeWarning: invalid value encountered in scalar subtract
  p = (xf - fulc) * q - (xf - nfc) * r
/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/scipy/optimize/_optimize.py:2359: RuntimeWarning: invalid value encountered in scalar subtract
  q = 2.0 * (q - r)

GAM fit summary
============================================================
Formula:    y ~ s(x, by='group_a') + s(x, by='group_b')
Family:     Gaussian(link='identity')
Inference:  GCV
Observations: 200
Coefficients: 21

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                  0.0723     0.0329      2.199    0.02901

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(x, by='group_a')         1.02      2    311.211    < 1e-16
  s(x, by='group_b')         1.09      2    475.682    < 1e-16

Total EDF:  3.10
Scale est:  0.103564
Deviance:   20.3914
Null dev:   105.4711
Dev. expl:  80.7%
GCV score:  0.105754
AIC:        117.17
BIC:        127.40

Each smooth gets its own EDF row in the summary. Both should have non-trivial EDF values, confirming that the curves are genuinely non-linear and differ between groups.

Partial effects

wk.partial_effects(model)

The partial effects panels show the two smooths on the same x-axis scale. Group A’s smooth should peak near the center while Group B’s troughs there. The opposite shapes the data generation encoded. A shared smooth would have missed this entirely.

Diagnostics

wk.check(model)

Residuals from a correctly specified by-variable model should show no remaining group structure. If you see a pattern tied to the group, consider whether additional terms or a higher k are needed.