# Determining spline basis dimension using Wood's statistical test

In Simon Wood's book Generalized Additive Models (2nd ed.) on page 243, he describes the following procedure for checking that the basis dimension is too small:

Fortunately informal checking that the basis dimension for a particular term is appropriate is quite easy. Suppose that a smooth term is a function of a covariate $$x$$ (which may be vector in general). Let $$\text{nei}(i)$$ denote the set of indices of the $$m$$ nearest neighbours of $$x_i$$ according to an appropriate measure $$\left\|x_i - x_j\right\|$$, and compute the mean absolute or squared difference between the deviance residual $$\epsilon_i$$ and $$\left\{\epsilon_j : j \in \text{nei}(i) \right\}$$. Average this difference over $$i$$ to get a single measure $$\Delta$$ of the difference between residuals and their neighbours. Now repeat the calculation several hundred times for randomly reshuffled residuals. If there is no residual pattern with respect to co variate $$x$$ then $$\Delta$$ should look like an ordinary observation from the distribution of $$\Delta$$ under random re-shuffling, but if the original $$\Delta$$ is unusally small this means that residuals appear positively correlated with their neighbours, suggesting that we may have under-smoothed. IIn that case if the EDF is also close to the basis dimension then it may be appropriate to increase the basis dimension. Note, however, that un-modelled residual auto-correlation and an incorrectly specified mean variance relationship can both lead to the same result, even when the basis dimension is perfectly reasonable. This residual randomization test is easily automated and computationally efficient (given an efficient nearest neighbour algorithm). It is implemented by default in gam.check in mgcv.

I have been using this procedure in mgcv and so far have found it unhelpful. It always seems to reject the null hypothesis, suggesting the basis dimension is too small, and increasing the basis dimension to large numbers leads to rough calculations, minimal change in the resulting fit, and this procedure still rejecting the null hypothesis. Eventually I am including too many basis functions and gam in mgcv says that it cannot be done due to too little data.

Is there any discussion of this procedure that goes beyond a paragraph in Wood's (excellent) book? Should this procedure just not be used? Is my data just insufficient for this type of modeling (I cannot share it or really talk about it) if this procedure is failing?

Generally, for models that are imperfect, which if one takes George Box's statement literally; "Essentially, all models are wrong but some are useful," we would then naturally assume that residuals are, as a rule, structured. That, in the 2D case of fitting a smooth curve to well behaved if noisy data means that the resulting independent variable ranked residuals have a significant polynomial that fits those residuals. We can borrow a term from imaging and call that model misregistration. If we then randomize the independent variable ranking of those residuals, we would decrease the goodness of fit of the polynomial structure discovered prior to randomization, which means that we will then be borrowing from the model's misregistration and converting that into less structured noise. Sometimes, pictures help. So, let's take a pretty good model and show that the residuals are structured in a wavy fashion

Now the difference between a very good model, and one that is not so good is that the relative misregistration is of lesser magnitude compared to the true noise in the data and greater for the not-so-good model. So, to answer the question, it seems to me that the only way to be able to randomize residuals and obtain an insignificant result from the test described, would be to have a very very good model to begin with. Now, I apologize if the description above looks a bit like "hand-waving," to borrow a term favored by physicists, but all I am trying for here is to outline what I think may be relevant, and is certainly not offered as a definitive proof.

• +1 however I can't quite tell whether your answer is intended as a critique of the method or not? Apr 19, 2022 at 9:41
• @GavinSimpson Unfortunately, unstated assumptions are sometimes made that lead to an unrecognized demonstration that those assumptions were made, which demonstration is then misinterpreted to mean something that cannot possibly be shown. I was not critiquing the method as such, I was just wondering out loud if it can be applied to anything usefully, and how one would know that. Personally, when I model something, I assume that the model is incorrect and work to make it "less bad."
– Carl
Apr 19, 2022 at 14:57

It is quite simple to construct a counter example where the test only rejects the null hypothesis with 1 - α probability:

library("gratia")
library("mgcv")

df <- data_sim("eg1")
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")
k.check(m) # gam.check() calls k.check()

r$> k.check(m) k' edf k-index p-value s(x0) 9 3.222496 0.9680923 0.2575 s(x1) 9 2.533375 0.9501894 0.1750 s(x2) 9 7.916005 1.0376331 0.7100 s(x3) 9 1.000894 0.9533817 0.2000  This is a pretty good model from the perspective of the features of this model; smooth functions whose wiggliness has been penalized to minimise the intergrated squared second derivative of the function over the range of the covariate. draw(m)  Here there is a lot more noise than in @Carl's example, so I'm not sure if they would describe this as a "very very good model" or not? :-) In @Carl's informative answer, the blue polynomial is clearly a poor model from the perspective of this test — the fitted polynomial is insufficiently complex enough to capture the non-linear structure in the data. A few more degrees of freedom in the basis would handle that. @Carl's model suggests the presence of autocorrelated noise. In which case it would not surprise me if the test of Wood would "fail". Note the scare quotes: I would suggest this is not a failure of the test (if the data are contaminated by autocorrelated noise) but rather of the model if you assumed independent observations. Smoothing autocorrelated data can be tricky; you could model such a process/data as 1. a very wiggly trend by increasing k, to capture the structure in @Carl's residuals for example, in which case you are saying that you are interested in modelling the time varying mean (level) and take these short runs of high or low values as representing that time-varying mean, 2. by assuming a simple stationary ARMA process (say an AR(1)) for the covariance structure of the residuals/observations. In which case are decomposing the data into a smoothly varying mean (level) plus some effect of previous observations - the autocorrelated noise. Technically 2. has three terms (plus the constant) $$y_t = \beta_0 + f(t) + \gamma_{t-1} + \varepsilon_t$$ where $$\gamma$$ in the AR(1) component and $$\varepsilon$$ is the white noise component. In option 2. we fit those terms separately, either as a fixed effect in the mean component of the model (via an explicit $$\beta_1 y_{t-1}$$ term for an AR(1) say) or via a stationary ARMA process in the residuals or covariance structure of the model). Option 1. says that $$f(t)$$ combines the smooth and the $$\gamma_{t-1}$$ bit into a single term: $$f(t) == s(t) + \gamma_{t-1}$$ because the increased k allows the fitted smooth to be less smooth and model the short-run dependence. It is easy to imagine many ways in which your data may be failing this independence assumption of the test; unmodelled spatial or temporal dependence, one or more missing covariates, etc. But as you can't tell us anything more about your data, one can only speculate. Also note that this test is a heuristic guide at best; it's not intended as some rigid decision maker that tells you to increase k. All it is really doing is telling you that neighbouring residuals are more similar to one another than under randomisation and it is up to you to figure out why that is. I certainly wouldn't be increasing k for smooths that would result in physically (biologically, etc) implausible fitted functions. • Actually, I don't think there is autocorrelation. I say that because the error is actually proportional and the number of waves corresponds to the number of parameters in the model, which I can see in the data by comparing models having different numbers of parameters, thus different degrees of freedom. Interesting, no? – Carl Apr 19, 2022 at 15:03 • BTW, the model I used is a gamma-Pareto type I convolution distribution (GPC), top left of Figure 1 from Comparison of the gamma-Pareto convolution with conventional methods.... As the GPC has no less than two shape parameters, it is not too dissimilar to a spline in that it does a fair job of fitting derivatives, and in addition and unlike a spline does a decent job of extrapolating physically beyond the range of the fit data... – Carl Apr 19, 2022 at 15:48 • cont... Although individual cases have visible structure, on average they don't especially, at least for the GPC model, see Figure 3, where as for worse models without shape parameters per se the residual structure for grouped cases is a lot more obvious, see Figure 4. On the off chance that anyone wants to actually use a GPC model for something, there is accelerated Mathematica code for that as supplementary material here. – Carl Apr 19, 2022 at 16:06 • Finally, one last (I hope) comment. If the noise is large compared to the signal, as you have shown above, the misregistration of the model, compared to the noise may be sufficiently small to allow for use of the test above, but is that a good model? For my purposes, a good model has an R$^2\$ of 0.99+, in which case the misregistration of the model and signal in the data are not negligible and become comparable to the noise magnitude. Indeed, that may be the answer the OP seeks, it may be that the test outlined above is not useful for that type of "good" model. BTW, +1.
– Carl
Apr 19, 2022 at 16:39
• 1) Autocorrelation; I was actually using the time idea to motivate one possible explanation for the test to always reject the null hypothesis. That said, if the structure in your residuals varies with the complexity of your fitted model, doesn't that suggest un-modelled structure or dependence in time? (knowing very little about your area of research I'll refrain from commenting further) Apr 20, 2022 at 8:47