#===================================
# Quasi-likelihood for binomial data
# Example: Rat data (section 8.2.4)
#==================================

# Read rats data from the web
data="http://www.stat.ufl.edu/~aa/glm/data/Rats.dat"
Rats=read.table(data,header=T)
head(Rats)

# Consider two groups: placebo or treated
Rats$placebo=ifelse(Rats$group==1,1,0)

# Fit binomial model
fit.binom=glm(s/n~placebo+h,weights=n,family=binomial,data=Rats)
summary(fit.binom)

# Compute overdispersion parameter phi:
n=dim(Rats)[1]
p=length(fit.binom$coefficients)
X2=sum(residuals(fit.binom,type="pearson")^2)  
phi=X2/(n-p)
phi

# The computations may be done directly using the quasi-family:
fit.quasibinom=glm(s/n~placebo+h,weights=n,
                   family=quasi(link="logit",variance="mu(1-mu)"),data=Rats)
summary(fit.quasibinom)


#========================================
# Sandwich estimator of covariance matrix
# (cf section 8.3.4)
#========================================

library(gee)
obs=1:dim(Rats)[1]

# Fit model using estimating equations:
fit.gee=gee(cbind(s,n-s)~placebo+h,id=obs, family=binomial,scale.fix=T,data=Rats)
summary(fit.gee)