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.w...@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.
______________________________________________
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.
