[R] mgcv: increasing basis dimension

[Simon Wood]

Hi Greg,

Recent mgcv versions use extended quasi-likelihood in place of the 
likelihood for (Laplace approx) REML with quasi families (e.g. McCullagh 
and Nelder, GLM book 2nd ed section 9.6): this fixes the problems with 
trying to use the quasi-likelihood directly with REML.

best,
Simon

On 14/02/12 10:42, Greg Dropkin wrote:
> thanks Simon
>
> I'll upgrade R to try t2. The data I'm actually analysing requires scaled
> Poisson so I don't think REML is an option.
>
> thanks
>
> Greg
>
> On 14/02/12 11:22 Simon Wood wrote:
>
> That's interesting. Playing with the example, it doesn't seem to be a
> local minimum. I think that this happens because, although the higher
> rank basis contains the lower rank basis, the penalty can not simply
> suppress all the extra components in the higher rank basis and recover
> exactly what the lower rank basis gave: it's forced to include some of
> the extra stuff, even if heavily penalized, and this is what is
> degrading the higher rank fit in this case.
>
> t2 tensor product smooths seem to be less susceptible to this effect,
> and for reasons I don't understand so does REML based smoothness
> selection (gam(...,method="REML"))
>
> best,
> Simon
>
>
>> hi
>>
>> Using a ts or tprs basis, I expected gcv to decrease when increasing the
>> basis dimension, as I thought this would minimise gcv over a larger
>> subspace. But gcv increased. Here's an example. thanks for any comments.
>>
>> greg
>>
>> #simulate some data
>> set.seed(0)
>> x1<-runif(500)
>> x2<-rnorm(500)
>> x3<-rpois(500,3)
>> d<-runif(500)
>> linp<--1+x1+0.5*x2+0.3*exp(-2*d)*sin(10*d)*x3
>> y<-rpois(500,exp(linp))
>> sum(y)
>>
>> library(mgcv)
>> #basis dimension k=5
>> m1<-gam(y~x1+x2+te(d,bs="ts")+te(x3,bs="ts")+te(d,x3,bs="ts"),family="poisson")
>>
>> #basis dimension k=10
>> m2<-gam(y~x1+x2+te(d,bs="ts",k=10)+te(x3,bs="ts",k=10)+te(d,x3,bs="ts",k=10),family="poisson")
>>
>> #gcv increased
>> m1$gcv
>> m2$gcv
>>
>> summary(m1)
>> summary(m2)
>>
>> gam.check(m1)
>> gam.check(m2)
>>
>>
>> #is this due to bs="ts"?
>>
>> #basis dimension k=5
>> m1tp<-gam(y~x1+x2+te(d,bs="tp")+te(x3,bs="tp")+te(d,x3,bs="tp"),family="poisson")
>>
>> #basis dimension k=10
>> m2tp<-gam(y~x1+x2+te(d,bs="tp",k=10)+te(x3,bs="tp",k=10)+te(d,x3,bs="tp",k=10),family="poisson")
>>
>> m1tp$gcv
>> m2tp$gcv
>>
>> #no
>>
>> summary(m1tp)
>> summary(m2tp)
>>
>> gam.check(m1tp)
>> gam.check(m2tp)
>>
>>
>
>
>
>


-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283