[R] help with gls

[Markus Jäntti]

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