[R] gamm design matrices

[Jarrod Hadfield]

Hi Simon,

Sorry, I spoke too soon. I now have the basic spline code working but  
with `by' variables I get an error in smoothCon:

Error in ExtractData(object, data, knots) :
   'names' attribute [1] must be the same length as the vector [0]

  and a segfault with smooth.construct (see example code below). I  
think it is something to do with by.done, but although by.done is  
mentioned as an item/object/attribute in the details section of  
?smoothCon and the value section of ?smooth.construct, there seems to  
be no additional documentation so I'm not sure what by.done is or  
should be.

Any help would be gratefully received.

Cheers,

Jarrod


nt<-10
nid<-100

X<-matrix(rnorm(nt*nid),nid,nt)
beta<-sin(1:nt)
y<-X%*%beta+rnorm(nid)

index<-matrix(rep(1:nt,each=nid),nid, nt)

dat<-data.frame(y=y, X=X, index=index)

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

   sm<-smoothCon(smooth.spec.object,data=dat, knots=NULL,absorb.cons=TRUE)
Error in ExtractData(object, data, knots) :
   'names' attribute [1] must be the same length as the vector [0]

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






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.