###########################################################################
############# Loesungen zu Aufgaben aus Tutorium 1 ########################
###########################################################################



# Daten einlesen ----------------------------------------------------------

# Arbeitsverzeichnis auf Ordner mit Daten setzen
#setwd("User/Daten")
# Datensatz einlesen
orthodont <- read.table("https://moodle.lmu.de/pluginfile.php/60324/course/section/15118/orthodont.txt", header = TRUE)




# Aufgabe 1 ---------------------------------------------------------------

### a) 
head(orthodont)
#long-Format, da für jede Messung eines Subjektes eine Zeile


### b)
set.seed(2016) # seed setzen zur Ziehung gleicher Zufallszahlen
orthodont$InBehandlung <- sample(0:1, nrow(orthodont), replace = TRUE)

### c)
#install.packages(nlme)
library(nlme)
?groupedData
orthodont2 <- groupedData(distance ~ age | Subject,
                          outer = ~ Sex, inner = ~ InBehandlung, data = orthodont)
# "inner": indicating covariates that are inner to the grouping factor
# "outer": indicating covariates that are outer to the grouping factor
### d)

# i)
plot(orthodont2)
plot(orthodont)

# ii)
plot(orthodont2, outer = TRUE) # alternativ: 'outer = ~ Sex'

# iii)
plot(orthodont2, inner = TRUE) # alternativ: 'inner = ~ InBehandlung'




# Aufgabe 2 ---------------------------------------------------------------

### a)
# siehe Mitschrift


### b)

# i)
library(nlme)
model_lme <- lme(distance ~ age*Sex, random = ~ 1 + age | Subject,
                 data = orthodont)
summary(model_lme)

# ii)
#install.packages(lme4)
library(lme4)
model_lmer <- lmer(distance ~ age*Sex + (1 + age| Subject),
                   data = orthodont)
summary(model_lmer)


### c)

#Vergleich der Fixed Effects
FE_lme <- fixef(model_lme)
FE_lmer <- fixef(model_lmer)

#Vergleich der Random Effects
REs_lme<- ranef(model_lme) # alternativ: model_lme$coefficients$random
REs_lmer <- ranef(model_lmer)$Subject

# Kovarianzmatrizen
getVarCov(model_lme) # Kovarianzmatrizen
cov2cor(getVarCov(model_lme)) # Umwandlung Korrelationsmatrix
VarCorr(model_lmer) # Varianzen & Korrelationen; alternativ: summary(model_lmer)$varcor

# Residualvarianz
attr(VarCorr(model_lmer), "sc")^2 
model_lme$sigma^2

## Vergleich der Modellergebnisse: keine Unterschiede
# Visualisierung des Vergleichs der Random Intercepts
plot(data.frame("lme" = REs_lme[,1], "lmer" = REs_lmer[,1]))

# i)

# 'type = "marginal"': Marginale Kovarianzmatrix des Response
varcov_lme <- getVarCov(model_lme, individuals = "M01", type = "marginal")
cov2cor(varcov_lme$M01)
# -> keine Diagonalmatrix, da hier die Random Effects mit in die Kovarianzmatrix einfliessen
# -> relativ hohe Korrelationen innerhalb der Subjekte
#    -> GLM mit Unabhaengigkeitsannahme waere hier z.B. unangebracht


# ii)

getVarCov(model_lme, individuals = "M01", type = "conditional")
getVarCov(model_lme, individuals = "F02", type = "conditional")

getVarCov(model_lme, individuals = "M01", type = "marginal")
getVarCov(model_lme, individuals = "F02", type = "marginal")

# -> gleiche Ergebnisse, da:
# 1) die Random Effects ueber alle Personen iid-verteilt sind, und
# 2) fuer die Kovarianzmatrix der Residuen hier gilt: Sigma_i = Sigma
#    (da eine solche Struktur der Residuen hier nicht explizit angenommen wurde)