Control Basis Dimension with k

Choose an appropriate basis dimension k and diagnose underfitting with wk.check().

The k argument in s(x, k=...) sets the maximum number of basis functions: the upper limit on how wiggly the smooth can be. The actual wiggliness is then controlled by the penalty. If k is too small, the smooth cannot represent the true shape no matter how the penalty is tuned. Use wk.check() to catch this: a k-index below 1 and a small p-value signal that k needs to be raised.

Load data

import whittaker as wk

data = wk.load_dataset("mcycle", as_frame=True)
data.head()

	times	accel
0	0.0	-11.075042
1	0.2	2.887586
2	0.4	3.909896
3	0.6	14.610615
4	2.2	-0.172051

Fit with k too small

Fit with k=4 (far too few basis functions for the complex shape in mcycle).

# Fit with k too small
model_small = wk.GAM("accel ~ s(times, k=4)").fit(data)

model_small.summary()

GAM fit summary
============================================================
Formula:    accel ~ s(times, k=4)
Family:     Gaussian(link='identity')
Inference:  GCV
Observations: 133
Coefficients: 4

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                -45.6924     2.5368    -18.012    < 1e-16

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(times, k=4)              2.99      3    287.957    < 1e-16

Total EDF:  3.99
Scale est:  855.890185
Deviance:   110416.4241
Null dev:   357878.4929
Dev. expl:  69.1%
GCV score:  882.376746
AIC:        1279.46
BIC:        1291.00

The EDF is close to the maximum allowed by k, which is a warning sign that the smooth is being constrained rather than smoothed.

Check the small-k model

wk.check(model_small)

Look at the k-index in the diagnostics table. A value below 1 means the smooth is running up against the basis dimension limit. The residual pattern in the plots will also show structure that the model cannot capture.

Fit with adequate k

Raise k to 15 to give the smooth enough room.

# Fit with adequate k
model_ok = wk.GAM("accel ~ s(times, k=15)").fit(data)

model_ok.summary()

GAM fit summary
============================================================
Formula:    accel ~ s(times, k=15)
Family:     Gaussian(link='identity')
Inference:  GCV
Observations: 133
Coefficients: 15

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                -45.6924     1.8434    -24.787    < 1e-16

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(times, k=15)             8.63      9    654.054    < 1e-16

Total EDF:  9.63
Scale est:  451.969153
Deviance:   55760.7100
Null dev:   357878.4929
Dev. expl:  84.4%
GCV score:  487.237757
AIC:        1200.18
BIC:        1228.00

The EDF should now be well below k-1, indicating the penalty (not the basis) is the binding constraint. This is the desired state.

Check the adequate-k model

wk.check(model_ok)

The k-index should now be above 1 and the residuals should be much cleaner. Once the k-index is comfortably above 1, further increasing k has little effect on the fitted curve (only on computation time).

Compare AIC

# Compare AIC across models
{"k=4 AIC": model_small.aic, "k=15 AIC": model_ok.aic}

{'k=4 AIC': 1279.464846374149, 'k=15 AIC': 1200.1754802922887}

The lower AIC for the adequate-k model confirms the fit improved. In practice, set k high enough that wk.check() gives a clean k-index, then leave it alone.