[R] help with gls

From: Zhenqiang Lu <lu25_at_stat.purdue.edu>
Date: Sat 06 Jan 2007 - 06:21:29 GMT

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



R-help@stat.math.ethz.ch 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. Received on Sat Jan 06 17:26:27 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sat 06 Jan 2007 - 12:30:26 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.