# ------------------------------------------------------------------------------
# -------- Übung zur Vorlesung Kategoriale Daten -------------------------------
# -------- R Code / Blatt 5                      -------------------------------
# ------------------------------------------------------------------------------



################################################################################
#
# -------- Aufgabe 19 ----------------------------------------------------------
#
################################################################################

library(aod)
data(orob2)
?orob2
orob2
attach(orob2)
response <- cbind(y,n-y)


# a)

# Binomial-Modell nur mit unterschiedlicher Gattung
Bin1 <- glm(response ~ seed, family=binomial)
summary(Bin1)

# Binomial-Modell mit unterschiedlicher Gattung und Wurzel-Extrakten
Bin2 <- glm(response ~ seed + root, family=binomial)
summary(Bin2)

# Binomial-Modell mit unterschiedlicher Gattung  und Wurzel-Extrakten 
# plus Wechselwirkung
Bin3 <- glm(response ~ seed*root, family=binomial)
summary(Bin3)

# Vergleich der Modelle mittels AIC:
AIC(Bin1, Bin2, Bin3)

# Vergleich der Modelle mittels Devianz-Analyse
anova(Bin1,Bin2,Bin3)
1-pchisq(Bin1$dev - Bin2$dev, 1)
1-pchisq(Bin2$dev - Bin3$dev, 1)
1-pchisq(Bin1$dev - Bin3$dev, 2)

# => Modell Bin3


# b)

# Zu beachten:
# - In dem Argument formula von betabin() muss die Zielgröße die Struktur 
#   cbind(y,n-y) haben (d.h. selbst wenn response=cbind(y,n-y) gilt, darf man 
#   nicht response schreiben) 
# - \phi im Output ist das \delta aus der Übung/Vorlesung

# Beta-Binomial-Modell nur mit unterschiedlicher Gattung
BetaBin1 <- betabin(cbind(y,n-y) ~ seed, random=~1, data=orob2)
summary(BetaBin1)

# Beta-Binomial-Modell mit unterschiedlicher Gattung und Wurzel-Extrakten
BetaBin2 <- betabin(cbind(y,n-y) ~ seed + root, random=~1, data=orob2)
summary(BetaBin2)

# Beta-Binomial-Modell mit unterschiedlicher Gattung  und Wurzel-Extrakten 
# plus Wechselwirkung
BetaBin3 <- betabin(cbind(y,n-y) ~ seed * root, random=~1, data=orob2)
summary(BetaBin3)

# Vergleich der Modelle mittels AIC:
AIC(BetaBin1, BetaBin2, BetaBin3)

# Vergleich der Modelle mittels Devianz-Analyse
anova(BetaBin1,BetaBin2,BetaBin3)
1-pchisq(BetaBin1@dev - BetaBin2@dev, 1)
1-pchisq(BetaBin2@dev - BetaBin3@dev, 1)
1-pchisq(BetaBin1@dev - BetaBin3@dev, 2)

# => Modell BetaBin3


# c)

# Quasi-Binomial-Modell nur mit unterschiedlicher Gattung
QuasiBin1 <- glm(response ~ seed, family=quasibinomial)
summary(QuasiBin1)

# Quasi-Binomial-Modell mit unterschiedlicher Gattung und Wurzel-Extrakten
QuasiBin2 <- glm(response ~ seed + root, family=quasibinomial)
summary(QuasiBin2)

# Quasi-Binomial-Modell mit unterschiedlicher Gattung  und Wurzel-Extrakten 
# plus Wechselwirkung
QuasiBin3 <- glm(response ~ seed*root, family=quasibinomial)
summary(QuasiBin3)

# Bei Quasi-Binomial-Modelle keine AIC- und Devianz-Bestimmung möglich, da keine
# Verteilung spezifiziert wird (Achtung! glm() gibt bei Quasi-Modell die 
# Devianzen des gewöhnlichen Binomial-Modells aus!!!)


# d) 

#-----------------------------------
# Vergleich der Outputs
#-----------------------------------

# Schätzer für Koeffizienten im Quasi-Binomial-Modell und im Binomial-Modell
# sind dieselben, da \phi in der Schätzgleichung rausfliegt

# Standardfehler im Quasi-Binomial-Modell größer als im Binomial-Modell
# -> siehe Zusammenhang mit \phi

# Standardfehler im Beta-Binomial-Modell größer als im Binomial-Modell
# -> Beachte: Nur Overdispersion-Modellierung im Beta-Binomial-Modell möglich


#-----------------------------------
# Möglichkeiten des Modellvergleichs
#-----------------------------------

# Vergleich der Binomial- und der Beta-Binomial-Modelle mit AIC möglich
# => Modell BetaBin3

# Vergleich der Binomial- und der Beta-Binomial-Modelle mit Devianzen wegen 
# hierarchischer Struktur möglich: 
# Binomial-Modell = Beta-Binomial-Modell bei \delta=0 (und damit \phi=1) 
1-pchisq(Bin1$dev - BetaBin1@dev, 1)
1-pchisq(Bin2$dev - BetaBin2@dev, 1)
1-pchisq(Bin3$dev - BetaBin3@dev, 1)
# => Modell Bin3

# Ein Vergleich mit den Quasi-Binomialmodellen über AIC bzw. Devianz ist nicht 
# möglich, da Vtlg nicht bekannt. 

# Vergleich sämtlicher Modelle über Kreuzvalidierung:
# Kriterien z.B.
#  - Trefferquote
#  - Kullback-Leibler-Distanz