###########################################################################
############# Loesungen zu Aufgaben aus Tutorium 4 ########################
###########################################################################



# Daten laden -------------------------------------------------------------

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

library(nlme)

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

source("https://moodle.lmu.de/pluginfile.php/60324/course/section/15118/Tutorium4_TransformierteResiduen.R")


### a)
model <- lme(Zufriedenheit ~ Einkommen + Jahr, random = ~ 1 | Subjekt, data = zufriedenheit)
plot(model)
# Was kann man aus der Grafik schlussfolgern?
# -> Nichts! Residuen sind korreliert und somit nicht interpretierbar!


### b)
# ...


### c)
r.star(model, plot = TRUE)
# Was kann man aus der Grafik schlussfolgern?
# -> Modell falsch spezifiert! Quadratischer Zeiteffekt waere evtl sinnvoll



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

###a)
?corClasses
# Compound symmetry: corCompSymm()
# Unstructured correlation: corSymm()
# AR(1)-Prozess: corAR1 (bei stetiger Zeitvariable: corCAR1)
# -> allgemeinere AR-Strukturen mittels corARMA umsetzbar
# Serielle exp-Korrelation: corExp()
# Serielle Gauss-Korrelation: corGaus()

###b)
lm1 <- lm(Zufriedenheit ~ Jahr + I(Jahr^2) + Einkommen+ Anzahl_Freunde, data = zufriedenheit)
summary(lm1)

#Abge�nderte Funktion aus �bung:
r <- resid(lm1) # get residuals from linear model

t<- zufriedenheit$Jahr # get time points
id <- zufriedenheit$Subjekt # get IDs

half.squared<- outer(r, r, FUN = function(r1, r2) (r1 - r2)^2/2) # half squared residuals for all observations
distances <- outer(t, t, FUN = function(t1, t2) abs(t1 - t2)) # absolute distances for all observations
samei <- outer(id, id, FUN = function(i1, i2) (i1 == i2)) # get index matrix specifying if same subject
nondiag <- 1 - diag(length(r)) # matrix with 0 on the diagonal and 1 for non-diagonal values as 0 on diagonal
u <- as.vector(distances[samei & nondiag]) 
v <- as.vector(half.squared[samei & nondiag])
totvar <- mean(half.squared[!samei]) # total variance = \sigma^2 + d^2+ \tau^2

# smooth and plot
smooth_mean <- lowess(u,v,iter=0) # note: iter = 0 such that not the median is approximated, see ?lowess
par(mar=c(5.1, 5.1, 4.1, 2.1))
plot(u, v, pch = 16, cex = 0.5, col = 8,ylim=c(0,totvar),xlim=c(0,max(u)+0.5),
     xlab="Abstand Jahre(u)",ylab=expression(paste(hat("v"),"(u)")),main="Empirsches Semi-Variogram")
lines(smooth_mean, lwd = 2,t="b")
abline(h = totvar, col = 2, lwd = 2, lty = 2)
# keine klare Struktur erkennbar --> serielle Korrelation vermutlich nicht n�tig


# Aufgabe 3 ---------------------------------------------------------------

# Hausman-Test-Funktion laden
source("https://moodle.lmu.de/pluginfile.php/60324/course/section/15118/Funktion_HausmanTest.R")

#Fit LMM
model_lme <- lme(Zufriedenheit ~ Jahr + I(Jahr^2) + Anzahl_Freunde + Einkommen,
               random = ~ 1 + Einkommen | Subjekt, data = zufriedenheit)

# Random effects assumption (REA) erfuellt?
# REA fuer alle metrischen Kovariablen testen
# (fuer kategoriale Kovariablen gibt es keinen Test, da diese im fixed
#  effects Modell nicht schaetzbar sind)
### Hausman-Test
# Beachten: Da neben dem Random Intercept auch ein Random Slope fuer Einkommen
#           im Modell ist muss im Fixed Effects Modell neben dem Haupteffekt
#           fuer die Subjekte auch eine Interaktion Subjekt:Einkommen aufgenommen
#           werden!
model_re <- model_lme
zufriedenheit$Subjekt <- factor(zufriedenheit$Subjekt)
model_fe <- lm(Zufriedenheit ~ Jahr + I(Jahr^2) + Anzahl_Freunde +
                 Einkommen*Subjekt, data = zufriedenheit)
# Verwendung der selbst geschriebenen Hausman-Funktion
fixed_effects <- c("Jahr","I(Jahr^2)","Anzahl_Freunde","Einkommen")
hausman.test(model_re, model_fe, fixed_effects)
# -> Hausman-Test wird nicht abgelehnt
#    -> eine Verletzung der REA kann nicht gesichert nachgewiesen werden
#       (die REA koennte aber trotzdem verletzt sein! Deswegen auch inhaltlich
#        ueberlegen, ob vielleicht etwas fuer eine Verletzung der Annahme spricht
#        (Bsp.: REA waere verletzt, wenn es einen Confounder gaebe, der nicht
#         im Modell ist, aber sowohl mit der Zufriedenheit als auch der Anzahl
#         an Freunden korreliert))
# -> Generell problematisch: Das fixed effects-Modell enthaelt nun 103
#    Parameter bei nur 600 Beobachtungen!