[R] mgcv: increasing basis dimension
Greg Dropkin
gregd at gn.apc.org
Tue Feb 14 11:42:56 CET 2012
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)
>
>
More information about the R-help
mailing list