Dear all

I am trying to use lme() to fit a spline mixed effects model with diagonal
random effects covariance matrix G. I would like to compare to obtain the
same results with these 2 alternatives

(1) considering a covariance matrix G, multiple of the identity sigma^2 * I
(2) considering a covariance matrix G, diagonal and equal to: sigma^2 *
D^-1, where D^-1 has 1/d_ii entries (all of them different)

###
## the code would be like this

X<-cbind(1,x)
Id <- factor(rep(1,length(y)))
X=model.matrix(y~X[,2])

# and some Z matrix

Case (1) is easy with

dat  <- groupedData( y ~ X[,-1] |  Id,data = data.frame(y,X,Z))
Z.block<-list(list(Id=pdIdent(~Z-1)))
Z.block<-unlist(Z.block,recursive=FALSE)

fit1 <- lme(y~X[,-1],data=data.fr,random=Z.block,method="REML")

For case (2), now I want to fit the same model, but with a
reparameterization of Z, i.e.

ZZ = Z%*%diag(d)^-0.5

I am trying to define this structure using pdDiag, but I could not get the
same results. Using this:

ZZ.block<-list(list(Id=pdDiag(log((1/d)^.5),~ZZ-1)))

I get

ZZ.block
$Id
Positive definite matrix structure of class pdDiag representing
            [,1]       [,2]       [,3]      [,4]      [,5]      [,6]
[,7]     [,8]   [,9]    [,10]    [,11]
 [1,] 0.06458012 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
0.0000000 0.000000 0.0000 0.000000  0.00000
 [2,] 0.00000000 0.07131829 0.00000000 0.0000000 0.0000000 0.0000000
0.0000000 0.000000 0.0000 0.000000  0.00000
 [3,] 0.00000000 0.00000000 0.08444643 0.0000000 0.0000000 0.0000000
0.0000000 0.000000 0.0000 0.000000  0.00000
 [4,] 0.00000000 0.00000000 0.00000000 0.1078388 0.0000000 0.0000000
0.0000000 0.000000 0.0000 0.000000  0.00000
 [5,] 0.00000000 0.00000000 0.00000000 0.0000000 0.1499137 0.0000000
0.0000000 0.000000 0.0000 0.000000  0.00000
 [6,] 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.2301864
0.0000000 0.000000 0.0000 0.000000  0.00000
 [7,] 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
0.3993904 0.000000 0.0000 0.000000  0.00000
 [8,] 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
0.0000000 0.813031 0.0000 0.000000  0.00000
 [9,] 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
0.0000000 0.000000 2.0765 0.000000  0.00000
[10,] 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
0.0000000 0.000000 0.0000 7.634768  0.00000
[11,] 0.00000000 0.00000000 0.00000000 0.0000000 0.0000000 0.0000000
0.0000000 0.000000 0.0000 0.000000 56.72517

where the diagonal elements are my 1/d_ii

ZZ.block<-unlist(ZZ.block,recursive=FALSE)
dat2 <- groupedData( y ~ X[,-1] | Id,data = data.frame(y,X,ZZ))

fit2 <- lme(y~X[,-1],data=dat2,random=ZZ.block,method="REML")

however fit2, does not give the same results as fit1

Any help on the correct use of pdDiag?

many thanks in advance

Dae-Jin

	[[alternative HTML version deleted]]


