################################################################################
# Exercise 2
# Bayesian LMM
################################################################################

library(nlme)
library(MCMCglmm)
?MCMCglmm

# Erzeuge Daten für 100m-Lauf
set.seed(123)
n1 = 10
n2 = n1
n = n1 + n2
M = 3
b = rep(rnorm(n1 + n2, 0, 0.5), each = M)
hist(b)
sex = c(rep(0,n1), rep(1,n2))
beta = c(rep(10.6, n1 * M), rep(11.9, n2 * M))
eps = rnorm(n, 0, 0.2)
times = beta + b + eps
athlete = rep(1:n, each = M)

plot(times ~ athlete, col = athlete, xlab = "LaeuferIn", ylab = "Zeit [sek]")
abline(v=10.5)

sprint <- data.frame(athlete=factor(athlete), sex=factor(rep(sex, each=M)), times)


#### (a) 
# Estimate random intercept model by REML
sprintREML <- lme(times ~  sex, random=~1|athlete, data=sprint, method="REML")
summary(sprintREML)


#### (b) 
# Estimate random intercept model using MCMC
sprintMCMC <- MCMCglmm(times ~  sex, random = ~athlete, data=sprint,
                 nitt=13000, thin=10, burnin=3000, pr=TRUE)
summary(sprintMCMC)
# "pr=TRUE": save posterior distribution of random effects

# # Check der autocorrelation of MCMC-Samples
# # increase "thin" and "burnin" if necessary
# autocorr(sprintMCMC$Sol, lag=1)
# # small values of autocorrelation


#### (c) fixed effects
# posterior mode of fixed effects
posterior.mode(sprintMCMC$Sol[,c("(Intercept)", "sex1")])

# posterior expectation of fixed effects
apply(sprintMCMC$Sol[,c("(Intercept)", "sex1")], MARGIN=2, FUN=mean)

# posterior median of fixed effects
apply(sprintMCMC$Sol[,c("(Intercept)", "sex1")], MARGIN=2, FUN=median)

# fixed effects of REML-estimation
fixef(sprintREML)

# true values
10.6 # Intercept
11.9 - 10.6 # sex1


#### (d) variance components

# variance components
# VCV  Posterior Distribution of (co)variance matrices
apply(sprintMCMC$VCV, MARGIN=2, FUN=mean) 


# variance components in REML model
getVarCov(sprintREML)
sprintREML$sigma^2

# true values
0.5^2   # var(b)
0.2^2   # sigma^2



#### (e) random effects

# posterior mode of random effects
posterior.mode(sprintMCMC$Sol[,-c(1:2)])

# posterior expectation of random effects
apply(sprintMCMC$Sol[,-c(1:2)], MARGIN=2, FUN=mean)

# posterior median of random effects
apply(sprintMCMC$Sol[,-c(1:2)], MARGIN=2, FUN=median)


# compare random effects between MCMC and REML
summary(apply(sprintMCMC$Sol, MARGIN=2, FUN=mean)[-c(1:2)]-random.effects(sprintREML))

# compare both with the true random effects
summary(apply(sprintMCMC$Sol, MARGIN=2, FUN=mean)[-c(1:2)]-b[seq(1, n*M, by=M)])
summary(random.effects(sprintREML)-b[seq(1, n*M, by=M)])



#### (f) 
# HPD-Intervalle
HPDinterval(sprintMCMC$Sol, prob = 0.95)

# # Quantiles
# apply(sprintMCMC$Sol, MARGIN=2, FUN=quantile, probs = c(0.025,0.975) )

# # posterior variance of effects
# apply(sprintMCMC$Sol, MARGIN=2, FUN=var)
# apply(sprintMCMC$VCV, MARGIN=2, FUN=var)


#### (g) 
# Markov-Chain and kernel density estimation
plot(sprintMCMC$Sol)     # fixed and random effects
plot(sprintMCMC$VCV)     # "units": estimation of error variance