GEE

### GEE ?
library(nlme)
lme(travel~1, random=~1|Rail, method="ML", data=Rail)

y<-Rail$travel
id<-Rail$Rail
time<-rep(1:3,6)

X<-rep(1,length(y))
mu<-solve(t(X)%*%X)%*%t(X)%*%y

Y<-cbind(y[time==1],y[time==2],y[time==3])
Res<-Y-matrix(rep(mu,18),6)
R<-cor(Res)
phi<-1/var(y)

V<-matrix(rep(1,9),3)
for(i in 1:3) for(j in 1:3) if(i!=j) V[i,j]=mean(c(R[1,2],R[1,3],R[2,3]))

solve(t(X)%*%solve(kronecker(V,diag(6)))%*%X)%*%t(X)%*%solve(kronecker(V,diag(6)))%*%y

alpha<-mean((Res[,1]*Res[,2]+Res[,2]*Res[,3])/2)/phi