###################################################
### 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(nlme)
m0 <- lme(distance~Sex*I(age-11), random=~1|Subject, data=Orthodont)
m1 <- update(m0, random=~I(age-11)|Subject)
(lr01 <- anova(m0,m1))


###################################################
### chunk number 3: 1b
###################################################
pchisqmix <- function(q, df = c(1, 2)) {
 sum(0.5 * pchisq(rep(q, 2), df))
}
1 - pchisqmix(lr01$L.Ratio[2], df = c(1, 2))
lr01$"p-value"[2]


###################################################
### chunk number 4: 1c
###################################################
n.sim <- 1000
lr01.sim <- simulate.lme(m0,m2=m1, nsim = n.sim)
print(plot(lr01.sim, df = c(1, 2)))
lrt01 <- 2 * (lr01.sim$alt$REML[, "logLik"] - lr01.sim$null$REML[, "logLik"])
(p01.sim <- mean(lrt01 > lr01$L.Ratio[2]))
1 - pchisqmix(lr01$L.Ratio[2], df = c(1, 2))
lr01$"p-value"[2]


###################################################
### chunk number 5: 1d1 eval=FALSE
###################################################
## m0 <- lme(distance~Sex*I(age-11), random=~1|Subject, data=Orthodont)
## m1 <- update(m0, random=~I(age-11)|Subject)
## m2 <- update(m0, random=list(Subject = pdDiag(~I(age-11))))
## summary(m0)
## summary(m1)
## summary(m2)


###################################################
### chunk number 6: 1d2 eval=FALSE
###################################################
## (lr21 <- anova(m2,m1))
## lr21.sim <- simulate.lme(m2,m2=m1, nsim = n.sim)
## print(plot(lr21.sim, df = c(1,2)))
## lrt21 <- 2 * (lr21.sim$alt$REML[, "logLik"] - lr21.sim$null$REML[, "logLik"])
## (p21.sim <- mean(lrt21 > lr21$L.Ratio[2]))
## 1 - pchisq(lr21$L.Ratio[2], df = 1)


###################################################
### chunk number 7: 1d3 eval=FALSE
###################################################
## (lr02 <- anova(m0,m2))
## lr02.sim <- simulate.lme(m0,m2=m2, nsim = n.sim)
## print(plot(lr02.sim, df = c(0,1)))
## lrt02 <- 2 * (lr02.sim$alt$REML[, "logLik"] - lr02.sim$null$REML[, "logLik"])
## (p02.sim <- mean(lrt02 > lr02$L.Ratio[2]))
## 1 - 0.5*pchisq(lr02$L.Ratio[2], df = 1) - 0.5


###################################################
### chunk number 8: 1dsol1
###################################################
m0 <- lme(distance~Sex*I(age-11), random=~1|Subject, data=Orthodont)
m1 <- update(m0, random=~I(age-11)|Subject)
m2 <- update(m0, random=list(Subject = pdDiag(~I(age-11))))
summary(m0)
summary(m1)
summary(m2)


###################################################
### chunk number 9: 1dsol2
###################################################
(lr21 <- anova(m2,m1))
lr21.sim <- simulate.lme(m2,m2=m1, nsim = n.sim)
print(plot(lr21.sim, df = c(1,2)))
lrt21 <- 2 * (lr21.sim$alt$REML[, "logLik"] - lr21.sim$null$REML[, "logLik"])
(p21.sim <- mean(lrt21 > lr21$L.Ratio[2]))


###################################################
### chunk number 10: 1dsol3
###################################################
(lr02 <- anova(m0,m2))
lr02.sim <- simulate.lme(m0,m2=m2, nsim = n.sim)
print(plot(lr02.sim, df = c(0,1)))
lrt02 <- 2 * (lr02.sim$alt$REML[, "logLik"] - lr02.sim$null$REML[, "logLik"])
(p02.sim <- mean(lrt02 > lr02$L.Ratio[2]))
1 - 0.5*pchisq(lr02$L.Ratio[2], df =1) - 0.5


###################################################
### chunk number 11: 2ab
###################################################
######################################################################
# Author: Michael Höhle <http://www.stat.uni-muenchen.de/~hoehle>
# Date: 7 November 2004
# ######################################################################

#Time series of German shares found in the Fahrmeir et. al
#exercise book.
wertpap <- c(7.51 , 7.42 , 6.76 , 5.89 , 5.95 , 5.35 , 5.51 , 6.13 ,
6.45 , 6.51 , 6.92, 6.95 , 6.77 , 6.86 , 6.95 , 6.66 , 6.26 , 6.18 ,
6.07 , 6.52 , 6.52 , 6.71)

#Make the plot
plot.ts(wertpap)


###################################################
### chunk number 12: 2ab2
###################################################
#Moving averages / gleitende durchschnitte
ma3 <- filter(wertpap,rep(1/3,3))
ma11<- filter(wertpap,rep(1/11,11))

#Plot the two results and the original series in the same plot.
plot.ts(cbind(wertpap,ma3,ma11),plot.type="single",col=c(1,2,3), ylab='wertpap')


###################################################
### chunk number 13: 2c
###################################################
######################################################################
#A median filter on -q,-q+1,..-1,0,1,...,q+1,q.
#Here it is not possible to use the filter function
#
#Params:
# x - the time series
# q - width of the window
#
#Returns:
# median filtered version of x. Note that the first q and the last q values
# are NA.
######################################################################
filter.med <- function(x,q=1) {
  #Start with q NA's
  res <- rep(NA,q)
  #Loop over all computable elements in the ts
  for (i in (q+1):(length(x)-q)) {
    res <- append(res,median(x[(i-q):(i+q)]))
  }
  #Add NA's
  res <- append(res,rep(NA,q))
  return(res)
}

#Apply median filters and plot the result(s)
me3  <- filter.med(wertpap,1)
me11 <- filter.med(wertpap,5)
plot.ts(cbind(wertpap,me3,me11),plot.type="single",col=c(1,2,3),ylab='wertpap')


###################################################
### chunk number 14: 3b
###################################################
######################################################################
# Fit a polynomial of degree 2 to the 5 values y at time -2,-1,0,1,2
#
# Params:
#  y - The data points to fit.
#
# Returns:
#  \hat{y}_0 = \hat{\beta}_0
######################################################################
fit.poly2 <- function(y) {
  n <- length(y)
  if(n!=5) stop('wrong length for y!')
  t <- c(-2,-1,0,1,2)

  m <- lm(y ~ 1 + t + I(t^2))
  yhat <- coef(m)[1]
  return(yhat)
}

######################################################################
#Apply the local polynomial filter manually. Basic version.
#
#Params:
# y - time series
######################################################################
filter.poly <- function(y,q=2) {
  #Initialize with NAs.
  res <- rep(NA,length(y))

  #Loop over all
  for (i in (q+1):(length(y)-q)) {
    yhat <- fit.poly2(y[(i-q):(i+q)])
    res[i] <- yhat
  }
  #End with NAs
  return(res)
}


######################################################################
# General polynomial fit function.
#
# Params:
#  y - data points to interpolate using a polynomial of the values
#      t-q,t-q+1,...,t+1,t,t+1,...,t+q-1,t+q around time t.
#  p - order of the poly to fit.
#
# Returns:
#  The prediction for \hat{y}_0 = \beta_0
######################################################################

fit.poly <- function(y,p,pred="t") {
  #Determine window width - note we don't catch errors when q is not
  #an integer!
  q <- (length(y)-1)/2
  t <- seq(-q,q)
  model <- "y~1"
  #Build model formula
  for (i in 1:p) {
    model <- paste(model," + I(t^",i,")",sep="")
  }
  #Fit polynomial using KQ.
  m <- lm(as.formula(model))
  #Predict the right thing, i.e. either y_t, y_{-q}:y_{-1} or
  #y_{+1}:y_{q}
  switch(pred,
         "t" =    {yhat = predict(m)[q+1]},
         "left" = {yhat = predict(m)[1:(q)]},
         "right"= {yhat = predict(m)[(q+2):(2*q+1)]})

  return(yhat)
}

######################################################################
#Apply the local polynomial filter manually including an edge
#correction.
#
#Params:
# y - time series
######################################################################
filter.poly.edge <- function(y,q=2,p=2) {
  #Start with border correction for the first values
  res <- fit.poly(y[1:(2*q+1)],p=p,pred="left")
  #Loop over all
  for (i in (q+1):(length(y)-q)) {
    yhat <- fit.poly(y[(i-q):(i+q)],p=2)
    res <- append(res,yhat)
  }
  #End with edge correction
  n <- length(y)
  res <- append(res,fit.poly(y[(n-2*q):n],p=p,pred="right"))
  #Done
  return(res)
}


#Polynomial filter
fP <- filter.poly(wertpap,q=2)
plot.ts(cbind(wertpap,fP),col=1:2,plot.type="single", ylab='wertpap')

#Filter solution equivalent.
D  <- 1/35*c(-3,12,17,12,-3)
fD <- filter(wertpap,D)

#Difference (small!)
fP-fD


###################################################
### chunk number 15: 3b2
###################################################
##Edge corrected version
fPe <- filter.poly.edge(wertpap,q=2,p=2)
plot.ts(cbind(wertpap,fPe),col=1:2,plot.type="single", ylab='wertpap')