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

Logistic Regression - estimation

### GLM

## logistic regression
respire<-read.csv("F:/publish/data analysis using R/data/respire.csv")
respire2<-read.csv("F:/publish/data analysis using R/data/respire2.csv")

out<-glm(outcome ~ treat, weights=count, family=binomial, data=respire)

logist.lm<-function(beta){
p0<-exp(beta[1])/(1+exp(beta[1]))
p1<-exp(beta[1]+beta[2])/(1+exp(beta[1]+beta[2]))
ll<-16*log(p0)+48*log(1-p0)+40*log(p1)+20*log(1-p1)
return(-ll)
}
optim(c(-1,2),logist.lm)

library(rootSolve)

model<-function(beta){
x<-(respire$treat=="test")
y<-respire$outcome
mu<-exp(beta[1]+beta[2]*x)/(1+exp(beta[1]+beta[2]*x))
W<-respire$count
F1<-sum(W*(y-mu))
F2<-sum(W*(y-mu)*x)
c(F1=F1,F2=F2)
}
multiroot(f=model,start=c(-1,2))$root