[R] mgcv: REML increase with new predictor

[Greg Dropkin]

hi Simon

this follows on from the example where gcv increased unexpectedly with
increasing basis dimension. I'm now looking at t2 tensor splines with
REML, and find that the REML score can increase when adding a new
predictor. Again, this seems odd.

thanks

Greg

library(mgcv)
set.seed(0)
#simulate some data
x1<-runif(2000)
x2<-rnorm(2000)
x3<-rpois(2000,3)
d<-runif(2000)
linp<--5+x1+0.3*x2+4*cos(x1*x2)
y<-rpois(2000,exp(linp))
#y does not depend on d
table(y)
sum(y)

m0<-gam(y~t2(x1)+t2(x2)+t2(x1,x2),family="quasipoisson",method="REML",scale=-1)
sc0<-m0$scale
sc0
m0<-gam(as.formula(m0),family="quasipoisson",method="REML",scale=sc0)
m1<-gam(update.formula(m0,~.+t2(d)),family="quasipoisson",method="REML",scale=sc0)

#dev has decreased slightly, though not significantly
m0$dev
m1$dev

#REML has increased, unexpected (to me)
m0$gcv
m1$gcv

#does this go away if using full=TRUE?

m0t<-gam(y~t2(x1,full=TRUE)+t2(x2,full=TRUE)+t2(x1,x2,full=TRUE),family="quasipoisson",method="REML",scale=-1)
sc0t<-m0t$scale
sc0t
m0t<-gam(as.formula(m0t),family="quasipoisson",method="REML",scale=sc0t)
m1t<-gam(update.formula(m0t,~.+t2(d,full=TRUE)),family="quasipoisson",method="REML",scale=sc0t)

#dev has still decreased slightly
m0t$dev
m1t$dev

#REML has still increased
m0t$gcv
m1t$gcv

#how about te rather than t2?

m0te<-gam(y~te(x1)+te(x2)+te(x1,x2),family="quasipoisson",method="REML",scale=-1)
sc0te<-m0t$scale
sc0te
m0te<-gam(as.formula(m0t),family="quasipoisson",method="REML",scale=sc0t)
m1te<-gam(update.formula(m0t,~.+te(d)),family="quasipoisson",method="REML",scale=sc0t)

#dev has still decreased slightly
m0te$dev
m1te$dev

#now REML decreases
m0te$gcv
m1te$gcv