###########################################################################
############# Loesungen zu Aufgaben aus Tutorium 5 ########################
###########################################################################


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

# Datensatz einlesen
zufriedenheit <- read.table("https://moodle.lmu.de/pluginfile.php/60324/course/section/15118/zufriedenheit.txt", header = TRUE, na.strings = ".")

library(nlme)


### a)

# siehe online-Loesung



### b)

#mAIC: 
# - beruht auf marginaler Likelihood
# - gebiased, wenn random effects Struktur in den verglichenen Modellen verschieden ist
#   --> eher zu empfehlen, wenn random effects in beiden Modellen gleich
# - bei Vergleich von Modellen mit verschiedenen fixed effects: ML-Schätzung!

#cAIC:
# - beruht auf konditionaler Likelihood
# - für Vergleich bzgl. der random effects Struktur entwickelt
# - prinzipiell auch zum Vergleich verschiedener fixed effects Strukturen geeignet

# mögliches Vorgehen:
# mAIC: Vergleich von Modellen mit unterschiedlichen fixed effects,
#       bei gleichen random effects
# cAIC: Vergleich von Modellen mit unterschiedlichen random effects,
#       bei gleichen fixed effects


### c)
model1_ml <- lme(Zufriedenheit ~ Einkommen + Jahr, random = ~ 1 | Subjekt,
                 data = zufriedenheit, control = lmeControl(opt = "optim"),
                 method = "ML")

## i)
# Entscheidung ueber Aufnahme eines fixed effects anhand des mAIC
model2_ml <- update(model1_ml, . ~ . + I(Jahr^2))
AIC(model1_ml)
AIC(model2_ml)
# -> Aufnahme des quadratischen Jahr-Effektes sinnvoll,
#    da mAIC(model2) < mAIC(model1)


## ii)
# Entscheidung ueber Aufnahme des Random Slopes anhand dem cAIC
model2_reml <- update(model2_ml, method = "REML")
model3_reml <- update(model2_reml, random = ~ 1 + Einkommen | Subjekt)

#Laden der Funktion cAIC
source("https://moodle.lmu.de/pluginfile.php/60324/course/section/15118/compute_cAIC.R")
compute_cAIC(model2_reml, B = 25)
compute_cAIC(model3_reml, B = 25)
# -> Aufnahme des Random Slopes fuer Einkommen sinnvoll, da cAIC(model3) < cAIC(model2)




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

# Daten laden
library(MASS)
data(epil)
head(epil)

# Binaere Zielgroesse erstellen
epil$y_binaer <- (epil$y > 3)*1


### a)
m_glmm <- glmmPQL(y_binaer ~ age + lbase + trt,
                 random = ~ 1 | subject,
                 family = binomial, data = epil)
summary(m_glmm)

# Konditionale Interpretation der Parameter:
# - treatment:
# Betrachtet man zwei Personen mit gleichem Random Intercept, gleichen
# Auspraegungen in den Kontrollvariablen, jedoch einem unterschiedlichen
# Treatment, so ist die Chance, dass die Person mit Treatment mehr als 3
# Anfaelle in 2 Wochen hat, im Vergleich dazu, dass sie hoechstens 3 Anfaelle
# hat, multiplikativ um exp(-0.789) = 0.45 veraendert im Vergleich zur Chance der
# Person mit Placebo. Die Treatment-Person hat also eine 55% geringere Chance
# mehr als 3 Anfaelle in 2 Wochen zu erleiden.
# - age:
# Wird eine bestimmte Person um 1 Jahr aelter, so steigt die Chance, dass die
# Person mehr als 3 Anfaelle in 2 Wochen hat, im Vergleich dazu, dass sie
# hoechstens 3 Anfaelle hat, multiplikativ um exp(0.035) = 1.04, also um 4%.

# Marginale Interpretationen:
# Nicht zulaessig, da Logit-Link verwendet wird!


### b)
library(gee)
m_gee <- gee(y_binaer ~ age + lbase + trt,
             id = subject, corstr="unstructured",
             family = binomial, data = epil)
summary(m_gee)

# Marginale Interpretation der Parameter:
# - treatment:
# Die mittlere Chance, dass Personen mehr als 3 Anfaelle in 2 Wochen haben,
# im Vergleich dazu, dass sie hoechstens 3 Anfaelle haben, ist bei Personen,
# welche das Treatment bekommen, multiplikativ um den Faktor exp(-0.850) = 0.43
# veraendert im Vergleich zur Chance von Personen, welche ein Placebo bekommen -
# bei konstant gehaltenen uebrigen Kovariablen.
# Personen, welche das Treatment bekommen, haben also eine um 57% geringere
# Chance mehr als 3 Anfaelle in 2 Wochen zu erleiden.
# - age: 0.04
# Die mittlere Chance, dass Personen mehr als 3 Anfaelle in 2 Wochen haben,
# im Vergleich dazu, dass sie hoechstens 3 Anfaelle haben, steigt mit
# jedem Jahr multiplikativ um exp(0.037) = 1.04, also um 4% - bei konstant
# gehaltenen uebrigen Kovariablen.