###########################################################################
############# Loesungen zu Aufgaben aus Tutorium 3 ########################
###########################################################################




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

# Datensatz einlesen
einkommen <- read.table("https://moodle.lmu.de/pluginfile.php/60324/course/section/15118/einkommen.txt", header = TRUE)
# Funktion zur Berechnung der robusten Kovarianzmatrix laden
source("Tutorium3_Funktion_RobusteKovarianz.R")

library(nlme)


### a)

model_1a <- lme(Einkommen ~ Alter, random = ~ 1 | Subjekt, data = einkommen)
summary(model_1a)

# Vergleich der Kovarianzmatrizen
(Cov_robust <- computeCovMatrix(model_1a, robust = TRUE))
(Cov_unrobust <- vcov(model_1a))
# Vergleich der Standardfehler der Betas
(SE_robust <- sqrt(diag(Cov_robust)))
(SE_unrobust <- sqrt(diag(Cov_unrobust)))

# -> Kovarianzmatrizen zeigen gewisse Unterschiede
#    (z.B. ist der robust geschaetzte Standardfehler des Alters deutlich
#     groesser als der unrobust geschaetzte; mit der unrobusten Kovarianz
#     wuerde man den Alterseffekt als signifikant betrachten, mit der robusten
#     Kovarianz nicht!)

# geschaetzte marginale Korrelationsmatrix
cov2cor(getVarCov(model_1a, type = "marginal")[[1]])
# -> in model_1a ist marginale Korrelation zwischen Beobachtungen an selbem Subjekt
#    konstant; keine gute Abbildung der wahren Korrelationsstruktur


### b)

model_1b <- lme(Einkommen ~ Alter, random = ~ 1 | Subjekt, data = einkommen,
                      correlation = corSymm(form = ~ 1 | Subjekt))
summary(model_1b)

# Vergleich der Kovarianzmatrizen
(Cov_robust <- computeCovMatrix(model_1b, robust = TRUE))
(Cov_unrobust <- vcov(model_1b))
# Vergleich der Standardfehler der Betas
(SE_robust <- sqrt(diag(Cov_robust)))
(SE_unrobust <- sqrt(diag(Cov_unrobust)))

# -> nur noch geringe Unterschiede zwischen den Kovarianzen
# -> es sollte auch auffallen, dass der Alterseffekt im model_1b deutlich
#    unterschiedlich zu dem in model_1a ist! Unter Benutzung der korrekten
#    Varianzstruktur ist er ausserdem klar signifikant

# geschaetzte marginale Korrelationsmatrix
cov2cor(getVarCov(model_1b, type = "marginal")[[1]])
# -> bessere Abbildung der wahren Korrelationsstruktur im Modell






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

data(BodyWeight)
head(BodyWeight)


library(nlme)

m1_reml <- lme(weight ~ Diet*Time, random = ~ 1 | Rat, data = BodyWeight)
m1_ml <- lme(weight ~ Diet*Time, random = ~ 1 | Rat, data = BodyWeight, method = "ML")


### a)

# anova.lme:
# Moeglichkeit 1: Man uebergibt nur 1 Modell
#    1.1: Man uebergibt kein L
#    -> F-Tests aller fixed effects (bzw. der durch 'Terms' angegebenen fixed effects)
#    1.2: Man uebergibt zusaetzlich L
#    -> F-Test der linearen Hypothese, die durch L aufgestellt wird
# Moeglichkeit 2: Man uebergibt mehrere Modelle
# -> Vergleich der Modelle mittels LQ-Tests


### b)

# H_0: beta_Diet2 = beta_Diet3
# H_1: beta_Diet2 != beta_Diet3

# Ausformuliert:
# H_0: Der mittlere Gewichtsunterschied bei Time = 0 zu einer Ratte mit Diet 1
#      ist für Ratten mit Diet 2 und Ratten mit Diet 3 gleich hoch

L1 <- c(0,1,-1,0,0,0)
anova(m1_ml, L = L1)

#Alternativ
L2 <- c("Diet2"=1, "Diet3"=-1)
anova(m1_ml, L = L2)

# Testentscheidung: H_0
# -> Die beiden Koeffizienten unterscheiden sich nicht signifikant


### c)

# H_0: beta_Diet2 = beta_Diet3 = beta_Diet2:Time = beta_Diet3:Time = 0
# H_1: Mindestens eines der obigen Betas ist ungleich 0

# Beim LQ-Test von Modellen mit genesteten fixed effects-Strukturen ist zu beachten:
# - ML-Schaetzung in beiden Modellen
# - gleiche random effects in beiden Modellen

m2_ml <- lme(weight ~ Time, random = ~ 1 | Rat, method = "ML",
             data = BodyWeight)
anova(m1_ml, m2_ml)

# Testentscheidung: H_1
# -> signifikanter Unterschied zwischen Modellen
# -> Entscheidung fuer das komplexere Modell m1



### d)

# H_0: d_22 = d_12 = 0
# H_1: d_22 != 0 und/oder d_12 != 0

# Beim LQ-Test von Modellen mit genesteten random effects-Strukturen ist zu beachten:
# - gleiche fixed effects in beiden Modellen
# - Verwendung der Mischverteilung, falls Varianzkomponenten == 0 getestet werden

m3_reml <- update(m1_reml, random = ~ 1 + Time | Rat)
test <- anova(m1_reml, m3_reml)
# Verwenden der Mischverteilung (da d_22 eine Varianz ist und 0 am Rand des Wertebereichs liegt)
(p <- 0.5*(1-pchisq(test$L.Ratio[2],df=1)) + 
  0.5*(1-pchisq(test$L.Ratio[2],df=2)))

# Testentscheidung: H_1
# -> signifikanter Unterschied zwischen Modellen
# -> Entscheidung fuer das komplexere Modell m3