########################################################################### ############# Loesungen zu Aufgaben aus Tutorium 2 ######################## ########################################################################### # Aufgabe 1 --------------------------------------------------------------- ### i) Daten einlesen --------------------------------------------------------- # Datensatz einlesen zufriedenheit <- read.table("https://moodle.lmu.de/pluginfile.php/60324/course/section/15118/zufriedenheit.txt", header = TRUE) ### ii) Untersuchung des Zeiteffekts # Zielgroesse ueber die Zeit plotten plot(zufriedenheit$Jahr, zufriedenheit$Zufriedenheit, col = "gray") lo <- loess(Zufriedenheit ~ Jahr, data = zufriedenheit) lines(2005:2016, predict(lo, newdata = data.frame("Jahr" = 2005:2016))) # -> Einfluss wohl eher nicht linear; evtl quadratisch ### iii) Modellfit # nlme-Paket fuer lme-Funktion laden library(nlme) ## 1.Versuch: Ausprobieren eines einfachen Modells mit linearem Zeiteffekt # und einem Random Intercept model_1 <- lme(Zufriedenheit ~ Jahr + Geschlecht + Anzahl_Freunde + Einkommen, random = ~ 1| Subjekt, data = zufriedenheit) summary(model_1) # -> linearer Zeiteffekt signifikant ### iv) Residuenplot plot(model_1) # -> Residuen weisen noch klare Struktur auf! Verwerfe 1.Versuch! ### v) Modellfit 2. Versuch ## 2.Versuch: 1.Versuch + quadratischen Term fuer die Zeit # (linearer Effekt bleibt aber trotzdem im Modell, # da Zeiteffekt nun als quadratisches Polynom aufgenommen werden soll!) model_2 <- update(model_1, . ~ . + I(Jahr^2)) summary(model_2) plot(model_2) # -> Quadratischer Term signifikant und Struktur aus Resiuden verschwunden! ### vi) ICC: Ist Random Intercept sinnvoll? (var_ri <- as.numeric(getVarCov(model_2))) (var_error <- model_2$sigma^2) (icc <- var_ri/(var_ri + var_error)) # -> ICC = 0.44: Random Intercept erklaert hier fast die Haelfte der Fehlervarianz! # -> Hier wuerde man den RI also klar im Modell behalten! ### vii) Sind Random slopes sinnvoll? # Ueberpruefung anhand einzelner Modelle pro Subjekt # (zeitkonstante Variable 'Geschlecht' kann hier nicht mitaufgenommen werden; # macht aber auch nichts, da ein Random Slope von zeitkonstanten Variablen # eh nicht sinnvoll sind, da diese Effekte immer in den Random Intercept wandern) lmlist <- lmList(Zufriedenheit ~ Jahr + I(Jahr^2) + Anzahl_Freunde + Einkommen | Subjekt, data = zufriedenheit) plot(intervals(lmlist)) # Mögliche subjektive Interpretation: # -> Die Grafik spricht...: # - fuer einen Random Intercept # - gegen einen Random Slope fuer die Zeitvariable # - gegen einen Random Slope bzgl. Anzahl_Freunde # - fuer einen Random Slope bzgl. Einkommen ### viii) # 3.Versuch: 2.Versuch + Random Slopes fuer Anzahl_Freunde und Einkommen model_3 <- update(model_2, random = ~ 1 + Einkommen | Subjekt) summary(model_3) plot(model_3) # Residuen weisen keine Struktur auf #Interpretation Einkommen: Steigt das Einkommen um 1000 Euro, so steigt die #Lebenszufriedenheit durchschnittlich um 0.12 Punkte (c.p.). ### ix) Annahme korrelierter Random Effects sinnvoll? cov2cor(getVarCov(model_3)) # -> mit -0.176 nur sehr schwache Korrelation der Random Effects; # Hier koennen wir die Annahme korrelierter Random Effects beibehalten, # da unser Datensatz gross genug ist, um den durch die Annahme # zusaetzlich zu schaetzenden Varianzparameter mit zu berechnen; # -> Bei einer kleinen Stichprobe koennte man auch unkorrelierte Random # Effects annehmen, um 1 Parameter weniger schaetzen zu muessen # Die Annahme kann dann durch 'pdDiag' in das lme-Modell eingebaut werden # (Beachten: dazu muss der Datensatz vorher als groupedData-Objekt # formatiert werden, damit die lme-Funktion weiss, was die ID-Variable ist) ### x) #Annahmen: #mean(zufriedenheit$Jahr)=2010.5 #mean(zufriedenheit$Anzahl_Freunde)=5.518 #Geschlecht:0 (weiblich) #Betrachte Subjekt 1,2,4 #globale Prädiktion x11() plot(c(0,16), c(fixef(model_3)["(Intercept)"] + mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"], fixef(model_3)["(Intercept)"] + mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"]+ 16*fixef(model_3)["Einkommen"]), type="l", xlab="Einkommen", ylab="Zufriedenheit") #Prädiktion Subjekt 1 lines(c(0,16),c(fixef(model_3)["(Intercept)"]+ranef(model_3)[1,"(Intercept)"]+ mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"] , fixef(model_3)["(Intercept)"]+ranef(model_3)[1,"(Intercept)"]+ mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"]+ 16*fixef(model_3)["Einkommen"]+16*ranef(model_3)[1,"Einkommen"]), col="red") #Prädiktion Subjekt 2 lines(c(0,16),c(fixef(model_3)["(Intercept)"]+ranef(model_3)[2,"(Intercept)"]+ mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"] , fixef(model_3)["(Intercept)"]+ranef(model_3)[2,"(Intercept)"]+ mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"]+ 16*fixef(model_3)["Einkommen"]+16*ranef(model_3)[2,"Einkommen"]), col="red") #Prädiktion Subjekt 4 lines(c(0,16),c(fixef(model_3)["(Intercept)"]+ranef(model_3)[4,"(Intercept)"]+ mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"] , fixef(model_3)["(Intercept)"]+ranef(model_3)[4,"(Intercept)"]+ mean(zufriedenheit$Jahr)*fixef(model_3)["Jahr"]+ mean(zufriedenheit$Jahr)^2*fixef(model_3)["I(Jahr^2)"]+ mean(zufriedenheit$Anzahl_Freunde)*fixef(model_3)["Anzahl_Freunde"]+ 16*fixef(model_3)["Einkommen"]+16*ranef(model_3)[4,"Einkommen"]), col="red") # Legende legend("topright", lty = c(1,2), col = c("black", "red"), lwd = 3, legend = c("globaler Effekt","Subjektspezifische Effekte"))