[R] gamm design matrices

Jarrod Hadfield j.hadfield at ed.ac.uk
Sun Mar 2 18:16:34 CET 2014


Hi All,

I would like to fit a varying coefficient model using MCMCglmm. To do  
this it is possible to reparameterise the smooth terms as a mixed  
model as Simon Wood neatly explains in his book.

Given the original design matrix W and penalisation matrix S, I would  
like to find the fixed (X) and random (Z) design matrices for the  
mixed model parameteristaion. X are the columns of W for which S has  
zero eigenvalues and Z is the random effect design matrix for which  
ZZ' = W*S*^{-1}W* where W* is W with the fixed effect columns dropped  
and S* is S with the fixed effect row/columns dropped.

Having e as the eigenvlaues of S* and L the associated eigenvectors  
then Z  can be formed as W*LE^{-1}  where E = Diag(sqrt(e)).

However, obtaining W and S from mgcv's smooth.construct and applying  
the above logic I can recover the X but not the Z that gamm is passing  
to nlme. I have posted an example below.

I have dredged the mgcv code but to no avail to see why I get  
differences. I want to fit the model to ordinal data, and I also have  
a correlated random effect structure (through a pedigree) hence the  
need to use MCMCglmm.

Any help would be greatly appreciated.

Cheers,

Jarrod




library(mgcv)

x<-rnorm(100)
y<-sin(x)+rnorm(100)

dat<-data.frame(y=y,x=x) # generate some data

smooth.spec.object<-interpret.gam(y~s(x,k=10))$smooth.spec[[1]]

sm<-smooth.construct(smooth.spec.object,data=dat, knots=NULL)

# get penalty matrix S = sm$S[[1]] and original design matrix W=sm$X

    Sed<-eigen(sm$S[[1]])
    Su<-Sed$vectors
    Sd<-Sed$values
    nonzeros <- which(Sd > sqrt(.Machine$double.eps))

    Xn<-sm$X[,-nonzeros]

    Zn<-sm$X%*%Su[,nonzeros]%*%diag(1/sqrt(Sd[nonzeros]))

# form X and Z

    m2.lme<-gamm(y~s(x,k=10), data=dat)

    Xo<-m2.lme$lme$data$X
    Zo<-m2.lme$lme$data$Xr

# fit gamm and extract X and Z used by lme

plot(Xo,Xn)  # OK
plot(Zo,Zn)  # not OK













-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.



More information about the R-help mailing list