2
$\begingroup$

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.
Gradient range [-0.04365533,0.04277738]
(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.
Gradient range [-0.04365533,0.04277738]
(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
$\endgroup$
  • 2
    $\begingroup$ What happens if you use set.seed(12345) before each call to gam.check()? $\endgroup$ – Isabella Ghement Oct 4 '18 at 1:52
2
$\begingroup$

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?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.