# [R] help with gls

Markus Jäntti markus.jantti at iki.fi
Sat Jan 6 13:16:08 CET 2007

```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},
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