# Example 13.2 in Diggle et al., 2002: analysis of Cow milk protein dataset

week <- c(1:19)
m <- matrix(scan("e:/vhm/vhm882/data/cowmilkp.dat",na.strings="0.00"),ncol=19,byrow=T)
#m <- matrix(scan("d:/data.avc/teaching/vhm882/13l/cowmilkp.dat",na.strings="0.00"),ncol=19,byrow=T)
group <- c(rep(1,25),rep(2,27),rep(3,27))
cow <- c(1:79)

# code for graphics, using nlme and lattice libraries
library(nlme)
cows <- balancedGrouped( y ~ week|cow, matrix(m,nrow=79,ncol=19,dimnames=list(cow,week)),
                           labels=list(y="Protein content in weekly milk samples"))
cows$Group <- as.factor(rep(group,rep(19,79)))

# profile plots
plot(cows, outer=~Group, aspect=2)
plot(cows[cows$Group==1,,], outer=~Group, aspect=1)

# mean plot
cows.means <- aggregate(cows$y,list(cows$week,cows$Group),mean,na.rm=TRUE)
names(cows.means)[1]<-"Week"
cows.means$Week <- as.numeric(levels(cows.means$Week))[cows.means$Week]
names(cows.means)[2]<-"Group"
library(lattice)
xyplot(x~Week|Group,cows.means,ylab="Mean Protein content in weekly milk samples")

# definition of binary outcomes
drop <- matrix(0,ncol=4,nrow=79)
drop[,1] <- as.numeric(is.na(m[,15]) & is.na(m[,16]) & is.na(m[,17]) & is.na(m[,19]))
drop[,2] <- 10*as.numeric(drop[,1]) + as.numeric(is.na(m[,16]) & is.na(m[,17]) & is.na(m[,19]))
drop[,3] <- 10*as.numeric(drop[,1] | drop[,2]) + as.numeric(is.na(m[,17]) & is.na(m[,19]))
drop[,4] <- 10*as.numeric(drop[,1] | drop[,2] | drop[,3]) + as.numeric(is.na(m[,19]))
dropweek <- c(15,16,17,19)
dcows <- balancedGrouped( z ~ dropweek|cow, matrix(drop,nrow=79,ncol=4,dimnames=list(cow,dropweek)))
is.na (dcows$z) <- dcows$z>10
dcows$Group <- as.factor(rep(group,rep(4,79)))
dcows$prev <- c(t(matrix(c(m[,14],m[,15],m[,16],m[,18]),ncol=4,nrow=79,byrow=F)))
dcows$dropweek <- as.factor(dcows$dropweek)

dcows.glm1 <- glm(z ~ Group:dropweek+Group:dropweek:prev, data=dcows, family=binomial)
summary(dcows.glm1)
dcows.glm1a <- glm(z ~ Group:dropweek+Group:prev+dropweek:prev, data=dcows, family=binomial)
summary(dcows.glm1a)
dcows.glm2 <- glm(z ~ Group:dropweek+Group:prev, data=dcows, family=binomial)
summary(dcows.glm2)
dcows.glm3 <- glm(z ~ Group:dropweek+dropweek:prev, data=dcows, family=binomial)
summary(dcows.glm3)
dcows.glm4 <- glm(z ~ Group:dropweek+prev, data=dcows, family=binomial)
summary(dcows.glm4)
dcows.glm5 <- glm(z ~ Group*dropweek, data=dcows, family=binomial)
summary(dcows.glm5)
dcows.glm6 <- glm(z ~ Group+dropweek+prev, data=dcows, family=binomial)
summary(dcows.glm6)
dcows.glm7 <- glm(z ~ Group+prev, data=dcows, family=binomial)
summary(dcows.glm7)
dcows.glm7a <- glm(z ~ dropweek+prev, data=dcows, family=binomial)
summary(dcows.glm7a)
dcows.glm8 <- glm(z ~ prev, data=dcows, family=binomial)
summary(dcows.glm8)

library(gee)
dcows.gee8 <- gee(z ~ prev, data=dcows, family=binomial, id=cow)
summary(dcows.gee8)
# also evidence of clustering => not MCAR randomness