Dear List:
I am having some trouble reproducing some GLS estimates using matrix
operations that I am not having with other R procedures. Here are some
sample data to see what I am doing along with all code:
mu<-c(100,150,200,250)
Sigma<-matrix(c(400,80,16,3.2,80,400,80,16,16,80,400,80,3.2,16,80,400),n
c=4)
sample.size<-100
temp <-
data.frame(ID=seq(1:sample.size),mvrnorm(n=sample.size,mu,Sigma))
long<-reshape(temp, idvar="ID",
varying=list(c("X1","X2","X3","X4")),
v.names=c("score.1"),direction='long')
long$time<-long$time-1
Now I run this through a gls as follows:
fm1 <- gls(score.1~I(time), data=long,
correlation=corAR1(form=~1|ID), method='ML')
Now, I am trying to replicate the gls using matrix techniques. Here is
my code for to solve (X' V^{-1} X) ^{-1} X' V^{-1} y for the point
estimates and (X' V^{-1} X) ^{-1} for the standard errors.
score<-long$score
X.mat<-model.matrix(score~time, long, row.names=F)
var.mat<-getVarCov(fm1)
I<-diag(sample.size)
V <- kronecker(I,var.mat)
pe<-solve(crossprod(X.mat,solve(V,X.mat)))%*%crossprod(X.mat,solve(V,sco
re))
varcov<-solve(crossprod(X.mat,solve(V,X.mat)))
This perfectly replicates the gls point estimates, but does not
replicate the standard errors. The reason I am concerned is that this
does not occur when I do this for a mixed model or for an OLS model.
For example, using the following code I can perfectly replicate the
point estimates and standard errors from lme()
fm1<-lme(score.1~time, random=~time|ID,long)
V<-extract.lme.cov(fm1,long)
pe<-solve(crossprod(X.mat,solve(V,X.mat)))%*%crossprod(X.mat,solve(V,sco
re))
varcov<-solve(crossprod(X.mat,solve(V,X.mat)))
Also, doing the following perfectly replicates the OLS
fm1<-lm(score.1~time,long)
mse<-var(fm1$residuals)
V<-diag(mse,sample.size*4)
pe<-solve(crossprod(X.mat,solve(V,X.mat)))%*%crossprod(X.mat,solve(V,sco
re))
varcov<-solve(crossprod(X.mat,solve(V,X.mat)))
If you execute this code and compare the matrix results to the internal
functions you will see the OLS and lme results are perfectly replicated.
The gls function is not replicated for the standard errors, but it is
for the point estimates. Is there any reason why the exact same matrix
algebra works for lme and for ols, but not for gls? Theoretically, it
should work for all of them?
Thanks,
Harold
Windows XP
Ver 2.0.1
[[alternative HTML version deleted]]
