# Exercise 9.7 in Davis, 2002; Orthodont data in Pinheiro & Bates
# Continuation of Examples 3.3, 4.1, 4.3 and 6.1

library(nlme)
data(Orthodont)
library(gee)
library(MASS)

# (a): independence working correlation
Orth.gee.ind <- gee(distance~age*Sex, data=Orthodont, family=gaussian, id=Subject, corstr="independence")
summary(Orth.gee.ind)
# Wald test combined for intercept and slope
bhat <- Orth.gee.ind[9]$coefficients
rvar <- Orth.gee.ind[20]$robust.variance
c.sex <- matrix(c(0,0,1,0,0,0,0,1),nrow=2,ncol=4,byrow=T)
wald.sex <- t(c.sex %*% bhat) %*% ginv(c.sex %*% rvar %*% t(c.sex)) %*% c.sex %*% bhat
wald.sex
pchisq(wald.sex,2,lower.tail=F)

# (b): results from random intercept model
Orth.linmix.cs <- lme(distance~age*Sex,random=~1e|Subject,data=Orthodont,method="REML")
summary(Orth.linmix.cs)
# estimates are almost identical, and SE's are close
# for results from linear growth curve analysis, see Example 4.3

# (c): other working correlation structures
Orth.gee.cs <- update(Orth.gee.ind, corstr="exchangeable")
summary(Orth.gee.cs)
# results identical
Orth.gee.unstr <- update(Orth.gee.ind, corstr="unstructured")
summary(Orth.gee.unstr)
# results similar

# (d): cubic model of age
Orth.gee2.ind <- update(Orth.gee.ind,distance~age*Sex+I(age^2)*Sex+I(age^3)*Sex)
summary(Orth.gee2.ind)
# Wald test for non-linearity
bhat <- Orth.gee2.ind[9]$coefficients
rvar <- Orth.gee2.ind[20]$robust.variance
c.nlin <- matrix(c(rep(0,3),1,rep(0,8),1,rep(0,9),1,rep(0,8),1),nrow=4,ncol=8,byrow=T)
wald.nlin <- t(c.nlin %*% bhat) %*% ginv(c.nlin %*% rvar %*% t(c.nlin)) %*% c.nlin %*% bhat
wald.nlin
pchisq(wald.nlin,4,lower.tail=F)