On Sat, 2007-01-06 at 01:21 -0500, Zhenqiang Lu wrote:> 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.
>
On the contrary, this has to be true.
Due to a well-known result (a recent is is Greene, 2003, sec 14.2.2, pp
343-344) when the covariates are the same in every equation, a
"seemingly unrelated regression" of this sort has GLS and OLS produce
identical results. The OLS estimate is exactly the arithmetic
average.
@Book{greene2003,
Author = {William H Greene},
Title = {Econometric Analysis},
Publisher = {Prentice-Hall International Ltd},
Address = {London},
Edition = {Fifth},
year = 2003,
}
> 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.
This is what is strange. You should probably provide example code of
this too.
>
> Would you please tell me why this happens?
> The following is my testing code.
The code had some typos in it (and you did not load nlme). I believe
this would be correct:
require(mvtnorm)
library(nlme)
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.a/(1-rho^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))
Regards,
Markus
> 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 at 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.
>
--
Markus Jantti
Abo Akademi University
markus.jantti at iki.fi
http://www.iki.fi/~mjantti