Skip to content

Specify Z matrix with lmer function

2 messages · Mark Lyman, Spencer Graves

Original
Mark Lyman
# 
Is there a way to specify a Z matrix using the lmer function, where the 
model is written as y = X*Beta + Z*u + e?

I am trying to reproduce smoothing methods illustrated in the paper 
"Smoothing with Mixed Model Software" my Long Ngo and M.P. Wand. 
published in the /Journal of Statistical Software/ in 2004 using the 
lme4 and Matrix packages. The code and data sets used can be found at 
http://www.jstatsoft.org/v09/i01/.

There original code did  not work for me without slight modifications 
here is the code that I used with my modifications noted.

x <- fossil$age
y <- 100000*fossil$strontium.ratio
knots <- seq(94,121,length=25)
n <- length(x)
X <- cbind(rep(1,n),x)
Z <- outer(x,knots,"-")
Z <- Z*(Z>0)
# I had to create the groupedData object with one group to fit the model 
I wanted
grp <- rep(1,n)
grp.dat<-groupedData(y~Z|grp)
fit <- lme(y~-1+X,random=pdIdent(~-1+Z),data=grp.dat)

I would like to know how I could fit this same model using the lmer 
function. Specifically can I specify a Z matrix in the same way as I do 
above in lme?

Thanks,
Mark
16 days later
Spencer Graves
# 
There is probably a way to do what you want, but I don't know how. 
You are to be commended for providing a self-contained example citing an 
interesting article from the Journal of Statistical Software and 
including the modification necessary to make the S-Plus "lme" call work 
in R.

	  I tried several things including the following:

	  fit.r1 <- lmer(y~-1+X+(-1+Z|grp))

	  This seemed to run but gave an answer different from that from lme. 
I also noted that the lmer documentation said the argument "start" was 
"a list of relative precision matrices for the random effects.  This has 
the same form as the slot '"Omega"' in a fitted model.  Only the upper 
triangle of these symmetric matrices should be stored."  This suggested 
I provide "start" with a list consiting of R from the QR decomposition 
of Z:

 > grp.dat.qr <- qr(grp.dat[2:26])
 > grp.dat.R <- qr.R(grp.dat.qr)
 > fit.r <- lmer(y~-1+X+(-1+Z|grp), start=list(grp=grp.dat.R))
Error in lmer(y ~ -1 + X + (-1 + Z | grp), start = list(grp = grp.dat.R)) :
	Leading 1 minor of Omega[[1]] not positive definite

	  Some of the diagonal elements of R were negative, so I just changed 
the sign and tried again:

 > grp.dat.R2 <- (diag(sign(diag(grp.dat.R)))
+                %*% grp.dat.R)
 > fit.r2 <- lmer(y~-1+X+(-1+Z|grp), start=list(grp=grp.dat.R2))
Error in lmer(y ~ -1 + X + (-1 + Z | grp), start = list(grp = 
grp.dat.R2)) :
	Leading 2 minor of Omega[[1]] not positive definite

	  If I were to do more with this, I'd first review all the 
documentation I could find on lmer including Doug Bates, "Fitting linear 
mixed models in R", R News, 5(1):  27-30, May 2005, and "Linear mixed 
model implementation in lmer", July 26, 2005, distributed with lme4 and 
stored on my hard drive under 
"~\R\R-2.2.0\library\lme4\doc\Implementation.pdf".  I might also consult 
Pinheiro and Bates (2000) Mixed-Effects Models in S and S-PLUS 
(Springer), which is my primary source for mixed models generally.  If I 
couldn't figure it out from there, I'd then try to work through the code 
line by line.  Since "lmer" calls "standardGeneric", it's not completely 
obvious how to get the code.  Moreover, "methods" won't get it, because 
"lmer" follows the S4 standard (if I understand correctly). 
Fortunately, 'showMethods("lmer")' produced the following:

Function "lmer":
formula = "formula"

       With this information, I then tried, 'getMethod("lmer", 
"formula")', which gave me the desired source code.  I could then copy 
it into a script file, walk through it line by line, and learn something.

	  hope this helps.
	  spencer graves
Mark Lyman wrote: