# Control Basis Dimension with k

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


``` python
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`).


``` python
# 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


``` python
wk.check(model_small)
```


<style>
  #altair-viz-a85060b95328457f9d886c7bff5e6682.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-a85060b95328457f9d886c7bff5e6682.vega-embed details,
  #altair-viz-a85060b95328457f9d886c7bff5e6682.vega-embed details summary {
    position: relative;
  }
</style>


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.


``` python
# 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


``` python
wk.check(model_ok)
```


<style>
  #altair-viz-05c5ba1e00344976bbf7e68f89e9468c.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-05c5ba1e00344976bbf7e68f89e9468c.vega-embed details,
  #altair-viz-05c5ba1e00344976bbf7e68f89e9468c.vega-embed details summary {
    position: relative;
  }
</style>


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


``` python
# 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.
