[R] help with gls

[Zhenqiang Lu]

Hello R-users,

I am using gls function in R to fit a model with certain correlation
structure.

The medol as:
fit.a<-gls(y~1,data=test.data,correlation=corAR1(form=~1|aa),method="ML")
mu<-summary(fit.a)$coefficient

With the toy data I made to test, the estimate of  mu is exactly equal to
the overall mean of y which can not be true.

But, if I make a toy data with y more than two replicates (for each level
of aa we have more than 2 abservations, for example N=4), the estimates of
mu will not be as the same as the overall mean of y.

Would you please tell me why this happens?
The following is my testing code.
Thanks.
James



require(mvtnorm)
rho=-0.5
N=2
sigma.a=2

Rho.a<-matrix(c(1,rho^(1:(N-1)))[outer(X=1:N,Y=1:N,function(x,y)
1+abs(x-y))],ncol=N)
Sigma.a<-sigma.a^2 * Rho/(1-rho.a^2)
y<-rep(0,N*20)
for (ii in 1:20)
 {y[((ii-1)*N+1):(ii*N)]<-rmvnorm(1, mean=rep(0,N), sigma=Sigma.a)
 }
test.data<-data.frame(y=y,aa=rep(1:20,each=N,len=N*20))
fit.a<-gls(y~1,data=test.data,correlation=corAR1(form=~1|aa),method="ML")
mu.a<-summary(fit.a)$coefficient
rho.a<-getVarCov(fit.a)[1,2]/getVarCov(fit.a)[1,1]
print(c(mean(y),mu.a))

______________________________________________________
Zhenqiang (James) Lu

Department of Statistics
Purdue University
West Lafayette, IN 47906
TEL: (765) 494-0027
FAX: (765) 494-0558