###################################################
### 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: input eval=FALSE
###################################################
## library(nlme)
## url <- "http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08/"
## fev<-groupedData(log.fev1~age|subject, read.table(paste(url,'download/fev1.txt',sep=''),header=T))
## 
## groups <- unique(fev$subject)
## xl<-range(fev$age)
## yl<-range(fev$log.fev1)
## for(i in 1:30)
## {
##    samp <- groups[(1+(i-1)*10):(i*10)]
##    d<-subset(fev,subject %in% samp)
##    print(xyplot(log.fev1~age,data=d,groups=subject,type='b', xlim=xl, ylim=yl))
##    Sys.sleep(.5)
## }
## 
## xl<-range(fev$height)
## for(i in 1:30)
## {
##    samp <- groups[(1+(i-1)*10):(i*10)]
##    d<-subset(fev,subject %in% samp)
##    print(xyplot(log.fev1~height,data=d, groups=subject,type='b',xlim=xl, ylim=yl))
##    Sys.sleep(.5)
## }


###################################################
### chunk number 3: input2
###################################################
library(nlme)
url <- "http://www.statistik.lmu.de/institut/lehrstuhl/semwiso/longitudinal_ss08/"
#fev<-groupedData(log.fev1~age|subject, read.table(paste(url,'download/fev1.txt'),header=T))
fev<-groupedData(log.fev1~age|subject, read.table('W:/lehre/longitudinal_sommer08/uebung/fev1.txt',header=T))


###################################################
### chunk number 4: a
###################################################
m00<-lme(log.fev1~ age + height,  random=~age|subject, data=fev, method='ML')
m10<-lme(log.fev1~ poly(age,2) + height,  random=~poly(age,2)|subject, data=fev, method='ML')
m20<-lme(log.fev1~ log(age) + height,  random=~log(age)|subject, data=fev, method='ML')
m01<-lme(log.fev1~ age + log(height),  random=~age|subject, data=fev, method='ML')
m11<-lme(log.fev1~ poly(age,2) + log(height),  random=~poly(age,2)|subject, data=fev, method='ML')
m21<-lme(log.fev1~ log(age) + log(height),  random=~log(age)|subject, data=fev, method='ML')
m02<-lme(log.fev1~ age + poly(height,2),  random=~age|subject, data=fev, method='ML')
m12<-lme(log.fev1~ poly(age,2) + poly(height,2),  random=~poly(age,2)|subject, data=fev, method='ML')
m22<-lme(log.fev1~ log(age) + poly(height,2),  random=~log(age)|subject, data=fev, method='ML')

anova(m22,m12,m02,m21,m11,m01,m20,m10,m00,test=F)


###################################################
### chunk number 5: a.2
###################################################
summary(m12)


###################################################
### chunk number 6: a.inter
###################################################
m12.inter <- update(m12, .~. + age:height)
anova(m12,m12.inter)


###################################################
### chunk number 7: a1
###################################################
print(plot(m12, subject ~ resid(.,type='n'), abline=0, cex=.5))


###################################################
### chunk number 8: a2
###################################################
print(plot(m12, resid(.)~ age | subject, abline=0, cex=.5, strip=F))


###################################################
### chunk number 9: a3
###################################################
print(plot(m12, resid(., type='p')~  age,  cex=.7))


###################################################
### chunk number 10: a4
###################################################
print(qqnorm(m12, ~resid(. ,type='n'), abline=c(0,1),  cex=.7))


###################################################
### chunk number 11: a42
###################################################
plot(density(resid(m12 ,type='n')),main='Dichteschätzung der normalisierten Residuen \n und Dichte der N(0,1) ')
curve(dnorm(x),from=-6,to=3,add=T,col=2)


###################################################
### chunk number 12: a5
###################################################
(Vg.12 <- Variogram(m12, form= ~age|subject) )
print(plot(Vg.12))


###################################################
### chunk number 13: b
###################################################
print(qqnorm(m12, ~ranef(.),  cex=.7))


###################################################
### chunk number 14: c
###################################################
try(m12.Exp <- update(m12, corr=corExp(form=~age|subject)))
try(m12.Gauss <- update(m12, corr=corGaus(form=~age|subject)))
m12.diag<-update(m12,random=pdDiag(~poly(age,2)))
m12.Exp <- update(m12.diag, corr=corExp(form=~age|subject))
m12.Gauss <- update(m12.diag, corr=corGaus(form=~age|subject))
anova(m12, m12.Exp, m12.Gauss,test=F)


###################################################
### chunk number 15: c1
###################################################
getVarCov(m12)
getVarCov(m12.diag)
getVarCov(m12.Exp)


###################################################
### chunk number 16: c2
###################################################
Vg.Exp <- plot(Variogram(m12.Exp, form= ~age|subject, resType='p'),ylim=c(.4,1.2))
Vg.12 <- plot(Variogram(m12, form= ~age|subject, resType='p'),ylim=c(.4,1.2))
print(Vg.Exp, split=c(1,1,1,2), more=TRUE)
print(Vg.12, split=c(1,2,1,2))


###################################################
### chunk number 17: c3
###################################################
res.serial<-plot(m12.Exp,subject ~ resid(. ,type='n'), xlim=c(-6,4),cex=.5,abline=0,
                 main='Normalisierte Residuen, mit Exponential-Korrelation')
res.noserial<-plot(m12,subject ~ resid(. ,type='n'), ,xlim=c(-6,4),cex=.5, abline=0,
                 main='ohne Exponential-Korrelation')
print(res.serial, split=c(1,1,1,2), more=TRUE)
print(res.noserial, split=c(1,2,1,2))


###################################################
### chunk number 18: 2a
###################################################
mahalanobis.residcheck<-function(m,p=.01)
{
 #get grouping factor
 grps <- getGroups(m)
 # get no. of subjects' observations
 ng <- table(grps)
    #get group levels
 grplevels  <- levels(grps)

 #raw residuals y - X beta
 r <- resid(m, level=0, type='response')
 # divide into groups
 r <- split(r,grps)

 for(i in grplevels)
 {
            #get marginal covariance Z_i D Z_i' + H_i + sigma^2 I_{n_i}
            #for ng=1 getVarCov doesn't work with correlation structures
            #so we don't calculate anything for single individuals here.
            #(constructing V by hand is complicated)
            if(ng[i]>1){
            V <- getVarCov(m,individual=i,type='marginal')
            # multiply group residuals with inverse matrix root of V
            r[[i]] <- backsolve(chol(V[[1]]),diag(ng[i]))%*%r[[i]]
            } else  r[[i]] <-NA
 }

 #calculate mahalanobis distance (= r'r) for each subject ...
 d <- sapply(r, crossprod)

 #... and calculate P(d_i > ChiSq(ng_i) ) for all i
 pvals <-  1 - pchisq(q=d, df=ng)

 res<-data.frame(N=as.numeric(ng),P=pvals, out = pvals < p)
 #remove NAs
 return(res[complete.cases(res),])
}


###################################################
### chunk number 19: 2b
###################################################
mah <- mahalanobis.residcheck(m12.Exp,.01)
(out<-rownames(mah)[mah$out])


###################################################
### chunk number 20: 2c
###################################################
m12.Exp.noout <- update(m12.Exp, data=subset(fev, !(subject %in% out)))

summary(m12.Exp.noout)
summary(m12.Exp)

r.m12.noout <-plot(m12.Exp.noout, subject~resid(.,type='n'),cex=.7, main = 'Ausreißer Subjekte entfernt',xlim=c(-7,5))
r.m12 <-plot(m12.Exp, subject~resid(.,type='n'),cex=.7, main = 'mit Ausreißern',xlim=c(-7,5))
print(r.m12.noout, split=c(1,1,1,2), more=TRUE)
print(r.m12, split=c(1,2,1,2))