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


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

# die Daten
library(MASS)
data(epil)
#help(epil)

# (a)

# GLM mit Poissonverteiltem y
glm1 <- glm(y ~ lbase + trt + lage + V4 + lbase:trt, family=poisson, data=epil)
summary(glm1)

# Nach den geschätzten Koeffizienten zu urteilen reduziert progabide die
# erwartete Anzahl an Anfällen, allerdings nur für weinger schwerwiegende Fälle
# (niedriges oder mäßiges lbase). Für hohes lbase ist der Zusammenhang dem 
# Anschein nach umgekehrt (positiver Koeffizient beim Wechselwirkungsterm!).
# Bei älteren Patienten ist die Anzahl der Anfälle offenbar höher, abenso
# (wie zu erwarten) bei einer großen Anzahl in der Baseline-Periode. In der 
# vierten Woche scheinen die Anfälle weniger zu werden.

# Problematisch an dem GLM ist jedoch: Die einzelnen Beobachtungen sind
# nicht unabhängig, da jeder Patient viermal eingeht. 
table(epil$subject,epil$period)
# Die Messungen an einer Person sind mit Sicherheit abhängig/positiv korreliert. 
# Daher sind insbesondere die geschätzten Standardfehler und somit die angegebenen 
# Teststatistiken nicht zuverlässig. 


# (b)

# Quasi-Poisson-Modell
glm2 <- glm(y ~ lbase + trt + lage + V4 + lbase:trt, family=quasipoisson, data=epil)
summary(glm2)

# Da nur ein Dispersionsparameter eingeführt wurde, bleiben die 
# Koeffizientenschätzer gleich. Der geschätzte Dispersionsparameter ist 
# allerdings (mit 4.4) deutlich größer als 1. Es liegt also Overdisperion vor
# (vermutlich auf Grund der positiven Korrelation zwischen den Beobachtungen).
# Daher sind die geschätzten Standardfehler nun deutlich größer, was 
# Auswirkungen auf die angegebene Teststatistiken hat (z.B. ist der Effekt von
# V4 nun nicht mehr signifikant).


# (c) ...


# (d)

# GEEs
library(gee)
#help(gee)

# mit working correlation "independence"
gee1 <- gee(y ~ lbase + trt + lage + V4 + lbase:trt, id=subject, family=poisson, 
corstr="independence", data=epil)
summary(gee1)

# mit equi-correlation, d.h. "exchangeable"
gee2 <- gee(y ~ lbase + trt + lage + V4 + lbase:trt, id=subject, family=poisson, 
corstr="exchangeable", data=epil)
summary(gee2)

# zusätzlich autoregressiv (erster Ordnung), d.h. mit "AR-M" (Mv=1)
gee3 <- gee(y ~ lbase + trt + lage + V4 + lbase:trt, id=subject, family=poisson, 
corstr="AR-M", Mv=1, data=epil)
summary(gee3)

# Bei "independence" ergeben sich für die Koeffizienten die gleichen Schätzwerte
# wie im GLM, bei "exchangable" zeigen sich hier (Spezialfall!) nur numerische
# Unterschiede. Auch bei "AR-1" sind sie qualtitativ ähnlich.

# Quantifizierung des Einflusses z.B. über
exp(gee3$coef["V4"])
# d.h. in der vierten Woche nimmt die Anzahl der Anfälle marginal um den 
# Faktor 0.86 ab.

# Die Spalte "Naive z" gibt den Standardfehler aus unter der Annahme, dass die 
# verwendete (working) covariance der wahren entspricht (mit Dispersion).
# "Robust z" verwendet die Schätzung über die Sandwich-Matrix.