Re: [R] help with gls

From: Markus Jäntti <markus.jantti_at_iki.fi>
Date: Sat 06 Jan 2007 - 12:16:08 GMT

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@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@iki.fi
http://www.iki.fi/~mjantti

______________________________________________
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 23:19:00 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 - 13:30:30 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.