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



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


# Einlesen der Daten
fakesoep <- read.table("fakesoep.dat",header=T)
attach(fakesoep)


# a)

# Man untersucht die Gestalt der Verteilung des Bruttoeinkommens empirisch,
# z.B. mit einem Histogramm, Kerndichteschätzer oder Boxplot.
par(mfrow=c(1,3))
hist(beink)
plot(density(beink))
boxplot(beink)

# Die Variable beink nimmt positive Werte auf metrischem Skalenniveau an, die
# (etwas) linkssteil bzw. rechtsschief verteilt sind. Dafür spricht auch die
# 5-Zahlen-Zusammenfassung aus

summary(beink) # Median < Mean

# d)

xseq <- seq(0.01,4,by=0.01)

# ein Beispiel 
a <- 2
plot(xseq,dgamma(xseq,shape=a,scale=1/a),ylab="",xlab="x",
main=paste("shape =",a),type="l")

# mit Funktion manipulate 
library(manipulate)

dens_gamma <- function(x,a){
  plot(xseq,dgamma(xseq,shape=a,scale=1/a),ylab="",xlab="x",
       main=paste("shape =",a),type="l")
}
manipulate(dens_gamma(x=xseq,a),a=slider(0,5,step=0.1))

# e)

# Natürlicher Link
gammaNat <- glm(beink ~ groesse + alter + dauer + verh + geschl + deutsch + abitur,
family=Gamma())
summary(gammaNat)

# Log Link
gammaLog <- glm(beink ~ groesse +alter + dauer + verh + geschl + deutsch + abitur,
family=Gamma(link=log))
summary(gammaLog)

# Vergleich der Koeffizienten 
cbind(gammaNat$coefficients, gammaLog$coefficients)

# Für eine tiefergehende inhaltliche Interpretation der Parameter wäre die
# Einbeziehung von Interaktionstermen (z.B. verh und geschl) wünschenswert.
# also etwa (Zugabe, keine Aufgabenstellung aber interessant!):
gammaLog2 <- update(gammaLog, . ~ . + verh:geschl)
summary(gammaLog2)

# (f)

par(mfrow=c(1,2))

# Plot: Beobachtete Werte gegen Geschätzte Werte
plot(beink,gammaLog$fitted,xlab="observed",ylab="fitted")
abline(0,1)
# Zusatz
points(beink[geschl==1 & abitur==0],gammaLog$fitted[geschl==1 & abitur==0],col=2)

# Plot: Geschätzte lineare Prädiktoren gegen Residuen 
plot(gammaLog$linear.predictors,residuals(gammaLog),xlab="linear predictors",ylab="residuals")
# Zusatz
points(gammaLog$linear.predictors[geschl==1 & abitur==0],gammaLog$residuals[geschl==1 & abitur==0],col=2)

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

# b) 

library("statmod")
xseq <- seq(0.01,4,by=0.01)

# ein Beispiel 
lambda <- 2
plot(xseq,dinvgauss(xseq,mean=1,dispersion=1/lambda),type="l",ylab="",xlab="x",
     main=paste("lambda =",lambda))

# mit Funktion manipulate 
library(manipulate)

dens_invgauss <- function(x,lambda){
  plot(xseq,dinvgauss(xseq,mean=1,dispersion=1/lambda),type="l",ylab="",xlab="x",
       main=paste("lambda =",lambda))
}
manipulate(dens_invgauss(x=xseq,lambda),lambda=slider(0,5,step=0.1))

# (c)

# Modell mit natürlichem Link:
help(inverse.gaussian)
invgNat <- glm(beink ~ groesse + alter + dauer + verh + geschl + deutsch + abitur,
family=inverse.gaussian())

# Der Schätzer mit natürlichem Link kann nicht berechnet werden, da R hier keine
# Startparameter für den Schätzalgorithmus findet, die die durch den natürlichen
# Link induzierten Restriktionen erfüllen.


# Modell mit log-Link:
invgLog <- glm(beink ~ groesse + alter + dauer + verh + geschl + deutsch + abitur,
family=inverse.gaussian(link=log))
#Genau wie im Gamma-Modell unterliegt der Parametervektor auch bei log-Link 
#keinen Restriktionen

# Summary
summary(invgLog)

# zum Vergleich
summary(gammaLog)


# Übergibt man als Startwerte die Parameterschätzer aus dem Gamma-Modell,
# beginnt R zwar mit der Berechnung, der Algorithmus konvergiert allerdings nicht:
invgNat <- glm(beink ~ groesse + alter + dauer + verh + geschl + deutsch + abitur,
family=inverse.gaussian(),start=gammaNat$coef)

warnings()