###################################################
### 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: 1a
###################################################
library(ts)

######################################################################
#Filter version
######################################################################
mysim.ar1 <- function(n,alpha,sigma) {
  burnin <- 1000
  x <- rnorm(n+burnin,sigma)
  x <- filter(x, alpha, method = "recursive")
  return(x[-c(1:burnin)])
}

######################################################################
#Manual version
######################################################################
mysim2.ar1 <- function(n,alpha,sigma) {
  burnin <- 1000
  x <- rnorm(n+burnin,sigma)
  for (i in 2:length(x)) {
    x[i] <- alpha*x[i-1] + x[i]
  }
  return(x[-c(1:burnin)])
}

#Direct calculation of y_n using random noise sequence x, x[1]=x_0
yn <- function(n,x) {
  alpha^n * x[1] + sum(alpha^(n-1:n)*x[2:(n+1)])
}


###################################################
### chunk number 3: 1b
###################################################
######################################################################
# Vergleich der zwei Arten die Reihe zu bekommen
###############################################################################
set.seed(1234)
x1 <- mysim.ar1(n=500,alpha=0.8,sigma=1)
set.seed(1234)
x2 <- mysim2.ar1(n=500,alpha=0.8,sigma=1)
#Are those two equal given a tolerence of 1e-12?
sum(x1-x2>1e-12) == 0
plot(ts(x1))


###################################################
### chunk number 4: 1c
###################################################
my.acf <- function(y,lagmax=100) {
 n <- length(y)
 if(n<(lagmax+1)){
    warning('max. Lag longer than time series.\n')
    cat('max. Lag set to ', lagmax <- floor(n/2), '.\n')
 }
 y<-scale(y,scale=F)
 M <- tcrossprod(y)
 r <- rep(NA,lagmax+1)
 for(i in 1:(lagmax+1)){
   r[i] <- sum(diag(M))
   M<-M[-1,]
 }
 return(cbind(lag=0:lagmax,acf = r/sum(y^2)))
}


###################################################
### chunk number 5: 1d
###################################################
x1 <- mysim.ar1(n=1000,alpha=0.8,sigma=1)
plot(my.acf(x1,20))
curve(0.8^x,from=0,to=20,add=T)


###################################################
### chunk number 6: 2
###################################################
require(nlme)
data(BodyWeight)
BW2<-BodyWeight

weight1 <- BW2$weight[BW2$Time==1]
weight1 <- rep(weight1,each=11)
BW2$gain <- 100*(BW2$weight-weight1)/weight1
BW2 <- BW2[BW2$Time!=1,]
BW2$Time <- BW2$Time-1

attr(BW2,'formula') <- gain ~ Time | Rat
attr(BW2,'units')$y <- '(%)'
attr(BW2,'labels')$y <- 'Body weight gain'

print(plot(BW2,  outer=~Diet))


###################################################
### chunk number 7: 2a
###################################################
m<-lme(gain ~ Time:Diet, random=~Time|Rat,  data=BW2)
summary(m)
print(plot(augPred(m)))


###################################################
### chunk number 8: 2a2
###################################################
print(plot(m,resid(.)~Time|Rat,abline=0))


###################################################
### chunk number 9: 2b
###################################################
(Vm <- Variogram(m, form=~Time, maxDist=42))
print(plot(Vm,ylim=c(0,1.1*max(Vm$variog)),main='Semi-Variogram'))


###################################################
### chunk number 10: 2c1
###################################################
mExp <- update(m, corr=corExp(form=~Time))
summary(mExp)


###################################################
### chunk number 11: 2c2
###################################################
mGauss <- update(m, corr=corGaus(form=~Time))
summary(mGauss)


###################################################
### chunk number 12: 2c3
###################################################
anova(m,mExp)
anova(m,mGauss)
anova(mGauss,mExp)


###################################################
### chunk number 13: 2c4
###################################################
VmE <- Variogram(mExp, form=~Time, maxDist=42)
print(plot(VmE,ylim=c(0,1.1*max(VmE$variog)),main='Semi-Variogram Exponential-Correlation'))


###################################################
### chunk number 14: 2c5
###################################################
VmG <- Variogram(mGauss, form=~Time, maxDist=42)
print(plot(VmG,ylim=c(0,1.1*max(VmG$variog)),main='Semi-Variogram Gauss-Correlation'))