#==============================
# GLMM and marginal models
# Example: Attitude to abortion
#==============================

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

# Creat dummy variables for situations 1 and 2:
abortion$z1=ifelse(abortion$situation==1,1,0)
abortion$z2=ifelse(abortion$situation==2,1,0)

# For models with a random intercept, we may use the glmmML-library 
library(glmmML)

# Fit model with random effect for person:
fit.glmmML=glmmML(response~gender+z1+z2, cluster=case, family=binomial, data=abortion,method="ghq", n.points=70, start.sigma=9) 
summary(fit.glmmML)


# Alternatively we may use the lme4-library
# (which allows for more than one random effect)
library(lme4)

# Fit the same model as above
fit.lme4=glmer(response~gender+z1+z2+(1|case), family=binomial,nAGQ=70,data=abortion) 
summary(fit.lme4)


# Then consider marginal models

# Use the gee-package
library(gee)

# Fit marginal model assuming "exchangeable correlation":
fit.exc=gee(response~gender+z1+z2,id=case, family=binomial, corstr="exchangeable", data=abortion)
summary(fit.exc)