# ------------------------------------------------------------------------------
# -------- Ü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


# b)

# Vorsicht! Soeben wurde die Randverteilung von beink untersucht, für 
# GLM ist aber die bedingte Verteilung geg. die Prädiktoren entscheidend!

# Nicht-negativer Träger und linkssteil trifft aber mit Sicherheit bzw. 
# vermutlich auch auf bedingte Verteilung zu, da Einkommensdaten vorliegen.

# Gamma-Verteilung:
# - nicht negativ,
# - linkssteil,
# - sehr flexibel durch shape-Parameter.


# c)
# Dichte zeichnen
# help(dgamma)

xseq <- seq(0.01,4,by=0.01)
par(mfrow=c(2,2))
a <- 0.5
plot(xseq,dgamma(xseq,shape=a,scale=1/a),ylab="",xlab="x",
main=paste("shape =",a),type="l")
a <- 1
plot(xseq,dgamma(xseq,shape=a,scale=1/a),ylab="",xlab="x",
main=paste("shape =",a),type="l")
a <- 2
plot(xseq,dgamma(xseq,shape=a,scale=1/a),ylab="",xlab="x",
main=paste("shape =",a),type="l")
a <- 5
plot(xseq,dgamma(xseq,shape=a,scale=1/a),ylab="",xlab="x",
main=paste("shape =",a),type="l")


# d)

# 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)

# Da man für beide Modelle die gleiche Verteilungsannahme tätigt, kann man die
# Anpassungsgüte anhand der Devianz, Dev=-2\phi\sum_i(l_i(\hat\mu_i)-l_i(y_i)),
# vergleichen. Da das Modell mit log-Link mit 502.50 den kleineren Wert aufweist
# (nat. Link: 507.22), ist dieses zu bevorzugen. Vergleich über AIC führt 
# (hier selbstverständlich) zum selben Ergebnis.

# Man beachte, dass hier der Dispersionsparameter gemäß der Formel
# \hat\phi=1/(n-p)\sum_i(y_i-\hat\mu_i)^2/\hat\mu_i^2 mitgeschätzt wird.

help(Gamma)
# Ein weiterer Vorteil des log-Links ist, dass die Erwartungswerte aufgrund von
# \mu=exp(\eta) nur positive Werte annehmen können. Das ist beim natürlichen
# (inversen) Link mit \mu=\eta^{-1} nicht der Fall, d.h. gegebenenfalls
# werden Restriktionen an den Parametervektor \beta erforderlich,
# um \mu>0 zu garantieren.

# Vorsicht bei der Interpretation der Parameter: bei Verwenden des inversen
# Links gilt \mu=\eta^{-1}, d.h. negative Werte für Koeffizienten von \beta
# indizieren einen steigenden Erwartungswert.

cbind(gammaNat$coefficients, gammaLog$coefficients)
# Die Effekte der Kovariablen auf das Bruttoeinkommen weisen in beiden
# Modellen (in fast allen Fällen) in die gleiche Richtung (positiv: groesse, 
# alter, dauer, deutsch, abitur; negativ: geschl), zu beachten ist, dass der 
# Effekt von verh beim natürlichen Link in die andere Richtung weist als beim 
# log-Link. Der Effekt ist allerdings auch nicht signifikant. Das selbe gilt für 
# alter und deutsche Staatsangehörigkeit. Bzgl. Tests sei aber auf die kommende  
# Übung verwiesen.

# 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)


# (e)

# Graphik-Parameter zurück
par(mfrow=c(1,1))

# Plot: Beobachtete 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)


# mit plot.glm
plot(gammaLog, which=1:5)

# 5 Plots werden von generischer plot-Fkt. ausgegeben (Diagnostik-Maße):
# - Residuen gegen gefittete Werte (mit Glätter)
# - QQ-Plot von standardisierten Devianzresiduen gegen theoret. Quantile
#   der N(0,1) (vgl. später)
# - sqrt(stand. Devianzresiduen) gegen gefittete Werte (mit Glätter)
# - Cook's Distance (vgl. später)
# - Residual gegen Leverage. Veranschaulicht Konturen mit gleichen Werten 
#   von Cook's Distance (vgl. später)

# extreme Beobachtungen werden markiert


# (f) (Bonus: Schrittweise Modellwahl basierend auf AIC mit "step")

# Die backward-Variablenselektion bzgl. des AIC ist in der Funktion step
# implementiert. Hierzu ist das glm-Objekt des vollen Modells (object) zu
# übergeben, sowie "backward" für das Argument method zu wählen.

# Funktionsweise: Sukzessive wird das Modell jeweils mit Weglassen einer
# Variable berechnet. Jene Variable, deren Weglassen das AIC am stärksten
# verringert, wird im nächsten Schritt entfernt. Abbruch, falls keine
# Verringerung des AIC mehr auftritt.

gammaNatsel <- step(gammaNat,method="backward")
summary(gammaNatsel)

gammaLogsel <- step(gammaLog,method="backward")
summary(gammaLogsel)

# Beim Modell mit natürlichem Link sowie beim log-Link wird das Modell ohne 
# deutsch, alter und verh gewählt.
# Nach wie vor ist das Modell mit Log-Link zu bevorzugen. Wie zuvor ist hier 
# die Devianz bzw. das AIC niedriger im Vergleich zum Modell mit natürlichem 
# Link.




# -------- Aufgabe 2 -----------------------------------------------------------
# (b)

# 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()