Determining spline basis dimension using Wood's statistical test

Question

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?

Answer 1

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.

Answer 2

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.