# Inconsistent mgcv gam.check results

I'm not sure if this question is appropriate for cross validated, but I'm not sure where else to post it.

I've built a simple model using the mgcv package.

a <- gam(x ~ s(y), method="REML", data=dat100k)


However, when I run gam.check(a), I get inconsistent output each time I run gam.check(a) with the same model. k is considered OK in the first gam.check(a) call (although, k is close to edf and p is not very large), and on the second call k is considered too low according to the p-value. Is this normal?

Here is output from running gam.check(a) twice:

gam.check(a)

Method: REML   Optimizer: outer newton
full convergence after 8 iterations.
(score 407913 & scale 214.1436).
Hessian positive definite, eigenvalue range [4.001815,49713.04].
Model rank =  10 / 10

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

k'  edf k-index p-value
s(y) 9.00 8.98       1    0.66

> gam.check(a)

Method: REML   Optimizer: outer newton
full convergence after 8 iterations.
(score 407913 & scale 214.1436).
Hessian positive definite, eigenvalue range [4.001815,49713.04].
Model rank =  10 / 10

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

k'  edf k-index p-value
s(y) 9.00 8.98    0.98    0.04 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• What happens if you use set.seed(12345) before each call to gam.check()? – Isabella Ghement Oct 4 '18 at 1:52

The issue is due to the basis dimension test used in gam.check() being based on permutations of model residuals. These permutations are computed using a pseudo random number generator; by design each time you call gam.check() (or directly k.check() itself), a different set of permutations are produced, which subtly alters the p-value of the permutation tests performed.

If you set the seed of the pseudo random number generator before calling gam.check() or k.check(), the results become consistent:

library(mgcv)
set.seed(0)
dat <- gamSim(1,n=200)
b <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat)

set.seed(1)
k.check(b)
set.seed(1)
k.check(b)


which produces:

> set.seed(1)
> k.check(b)
k'      edf   k-index p-value
s(x0)  9 2.318172 0.9959628  0.4950
s(x1)  9 2.305695 0.9693887  0.3225
s(x2)  9 7.654740 0.9605490  0.2825
s(x3)  9 1.232618 1.0372831  0.6850
> set.seed(1)
> k.check(b)
k'      edf   k-index p-value
s(x0)  9 2.318172 0.9959628  0.4950
s(x1)  9 2.305695 0.9693887  0.3225
s(x2)  9 7.654740 0.9605490  0.2825
s(x3)  9 1.232618 1.0372831  0.6850


If you do a greater number of permutations than the default n.rep = 400, then all else equal the p-value should start to stabilise, but this comes at much greater computational cost.

I would also tret the output from this test (certainly the p-value) as just a guide. In your instance, the EDF of the smooth is almost at it's maximum possible value. Regardless of what the p-value of the simple, heuristic test is, I'd want to increase k (double it to k = 20 in the s(...)` call in the formula) and see if the estimated smooth changes at all; does it use a larger EDF than the original fit?