# Based on Code by Fabian Scheipl
###############################################################################

################################################################################
# Exercise 1
# LMM
################################################################################

library(MASS)
library(lattice)
library(nlme)

ex <- read.table("amsterdam.txt", header=TRUE)

head(ex) 

ex$smoker <- factor(ex$smoker, labels=c("non-smoker","smoker"))
ex$gender <- factor(ex$gender, labels=c("female", "male"))

###### (1a) 
## Deskription

## Spaghettiplot der TC-Messungen fuer die vier moeglichen Kombinationen von smoker und gender
xyplot(tc ~ t | smoker*gender, groups=id, data=ex, type="b")

## Boxplots der TC-Messungen ueber die Zeit fuer die vier moeglichen Kombinationen von smoker und gender
bwplot(tc ~ as.factor(t)| smoker*gender, groups=id, data=ex)

## Mittelwerte und Standardabweichungen der TC-Messungen fuer die beiden Geschlechter
t <- unique(ex$t)
means <- with(ex,tapply(tc,list(t,gender),mean))
sds <- with(ex,tapply(tc,list(t,gender),sd))
res <- cbind(means[,1], sds[,1], means[,2], sds[,2])
colnames(res) <- c("mean(F)","sd(F)","mean(M)","sd(M)")
res


###### (1b)
## Modellierung

## LMM mit zeitlichem Trend sowie Random-Intercept und Random-Slope
m.slope.lme <- lme(tc ~ 1 + t, random= ~ 1+t|id, data=ex)
summary(m.slope.lme)

## wichtige Parameter
fixef(m.slope.lme)
m.slope.lme$sigma^2
getVarCov(m.slope.lme)

# plot der vorhergesagten Werte
plot(m.slope.lme, form=fitted(.) ~ t, type="b")



###### (1c)

## marginale Kovarianzmatrix der Beobachtungen eines Individuums
t <- unique(ex$t)
Zi <- cbind(1, t)
D <- getVarCov(m.slope.lme)
R <- diag(m.slope.lme$sigma^2, 6)
(Vi <- Zi %*% D %*% t(Zi) + R)

getVarCov(m.slope.lme, individuals="1", type="marginal")

## Kovarianzmatrix fuer den gesamten Datensatz ist blockdiagonale Matrix,
## mit Kovarianzmatrizen der Individuen auf der Diagonale
## erzeuge blockdiagonale Matrix fuer n=2 Individuen
require(Matrix)
temp <- list(Vi = Vi)
listVi <- rep(temp, 2) # Liste mit 2 Kovarianzmatrizen Vi 
bdiag(listVi) # Kovarianzmatrix fuer 2 Subjekte

# Korrelationsmatrix fuer ein Subjekt
cov2cor(Vi)

# ## Extra: Residuencheck
# plot(m.slope.lme)
# xyplot(resid(m.slope.lme) ~ t|smoker*gender, groups=id, data=ex, type="l")
# 
# ## Extra: Modell mit zusaetzlichen festen linearen Effekten
# m.full.lme <- update(m.slope.lme,.~. + smoker + gender*t + fitness + bfatness)
# summary(m.full.lme)
# xyplot(resid(m.full.lme) ~ t|smoker*gender, groups=id, data=ex, type="l")



###### (1d) Zusatzaufgabe 
## positive und negative KOrrelation von random intercept und random slope im LMM
cor.pos <- matrix(c(1,0.9,0.9,1),2,2)
cor.neg <- matrix(c(1,-0.9,-0.9,1),2,2)
library(mvtnorm)
# Erzeuge Plot fuer Geraden aus Intercepts und Slopes
# zufaellige Effekte: (b_i1, b_i1)' ~ N( (0, 0)',  corM), 
# feste Effekte: beta_0 = 5, beta_1 = 2
do.plot <- function(corM, main) {
  # Simuliere N( (5, 2)',  corM) 
  res <- mvrnorm(n=10, mu=c(5, 2), Sigma=corM)
  # Plotte Effekte 
  plot(NA, xlim=c(0,5), ylim=c(3,10), type="n", xlab="t", ylab="", main=main)
  for (i in 1:nrow(res)) abline(coef=res[i,])
  abline(coef=c(5,2), lty=1, lwd=4, col='red') # Zeichne feste Effekte ein
}
x11(w=12, h=6)
par(mfcol=c(1,2))
set.seed(123)
do.plot(cor.pos, "Positive Korrelation")
do.plot(cor.neg, "Negative Korrelation")


###### (1e) smooth effect in time for males and females
require(mgcv)
m.full.amm <- gamm(tc ~ 1 + s(t, by=gender, k=5) + bfatness + smoker, 
                   random= list(id=~1+t), data=ex)
summary(m.full.amm$gam)
summary(m.full.amm$lme)

plot(m.full.amm$gam, pages=1)

# plot der Residuen
xyplot(resid(m.full.amm$lme) ~ t|smoker*gender, groups=id, data=ex, type="l")

# plot der Prognosen
xyplot(fitted(m.full.amm$lme) ~ t|smoker*gender, groups=id, data=ex, type="l")


# ###### (1f) Zusatzaufgabe
# 
# # Lineares Modell fuer jede Person
# ex <- groupedData( tc ~ t| id, ex)
# m.slope.lm <- lmList(ex)
# # unpenCoefs <- coef(m.slope.lm)
# # penCoefs <- coef(m.slope.lme)
# # plot(compareFits(unpenCoefs, penCoefs))
# 
# unpenCoefs <- data.frame(coef(m.slope.lm))
# penCoefs <- data.frame(coef(m.slope.lme))[rownames(unpenCoefs),]
# plot(unpenCoefs, col=1, pch=19, cex=1, ylab="slope", xlab="intercept")
# points(penCoefs, col=4, pch=19, cex=1)
# points(fixef(m.slope.lme)[1], fixef(m.slope.lme)[2], pch=3, cex=3, col=6, lwd=2)
# arrows(unpenCoefs[,1],unpenCoefs[,2], penCoefs[,1], penCoefs[,2], 
#        length=.1, angle=20, col="grey")
# legend("topright", legend=c("mixed model", "individual models"), col=c(4, 1), pch=19)
# 
# 
# # ## alternativ: lineares Modell mit Effekten fuer jede Person, 
# # ## bfatness und smoker sind zeitvariierend
# # m.fix <- gam(tc ~ -1 + id + id:t + s(t, by=gender, k=5) + bfatness + smoker, data=ex)
# # summary(m.fix)