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