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.
