###################################################
### chunk number 1: prelim
###################################################
rm(list=ls())
try(detach(package:nlme))
try(detach(package:lme4))
setwd('W:/lehre/longitudinal_sommer08/uebung')
options(width = 80, prompt = "R>")


###################################################
### chunk number 2: 1data eval=FALSE
###################################################
## url <- "http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08/"
## schizo <- read.table(paste(url,"download/fev1.txt", sep = ""), header = T)


###################################################
### chunk number 3: 1data2
###################################################
schizo <- read.table('schizo.txt')


###################################################
### chunk number 4: 1a
###################################################
(mean.early <- with(subset(schizo, onset=='< 20 yrs'),
                    tapply(unclass(disorder)-1,month,mean, na.rm=T)
                    ))
(n.early <- with(subset(schizo, onset=='< 20 yrs' & !is.na(disorder)),
                    tapply(disorder,month,length)
                    ))

(mean.late <- with(subset(schizo, onset=='> 20 yrs'),
                    tapply(unclass(disorder)-1,month,mean, na.rm=T)
                    ))
(n.late <-with(subset(schizo, onset=='> 20 yrs' & !is.na(disorder)),
                    tapply(disorder,month,length)
                    ))

plot(c(0,2,6,8,10), mean.early, type='b', ylim=c(0,1), cex=n.early/10,
     main='Relative Häufigkeiten')
lines(c(0,2,6,8,10), mean.late, type='b', cex=n.late/10, col=2)
legend('topright', col=1:2, lwd=1, legend= c('early onset','late onset'))


###################################################
### chunk number 5: 1b1
###################################################
library(lme4)
library(MASS)
library(glmmML)

glm.fixed <- glm(disorder ~ onset*month, family=binomial, data=schizo)
glmm.laplace<-lmer(disorder ~ onset*month + (1|subject), family=binomial, data=schizo, scale=F)
glmm.PQL<-lmer(disorder ~ onset*month + (1|subject), family=binomial, data=schizo, method='PQL')
glmm.PQL2 <- glmmPQL(disorder ~ onset*month, random = ~1|subject, family=binomial, data=schizo)

#influence of number of quadrature points?
glmm.GQ5 <- glmmML(I(unclass(disorder)-1) ~ onset*month, cluster=subject, family=binomial, data=schizo, method='ghq', n.points=5)
glmm.GQ30 <- glmmML(disorder ~ onset*month, cluster=subject, family=binomial, data=schizo, method='ghq', n.points=30)

coef.list <- list(fixed=summary(glm.fixed)$coefficients[,1:2],
                  Laplace=summary(glmm.laplace)@coefs[,1:2],
                  PQL.lme4=summary(glmm.PQL)@coefs[,1:2],
                  PQL=summary(glmm.PQL2)$tTable[,1:2],
                  GQ=cbind(glmm.GQ30$coefficients,glmm.GQ30$coef.sd))


plot.coef <- function(coef.list, se.factor=1.96)
#coef.list: list of matrices with coefficient values and se's for a number of models
#se.factor: show estimate +/- se.factor * se
{
  m <- length(coef.list)
  p <- sapply(coef.list, NROW)
  if(min(p) == max(p)) p <- p[1]
  else stop('different no. of parameters in models')
  par(mfrow=c(ceiling(p/2),2))
  nam <- rownames(coef.list[[1]])

  tab<-matrix(NA,nrow=m,ncol=3)
  for(i in 1:p)
  {
    for(j in 1:m)
    {
      tab[j,2] <- coef.list[[j]][i,1]
      se <- coef.list[[j]][i,2]
      tab[j,c(1,3)] <- tab[j,2] + c(-1,1)* se.factor * se
    }
    plot(tab[,2],1:m,main=nam[i],col='white',xlim=range(tab),yaxt='n',xlab='',ylab='')
    abline(h=1:m,col='grey')
    abline(v=0,col='darkgrey')
    points(tab[,2],1:m,pch=19)
    segments(tab[,1], 1:m, tab[,3], 1:m)
    points(tab[,1],1:m,pch='(')
    points(tab[,3],1:m,pch=')')
    axis(2,at=1:m,labels=names(coef.list),las=2)
  }
}

plot.coef(coef.list)


###################################################
### chunk number 6: 1b2
###################################################
deviance(glmm.PQL)
deviance(glmm.laplace)
deviance(glmm.GQ30)