[R] gamm design matrices

[Simon Wood]

Dear Jarrod,

I've only managed to have a very quick look at this, but I wonder if the 
difference is to do with identifiability constraints? It looks to me as 
if your smooth.construct call does not have sum to zero constraints 
applied to the smooth, but by default 'gamm' will impose these (even 
when the model has only one smooth). 'gamm' then includes an explicit 
intercept term, which would explain why the model matrix dimensions 
still look ok.

If you mimic gamm by calling smoothCon rather than smooth.construct and 
asking for constraints to be absorbed, and then add the model matrix 
column for the intercept back in, then it should be possible to get the 
model matrices to match.

best,

Simon

On 02/03/14 17:16, Jarrod Hadfield wrote:
> 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
>
>
>
>
>
>
>
>
>
>
>
>
>


-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283