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.
--
______________________________________________
R-help@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.
