###############################################################################
## Inhalt:
## Fortsetzung von Kapitel 7.1 (Regressionsmodelle) des Kurses von 
## Ruckdeschel & Kohl
###############################################################################


###############################################################################
## Ausblick: Mixed-Effect Modelle
## Lesen Sie "LineareMixedEffectModels.pdf"
###############################################################################

###########################################################
## Lineare (und nicht lineare) mixed-effect Modelle sind in den Paketen "nlme"
## und "lme4" implementiert.
## Referenz: Pinheiro, J. und Bates D. (2000). Mixed-Effects Models in S and 
##           S-Plus. Springer Verlag.
###########################################################
library(nlme)


###########################################################
## Wir betrachten den Datensatz "Orthodont" aus dem Paket "nlme"
data(Orthodont)
?Orthodont

## Es handelt sich dabei um ein Objekt der S3-Klasse "nfnGroupedData". Diese
## Klasse ist abgeleitet von "nfGroupedData" -> "groupedData" -> "data.frame".
class(Orthodont)

## Es gibt für diese Klassen spezielle plot-Methoden
plot(Orthodont)

## oder auch
plot(Orthodont, outer = TRUE)

## Für weitere Einzelheiten siehe
?plot.nfnGroupedData


#######################################
## Wir wählen die weiblichen Jugendlichen aus.
OrthoFem <- Orthodont[Orthodont$Sex == "Female",]

## Wir plotten nur die Daten der weiblichen Jugendlichen
plot(OrthoFem, outer = TRUE)

## Wir untersuchen die gemessene Distanz in Abhängigkeit vom Alter.

## Im ersten Schritt passen wir für jede der weiblichen Jugendlichen ein einfaches
## lineares Modell an. Hierfür können wir die Funktion "lmList" aus dem Paket
## "nlme" verwenden. Zusätzlich verschieben wir die Daten um 8; d.h., der Intercept
## spiegelt das Alter zum Beginn der Messungen wider.
lm1 <- lmList(distance ~ I(age-8), data = OrthoFem)
summary(lm1)

## Wir stellen die geschätzten Parameter inklusive entsprechender 95%
## Konfidenzintervalle graphisch dar
plot(intervals(lm1))

## Wir sehen das das Wachstum, d.h., die Steigung der Geraden bei allen weiblichen
## Jugendlichen sehr ähnlich ist. Jedoch weichen die Startpunkte voneinander ab.

## Dies modelieren wir jetzt durch einen sog. "random effect".
## Im zweiten Schritt verwenden wir nun die Daten aller weiblichen Jugendlichen 
## simultan und verwenden ein lineares mixed-effect Modell, wobei wir als 
## "random effect" das entsprechende Subjekt wählen.
lme1 <- lme(distance ~ I(age-8), data = OrthoFem, random = ~1)
summary(lme1)

## Wir erhalten also einen Intercept von ca. 21.20 und zusätzlich einen random
## Effekt mit Verteilung N(0, 2.068^2). Dieses Modell erklärt sehr gut die
## unterschiedlichen Werte zum Beginn der Messungen im Alter von 8 Jahren.


#######################################
## Allgemein: 
## Ein "fixed effect" ist eine Größe/Variable, die bei der Wiederholung des 
## Experiments unabhängig von der konkreten Stichprobe immer wieder so gemessen 
## würde. Hier das Wachstum.
##
## Ein "random effect" ist eine Größe/Variable, die von der Zufallsauswahl abhängt.
## Hier: Jedes Mädchen weist im Alter von 8 Jahren eine individuelle Distanz vor,
## die nicht für alle Mädchen einheitlich gleich groß ist.
#######################################


## Wir wiederholen die Analyse mit dem Paket "lme4".
library(lme4)
lme2 <- lmer(distance ~ I(age-8) + (1|Subject), data = OrthoFem)
summary(lme2)

## bzw. soll etwas schneller sein
lme3 <- lmer2(distance ~ I(age-8) + (1|Subject), data = OrthoFem)
summary(lme3)


#######################################
## Die Ausgabe bei "lmer" bzw. "lmer2" ist offenbar leicht anders als bei "lme"
## beinhaltet aber dieselbe Information. 
## Die Funktion "lmer" bzw. "lmer2" ist allgemeiner als "lme" und erlaubt auch
## die Betrachtung sog. generalisiert linearer mixed-effect Modelle.
## Was man unter generalisiert linearen Modellen versteht, werden wir im Folgenden 
## kennenlernen.
#######################################


###############################################################################
## Generalisiert Lineare Modelle
## Lesen Sie bitte Abschnitt 7.1.2 im R-Kurs und
## http://en.wikipedia.org/wiki/Generalized_linear_model
###############################################################################

###########################################################
## Binomiale Daten 
## Daten-Quelle: p. 353 in Altman (1991). Practical Statistics for Medical
## Research. Chapman & Hall.
## vgl. auch Abschnitt 11.2 in P. Dalgaard (2002). Introductory Statistics 
## with R. Springer Verlag.
###########################################################
no.yes <- c("No", "Yes")
smoking <- gl(2, 1, 8, no.yes)

## Die Funktion "gl" (= generate levels) erzeugt einen Vektor vom Typ Faktor 
## durch einen Aufruf von "rep" und "factor"
gl

## Wir erhalten
smoking

## Wir erzeugen die weiteren Variablen
obesity <- gl(2, 2, 8, no.yes) # Fettleibigkeit
snoring <- gl(2, 4, 8, no.yes)
n.tot <- c(60, 17, 8, 2, 187, 85, 51, 23)
n.hyp <- c(5, 2, 1, 0, 35, 13, 15, 8) # Bluthochdruck
n.no.hyp <- n.tot - n.hyp

## Wir erhalten also den folgenden Datensatz
data.frame(smoking, obesity, snoring, n.hyp, n.no.hyp)

## Wir wollen untersuchen, ob die Variablen "smoking", "obesity" und "snoring"
## Risikofaktoren für Bluthochdruck sind.

## Als Linkfunktion verwenden wir die kanonische Linkfunktion "logit".
## Um das logistische Regressionsmodell zu berechnen, verwenden wir die 
## Funktion "glm". Das entsprechende Modell läßt sich auf zweierlei Arten 
## formulieren.

#######################################
## 1. Möglichkeit: Matrix mit Spalten "Erfolge / Misserfolge"
hyp.mat <- cbind(n.hyp, n.no.hyp)
fit1 <- glm(hyp.mat ~ smoking + obesity + snoring, family = binomial("logit"))

## Wir erhalten
fit1

## Bzw. ausführlicher
summary(fit1)

#######################################
## 2. Möglichkeit: verwende relative Häufigkeiten ("Erfolgswahrscheinlichkeiten")
## Erfordert zusätzlich den weights-Parameter, damit klar ist auf der Basis
## wievieler Beobachtungen die "Erfolgswahrscheinlichkeiten" berechnet wurden.
hyp.rel <- n.hyp/n.tot
fit2 <- glm(hyp.rel ~ smoking + obesity + snoring, family = binomial("logit"),
            weights = n.tot)
summary(fit2)

## Die Variable "smoking" ist also nicht signifikant.
## Wir sehen uns zusätzlich noch die Korrelation zwischen den Schätzungen an.
summary(fit2, corr = TRUE)

## Die Korrelation zwischen den Variablen ist sehr klein und wir können daher
## erwarten, dass die Signifikant für "obesity" und "snoring" erhalten bleibt,
## auch wenn wir die Variable "smoking" entfernen.
fit3 <- glm(hyp.mat ~ obesity + snoring, family = binomial("logit"))
summary(fit3)


#######################################
## Eine andere Möglichkeit der Modellwahl ist die Verwendung der Devianz, welche
## asymptotisch Chi^2 verteilt ist.
anova(fit1, test = "Chisq")

## Bei dieser Analyse ist aber nun entscheidend, in welcher Reihenfolge die
## Variablen in der Formel auftreten. Es wird also sukzessive getestet.
## Nullmodell vs. Nullmodell + smoking
## Nullmodell + smoking vs. Nullmodell + smoking + obesity
## Nullmodell + smoking + obesity vs. Nullmodell + smoking + obesity + snoring

## Die verbleibenden Möglichkeiten sind:
anova(glm(hyp.mat ~ smoking + snoring + obesity, family = binomial("logit")), test = "Chisq")
anova(glm(hyp.mat ~ obesity + smoking + snoring, family = binomial("logit")), test = "Chisq")
anova(glm(hyp.mat ~ snoring + smoking + obesity, family = binomial("logit")), test = "Chisq")
anova(glm(hyp.mat ~ obesity + snoring + smoking, family = binomial("logit")), test = "Chisq")
anova(glm(hyp.mat ~ snoring + obesity + smoking, family = binomial("logit")), test = "Chisq")

## Wir sind besonders an den letzten beiden Devianz-Analysen interessiert, da 
## diese zeigen, dass die Variable "smoking" keinen signifikanten Beitrag leistet.


#######################################
## Ein weitere Möglichkeit der Modellwahl besteht in der Verwendung des 
## AIC-Kriteriums. Wir verwenden die Funktion "stepAIC" aus dem Paket "MASS".
library(MASS)
fit4 <- stepAIC(fit1)
summary(fit4)

## Erneut erhalten wir das Modell "hyp.mat ~ obesity + snoring".


#######################################
## Was nun die Mediziner interessiert ist die odds-Ratio.
## Konkret also zum Beispiel: Das Chancenverhältnis für Patienten mit und ohne
## Fettleibigkeit an Bluthochdruck zu leiden.
## Die Berechnung übernimmt die Funktion "logistic.display" aus dem Paket "epicalc".
library(epicalc)
logistic.display(fit3)

## Die Personen, die unter Fettleibigkeit leiden, haben also ein (geschätztes) 
## 2-mal höheres Risiko (signifikant verschieden zu 1) an hohen Blutdruck zu leiden 
## wie Personen ohne Fettleibigkeit.
## Ähnlich haben Personen, die Schnarchen ein (geschätztes) 2.3-mal höheres Risiko
## (signifikant verschieden zu 1).

## Zum Vergleich
logistic.display(fit1)

## Die odds-Ratio für den Faktor "smoking" ist also nicht signifikant verschieden
## von 1; d.h., Rauchen stellt im Fall der vorliegenden Daten keinen Risikofaktor 
## für Bluthochdruck dar.


#######################################
## Im Fall solcher tabelarischer Daten ist es naheliegend, das Modell dadurch
## zu überprüfen, dass man die gefitteten Werte mit den tatsächlichen Werten
## vergleicht.
data.frame(fit = fitted(fit3)*n.tot, n.hyp, n.tot)


###########################################################
## Binomiale Daten
###########################################################
## Die andere Möglichkeit ist, dass die binomialen Daten in Rohform, d.h., 
## nicht zusammengefasst vorliegen. Ein Beispiel ist der folgende Datensatz,
## mit dessen Hilfe Riskofaktoren für ein niedriges Geburtsgewicht bestimmt
## werden sollten; vgl. p. 194 in Venables, W. N. and Ripley, B. D. (2002). 
## Modern Applied Statistics with S. Fourth edition. Springer.
data(birthwt)
?birthwt

str(birthwt)
## Wir sehen, dass die Variablen einheitlich als "integer" Variablen kodiert
## sind. Wir ändern die Kodierung (vgl. example(birthwt))
attach(birthwt)
race <- factor(race, labels = c("white", "black", "other"))
ptd <- factor(ptl > 0)
ftv <- factor(ftv)
levels(ftv)[-(1:2)] <- "2+"
bwt <- data.frame(low = factor(low), age, lwt, race, smoke = (smoke > 0), ptd, 
                  ht = (ht > 0), ui = (ui > 0), ftv)
detach("birthwt")

## Wir erhalten
bwt

## Wir passen nun ein generalisiert lineares Modell an die Daten an, wobei wir
## dieses Mal mit dem "probit" Link arbeiten.
fit1 <- glm(low ~ ., family = binomial("probit"), data = bwt)
summary(fit1)

## Wir verwenden die Funktion stepAIC, um ein geeignetes Teilmodell auszuwählen.
fit2 <- stepAIC(fit1)
summary(fit2)

#######################################
## Für eine weitere tiefergehende Analyse dieser Daten siehe p. 194ff in
## Venables and Ripley (2002).
#######################################


###########################################################
## Multinomiale Daten.
## vgl. Venables and Ripley (2002), p. 203.
###########################################################
data(housing)
?housing

## Wir verwenden die Funktion "multinom" aus dem Paket "nnet", um das entsprechende
## multinomial Modell anzupassen.
library(nnet)

#######################################
## Wir untersuchen die Zufriedenheit der Hausbesitzer mit den aktuellen 
## Wohnbedingungen.
fit1 <- multinom(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
fit1

## etwas umfangreichere Ausgabe
summary(fit1)


#######################################
## Wir lassen nun zusätzlich Interaktionen zwischen den Variablen zu.
fit2 <- multinom(Sat ~ Infl*Type*Cont, weights = Freq, data = housing)
summary(fit2)


## Wir vergleichen die beiden Fits mittels einer Devianzanalyse.
anova(fit1, fit2)

## Das komplizierte führt also zu keiner signifikanten Reduktion der Devianz.
## Somit ist das einfachere Modell vorzuziehen.


#######################################
## Ein alternative Möglichkeit für diesen Datensatz bestünde darin ein 
## Poisson-Modell zu verwenden; vgl. p. 199 ff in Venables and Ripley (2002).
#######################################


###########################################################
## Poisson-verteilte Daten
## Der folgende Datensatz behandelt die Fehltage von High School Schülern.
## Es handelt sich dabei um eine STATA-Datensatz, den wir aber ohne Probleme
## mittels der Funktion "read.dta" aus dem Paket "foreign" einlesen können.
###########################################################
library(foreign)
HS <- read.dta(file = "poissonreg.dta")

## Wir betrachten die Abwesenheitstage "daysabs" und untersuchen den Einfluss
## von Geschlecht ("male") sowie den Ergebnissen im Mathematik- ("math") und 
## Sprachentest ("langarts").

## Wir wandeln die Variable "male" in eine logische Variable um.
HS$male <- as.logical(HS$male)

## Wir passen nun das entsprechende Poisson Regressionsmodell an.
fit1 <- glm(daysabs ~ male + math + langarts, family = poisson, data = HS)
summary(fit1)

## Wir betrachten die Devianz.
anova(fit1)

## Wie verhält es sich mit dem AIC-Kriterium
stepAIC(fit1)