[R] gamm(): nested tensor product smooths

[Simon Wood]

> I'd like to compare tests based on the mixed model representation of
> additive models, testing among others
>
> y=f(x1)+f(x2) vs y=f(x1)+f(x2)+f(x1,x2)
> (testing for additivity)
>
> In mixed model representation, where X represents the unpenalized part of
> the spline functions and Z the "wiggly" parts, this would be:
>
> y=X%*%beta+ Z_1%*%b_1+ Z_2%*%b_2
> vs
> y=X%*%beta+ Z_1%*%b_1+ Z_2%*%b_2 + Z_12 %*% b_12
>
> where b are random effect vectors and the hypothesis to be tested is
>
> H_0: Var(b_12)=0 (<=> b_12_i == 0 for all i)
>
> the problem:
> gamm() does not seem to support the use of nested tensor product splines,
> does anybody know how to work around this?

`gamm' doesn't support nested smooths (see `WARNINGS' in ?gamm). Actually it's 
possible that this restriction could be removed now. It was initially 
difficult to handle because of the way the side constraints were handled 
symbolically - this has been changed to a more general scheme to allow side 
conditions to be applied correctly for any user supplied smooth, so it may be 
that it's not to difficult to upgrade gamm - I'll take a look.

>
> example code:  (you'll need to source the P-spline constructor from
> ?p.spline beforehand)
>
> ###########
> test1<-function(x,z,sx=0.3,sz=0.4)
>      { (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
>        0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
>      }
>  n<-400
>  x<-runif(n);z<-runif(n);
>  f <- test1(x,z)
>  y <- f + rnorm(n)*0.1
>
> 
> b<-gam(y~s(x,bs="ps",m=c(3,2))+s(z,bs="ps",m=c(3,2))+te(x,z,bs=c("ps","ps")
>,m=c(3,1),mp=F))
>
> ##works fine
>
> b <-
> gamm(y~s(x,bs="ps",m=c(3,2),k=10))+s(z,bs="ps",m=c(3,2),k=10)+te(x,z,bs=c("
>ps","ps"),m=c(3,1),k=c(5,5)))
>
> # : Singularity in backsolve at level 0, block 1
>
> b <- gamm(y~te(x,z,bs=c("ps","ps"),m=c(3,0),k=c(5,5))) #works fine
>
> # Additionallly, i'd like to work with
> # just one smoothness parameter
> # for the interaction, but mp=F doesn't work either:
>
> b <-
> gamm(y~s(x,bs="ps",m=c(3,2),k=10)+s(z,bs="ps",m=c(3,2),k=10)+te(x,z,bs=c("p
>s","ps"),m=c(3,1),k=c(5,5),mp=F))
>
> # Tensor product penalty rank appears to be too low
- This is a bug - the error's triggered by a buggy sanity check. I'll fix it.

best,
Simon

>
> # which i don't really understand, since the penalty matrix for the
> # p-spline should be
>
> S<-t(diff(diag(5),diff=1))%*%(diff(diag(5),diff=1))
> # penalizing deviations from constant function
> # for the tensor product spline --> m=c(3,1)
> S
> #     [,1] [,2] [,3] [,4] [,5]
> #[1,]    1   -1    0    0    0
> #[2,]   -1    2   -1    0    0
> #[3,]    0   -1    2   -1    0
> #[4,]    0    0   -1    2   -1
> #[5,]    0    0    0   -1    1
>
> #and
> tensor.prod.penalties(list(S,S))
>
> #which should be the penalties for the tensor prod smooth,
> # has rank 20-  that's not enough ?
>
> sum(eigen(S[[2]])$values>1e-10)
> #> [1] 20
>
> ######
>
> I hope some of the experts out there can help me out- any hints would be
> appreciated.... the problem is not caused by the p-splines, it also
> douesn't work with bs="tp".
>
> thank you for your time.
> --
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html and provide commented, minimal,
> self-contained, reproducible code.

-- 
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603  www.maths.bath.ac.uk/~sw283