# ------------------------------------------------------------------------------
# -------- Übung zur Vorlesung Generalisierte Regressionsmodelle ---------------
# -------- R Code / Blatt 12                                     ---------------
# ------------------------------------------------------------------------------


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

# Einlesen der Daten
fakesoep <- read.table("fakesoep.dat",header=T)
# Package mgcv laden
library(mgcv)

# (a)

# Datensatz aufteilen
fakesoepW <- fakesoep[fakesoep$geschl == 1,]
fakesoepM <- fakesoep[fakesoep$geschl == 0,]

# In Blatt 11 wurde ein Additives Modell (AM) betrachtet. Nun haben wir einen
# Gamma-verteilten Response, also brauchen wir Generalisierte Additive Modelle
# (GAM)!

# In der gam()-Funktion können dieselben Link-Funktionen sowie Families
# spezifiziert werden wie in glm().
modelW <- gam(beink ~ s(groesse,bs="cr") + s(alter,bs="cr") + s(dauer,bs="cr")
+ verh + deutsch + abitur, family=Gamma(link="log"), data=fakesoepW)
modelM <- gam(beink ~ s(groesse,bs="cr") + s(alter,bs="cr") + s(dauer,bs="cr")
+ verh + deutsch + abitur, family=Gamma(link="log"), data=fakesoepM)

# Ergebnis betrachten
summary(modelW)
summary(modelM)
# Nach den geschätzen Modellen wirkt sich bei Frauen eine Hochzeit im Mittel
# offenbar negativ auf das Einkommen aus, bei Männern hingegen positiv. Der
# Einfluss der Variablen Abitur ist wie erwartet.

# Die glatten Kurven können durch Plots dargestellt werden.
# Der Einfluss von Körpergröße auf Einkommen erscheint (bei Frauen und
# Männern) linear.
plot(modelW)
plot(modelM)


# Prognose für eine verheiratete, 40 jährige, 1,70m große deutsche Frau
# ohne Abitur, die seit 15 Jahren demselben Berieb angehört.
# Daten als data.frame spezifizieren
newdata <- data.frame(groesse=170,alter=40,dauer=15,verh=1,deutsch=1,abitur=0)
# predict() Funktion verwenden
predict(modelW,newdata=newdata,type="response")




# (b)

# Um die beiden obigen Modelle in einem Modell zusammenzufassen,
# muss die Variable Geschlecht sowie Interaktionen aller anderen Variablen
# mit Geschlecht ins Modell aufgenommen werden. Hierfür müssen zunächst
# Geschlechts-Dummies erzeugt werden
geschl1 <- fakesoep$geschl
geschl0 <- 1-geschl1

modelAll <- gam(beink ~ s(groesse,bs="cr",by=geschl0) +
s(alter,bs="cr",by=geschl0) + s(dauer,bs="cr",by=geschl0) +
s(groesse,bs="cr",by=geschl1) +
s(alter,bs="cr",by=geschl1) + s(dauer,bs="cr",by=geschl1) +
geschl + verh + deutsch + abitur + verh:geschl + deutsch:geschl + abitur:geschl,
family=Gamma(link="log"),data=fakesoep)

# An den Parameter-Schätzern der parametrisch aufgenommenen Variablen erkennt
# man die Äquivalenz der Modelle (abgesehen von kleineren Ungenauigkeiten)
summary(modelAll)
summary(modelW)
summary(modelM)

# Auch das Ergebnis der glatten Modellierung ist ähnlich.
plot(modelAll)

## gemeinsames Modell ohne Trennung nach Geschlecht:
modelWM <- gam(beink ~ s(groesse,bs="cr") + s(alter,bs="cr") + s(dauer,bs="cr")
+ verh + deutsch + abitur, family=Gamma(link="log"), data=fakesoep)

## Vergleich
par(mfrow=c(3,3))

plot(modelW, main="Geschl=W")
plot(modelM, main="Geschl=M")
# für den gesamten Datensatz (Frauen und Männer zusammen)
plot(modelWM, main="Geschl=W&M")

## Vergleich des Effekts der Größe zwischen gemeinsamen Modell und den beiden
## nach Geschlechtern getrennten Modellen:
par(mfrow=c(2,2))
plot(modelW, select=1)
plot(modelM, select=1)
plot(modelWM, select=1)

## der Effekt der Größe auf das Einkommen ist sowohl für Männer als auch für
## Frauen eine Gerade, aber das allgemeine Niveau dieses Größeneffekts unter-
## scheidet sich zwischen den Geschlechtern, er ist bei Männern höher.
## Diese Niveauunterschiede sieht man in den einzelnen Plots wegen der Zentrierung
## der Kurven um Null herum nicht.
## Die nonparametrisch geschätzte Kurve in der dritten Grafik zeigt, wie
## der Schätzer versucht, von der niedrigeren zur höheren Gerade zu springen.
## In den extremen Größenbereichen, in denen fast nur eines der Geschlechter
## vorkommt, reproduziert der Schätzer den Verlauf der jeweiligen Gerade.
## Fazit: geschätzte Kurven in (G)AMs sind empfindlich gegenüber Interaktionen
## mit übrigen Einflussgrößen. (Hier dem Geschlecht.)


# (c)

# Eine Alternative zur Gamma-Verteilung ist die inverse Gauß-Verteilung
# (vgl. Blatt 4)
modelAllinvG <- gam(beink ~ s(groesse,bs="cr",by=geschl0) +
s(alter,bs="cr",by=geschl0) + s(dauer,bs="cr",by=geschl0) +
s(groesse,bs="cr",by=geschl1) +
s(alter,bs="cr",by=geschl1) + s(dauer,bs="cr",by=geschl1) +
geschl + verh + deutsch + abitur + verh:geschl + deutsch:geschl + abitur:geschl,
family=inverse.gaussian(link="log"),data=fakesoep)

# Die Ergebnisse sind sehr ähnlich.
summary(modelAllinvG)
plot(modelAllinvG)


# Dignostik-Plot können mittel gam.check() erzeugt werden.
gam.check(modelAll)
gam.check(modelAllinvG)
# Offenbar ist die Anpassung der (Devianz-)Residuen an die Normalverteilung
# im Gamma-Modell besser, weshalb jenes dem Invers-Gauß-Modell vorzuziehen
# ist.


# (d)

# Wir betrachten aus Gründen der Übersichtlichkeit nun Frauen und Männer
# wieder getrennt und schätzen zum Vergleich GLMs (vgl. Blatt 4).
glmW <- gam(beink ~ groesse + alter + dauer
+ verh + deutsch + abitur,family=Gamma(link=log),data=fakesoepW)
glmM <- gam(beink ~ groesse + alter + dauer
+ verh + deutsch + abitur,family=Gamma(link=log),data=fakesoepM)

# Der Vergleich mit den GAMs aus (a) kann mittel LR-Tests erfolgen
anova(modelW,glmW,test="F")
anova(modelM,glmM,test="F")
# Offenbar sind GAMs in beiden Fällen zur Modellierung wesentlich besser
# geeignet.

# Insbesondere bei der Kovariablen Körpergröße stellt sich aber die Frage,
# ob der Einfluss nonparametrisch modelliert werden sollte.
# Nach den Plots in (a) zu urteilen, ist die glatte Modellierung hier
# nicht notwendig.

modelW_k <- gam(beink ~ groesse + s(alter,bs="cr") + s(dauer,bs="cr")
+ verh + deutsch + abitur,family=Gamma(link=log),data=fakesoepW)
modelM_k <- gam(beink ~ groesse + s(alter,bs="cr") + s(dauer,bs="cr")
+ verh + deutsch + abitur,family=Gamma(link=log),data=fakesoepM)

# Paradoxerweise wird aber bei Verwendung von Tests die nonparametrische
# Modellierung der parametrischen vorgezogen, obwohl (bzw. weil) der
# geschätzte Einfluss (fast) linear ist.
anova(modelW,modelW_k,test="F")
anova(modelM,modelM_k,test="F")




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

# Veranschaulichung von lokaler linearer Regression

# Epanechnikov-Kern
epan <- function(x,hvalue)
  {
    ifelse(abs(x/hvalue)<=1,(3/4)*(1-(x/hvalue)^2), 0)/hvalue
  }

xseq <- seq(-3,3,by=0.01)
plot(xseq,epan(xseq,0.5),type="l",xlab="x",ylab="K(x/h)/h",lwd=2)
lines(xseq,epan(xseq,1),col=2,lwd=2)
lines(xseq,epan(xseq,2),col=3,lwd=2)


# Daten
x <- 3*runif(100)
fx <- function(x){2 - (x-1.5)^2 + 3*sqrt(x)}
y <- fx(x) + rnorm(100)
plot(x,y,xlab="x",ylab="y")
xseq <- seq(0,3,by=0.01)
lines(xseq,fx(xseq),col=3)


# neues x
xnew <- 1
abline(v=xnew,lty=2)

# Gewichte
w <- epan(x-xnew,0.5)
points(x,w,type="h")

# lokale lineare Regression
xregr <- x-xnew
lmlokal <- lm(y ~ xregr, weights=w)
points(xnew,lmlokal$coef[1],col=2,lwd=2)
lines(xnew+c(-0.5,0.5),lmlokal$coef[1]+ lmlokal$coef[2]*c(-0.5,0.5))


# für eine Kurve Grid von neuen x-Werten
xnewseq <- seq(0,3,by=0.01)
fit <- numeric(length(xnewseq))
i <- 1
for (xnew in xnewseq)
  {
    w <- epan(x-xnew,0.5)
    xregr <- x-xnew
    lmlokal <- lm(y ~ xregr, weights=w)
    fit[i] <- lmlokal$coef[1]
    i <- i+1
  }

plot(x,y,xlab="x",ylab="y")
lines(xnewseq,fx(xnewseq),col=3)
lines(xnewseq,fit,col=2)