[R] gamm(): nested tensor product smooths

[Fabian Scheipl]

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?

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("ps","ps"),m=c(3,1),k=c(5,5),mp=F))

# Tensor product penalty rank appears to be too low

# 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.
--