[R] gamm design matrices

[Jarrod Hadfield]

Hi Simon,

That's perfect - the matrices now agree.

Thanks,

Jarrod

Quoting Simon Wood <s.wood at bath.ac.uk> on Mon, 03 Mar 2014 13:35:04 +0000:

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



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