# @Author: andreas.bender@stat.uni-muenchen.de
# @Date:   2015-05-21 11:57:12
# @Last Modified by:   andreas.bender@stat.uni-muenchen.de
# @Last Modified time: 2015-05-29 21:35:34

library(effects)
library(ggplot2)
theme_set(theme_bw())

## ----Aufgabe2------------------------------------------------------------
## angelehnt an Everitt, Hothorn: "A Handbook of Statistical Analyses using R", Ch. 4

## einlesen der Daten
weight <- read.table("http://www.statistik.lmu.de/~abender/unterlagen/limo14/Datensatz4.txt", 
  header = TRUE)

# graphische Darstellung: Abweichung der Mittelwerte vom Gesamtmittelwert
# (Betrachtung der einzelnen Faktoren, nicht ihrer Kombination)
plot.design(weight)
mean(weight$weightgain)
tapply(weight$weightgain, weight$source, mean)
#aggregate(weightgain ~ source, data=weight, mean)
tapply(weight$weightgain, weight$amount, mean)
#aggregate(weightgain ~ source, data=weight, mean)

# Berechnung von Mittelwert und Standardabweichung fuer jede Kombination
tapply(weight$weightgain, list(weight$source, weight$amount), mean)
tapply(weight$weightgain, list(weight$source, weight$amount), sd)
# aggregate(weightgain ~ source*amount, data=weight, mean)
# aggregate(weightgain ~ source*amount, data=weight, sd)


## zweifaktorielle Anova (mit Interaktion)
# per default faktoren referenz-kodiert
wg.aov <- aov(weightgain ~ source + amount + source:amount, data=weight)

# Designmatrix
head(model.matrix(wg.aov))
cbind(model.matrix(wg.aov), weight[, 1:2])

summary(wg.aov)

# graphische Darstellung: Interaktion
# interaction.plot(weight$amount, weight$source, weight$weightgain, legend=FALSE, 
#   las = 1)
# legend(1.5, 95, legend = levels(weight$source), title = "weight$source", 
# 	lty = c(2,1), bty="n")
ggplot(aggregate(weightgain ~ source + amount, mean, data=weight)) + 
  aes(x = amount, y = weightgain, lty=source, group=source) + 
  geom_point() + 
  geom_line()

# Parameterschaetzer
coef(wg.aov)
coef(lm.sa <- lm(weightgain ~ source * amount, data = weight))
## alternative visualisierung
library(effects)
plot(allEffects(lm.sa))

## Modell mit Effektkodierung
wg.aov.effect <- aov(weightgain ~ source + amount + source:amount, data=weight, 
    contrasts=list(source=contr.sum, amount=contr.sum))

# Designmatrix
head(model.matrix(wg.aov.effect))

summary(wg.aov.effect)
# (keine Veraenderung)

# Parameterschaetzer
coef(wg.aov.effect)
## mit lm
lm.effect <- lm(weightgain ~ source*amount, data = weight,
  contrasts = list(source = contr.sum, amount = contr.sum))
coef(lm.effect)



## --------- Aufgabe 4 ----------------------------------------------------
## ----Aufgabe4-a----------------------------------------------------------
library(faraway)
teengamb$sexF <- factor(teengamb$sex, levels = c(0,1), labels = c("m", "w"))
lm.gamb <- lm(gamble ~ income + sexF + status + verbal, data = teengamb)
coef.gamb <- coef(lm.gamb)
summary(lm.gamb)

# visualisierung durch Koeffizienten-Plot
library(sjPlot)
sjp.lm(lm.gamb)

# Geschlechtseffekt ist betragsmäßig deutlich größer und hat eine höhere Varianz. 
# Aber: Parameterschätzer beziehen sich auf Veränderung um eine Einheit, d.h. 
# auch Minima/Maxima bzw. Streuung (Standardabweichung) der Kovariablen relevant. 
# Vgl. Blatt05/Standardisierung

##### Zusatz : 
## Effektkodierung
lm.gamb.effect <- update(lm.gamb, contrasts = list(sexF = "contr.sum"))
### Vorsicht bei Interpretation, Geschlechtseffekt bezieht sich nun auf die 
## Kategorie-Ausprägung "m" (männlich)
summary(lm.gamb.effect)
### Erwartungswerte in den jeweiligen Gruppen (m/w), wenn alle Variablen auf
## Mittelwerte gesetzt werden
df.effect <- aggregate(cbind(income, verbal, status) ~ sexF, mean, 
 data = teengamb)
predict(lm.gamb.effect, newdata = df.effect)
## entspricht hier den Mittelwerten der Zielvariable in den beiden Gruppen
aggregate(gamble ~ sexF, mean, data = teengamb)
sum(coef(lm.gamb.effect)*c(1,0,0,0,0)) - sum(coef(lm.gamb.effect)*c(1,0,1,0,0))

#Interpretation:
#Im Vergleich zum mittleren Erwartungswert (wenn alle metrischen Kovars Null), 
#erwartet man bei männlichen Teenagern im Durchschnitt um 
coef(lm.gamb.effect)["sexF1"]
# erhöhte Glückspielausgaben.
#### Ende Zusatz

## Zweifachinteraktionen der metrischen Variablen mit sex:
## ----Aufgabe4-b----------------------------------------------------------
lm.gamb.inter <- lm(gamble ~ (income + status + verbal)*sexF, data = teengamb, 
  x = TRUE)
summary(lm.gamb.inter)
## sub modell
lm.gamb.sub <- lm(gamble ~ income*sexF, data = teengamb)
summary(lm.gamb.sub)
# man beachte, dass im Submodell der Geschlechtseffekt nicht signifikant ist, 
# allerdings werden Haupteffekte (sexF) normalerweise im Modell belassen, wenn 
# hierarchisch untergeordnete Terme (sexF:income) im Modell enthalten sind
anova(lm.gamb.inter, lm.gamb.sub)
## H_0 kann nicht verworfen werden


# Graphische Darstellung der Interaktion: 
# Dargestellt ist der Einfluss des Einkommens gegeben Geschlecht. 
## ----Aufgabe2-c
# mit Hilfe des effects packages
library(effects)
effect.inter <- effect("income*sexF", mod = lm.gamb.sub)
plot(effect.inter)
plot(allEffects(lm.gamb.sub))

## mit ggplot2
ggplot(teengamb, aes(x = income, y = gamble)) + 
  geom_point() + facet_wrap(~ sexF) + 
  geom_smooth(method="lm")
# oder 
ggplot(teengamb, aes(x = income, y = gamble, pch = sexF)) + 
  geom_point() + 
  geom_smooth(method="lm")


# Verläufe nicht parallel, d.h. grafischer Hinweis, dass eine Interaktion 
# notwendig ist. 
# Deutlich niedrigere Ausgaben bei Frauen und offenbar unabhängig vom Einkommen
# Bei Männern ähnlich niedrige Ausgaben wie Frauen bei niedrigem Einkommen
# (Wenn man kein Geld hat, kann man auch keins ausgeben), aber steigende Ausgaben 
# mit steigendem Einkommen, im Ggs. zu Frauen
# Haupteffekte isoliert oft schwierig bis unmöglich zu interpretieren
# Stattdessen: 
## Erwartugswertvergleiche: 
(effect.sub <- summary(effect.inter)$effect)

## Im Vgl. zu einem weiblichen Teenager mit einem Einkommen 
# von $2$Pfund/Woche hat ein männlicher Teenager mit gleichem Einkommen im 
# Durchschnitt 
effect.sub[1,1] - effect.sub[1,2]
# Pfund/Jahr höhere Ausgaben für Glückspiel.
## Im Vgl. zu einem weiblichen Teenager mit einem Einkommen 
# von $10$Pfund/Woche hat ein männlicher Teenager mit gleichem Einkommen im 
# Durchschnitt 
effect.sub[5,1] - effect.sub[5,2]
# Pfund/Jahr höhere Ausgaben für Glückspiel.
## Männlich, 15Pfund Einkommen vs. weiblich, 10 Pfund Einkommen
# $E(y|sex=``m'', x_income = 15) - E(y|sex=``w'', x_income} = 10) =
sum(coef(lm.gamb.sub)[c(2, 2, 3, 4)]*c(15, -10, -1, -10))


# Manchmal sieht man zweifach-Interaktionen auch schon bei visueller Betrachtung.
ggplot(teengamb, aes(x = income, y = gamble, pch = sexF, color = sexF)) + 
  geom_point()