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


# -------- Aufgabe 7 ----------------------------------------------------------

# Einlesen der Daten
leafblotch <- read.table("leafblotch.dat", header = TRUE)
# Struktur der Daten
str(leafblotch)

# Die Kovariablen site und variety werden nun als Faktor-Variablen definiert
# (d.h. bei der Modellierung identifiziert sie R als Dummy-Variablen mit der
# ersten (!) Kategorie als Referenz)
leafblotch$sitef <- as.factor(leafblotch$site)
leafblotch$varietyf <- as.factor(leafblotch$variety)

# a)

# Das verlangte Modell kann man in R auf zweierlei Arten rechnen, mit
# family=quasibinomial (hier wird automatisch von der Varianzfunktion
# v(\mu)=\mu(1-\mu), Modellierung der Dispersion) ausgegangen:
quasiglm1a <- glm(blotch ~ sitef + varietyf, data = leafblotch, 
	family = quasibinomial(link = "logit"))
summary(quasiglm1a)

# oder mit family = quasi. Diese Funktion ist allgemeiner (funktioniert für
# verschiedene Link-Funktionen, z.B. auch "log", vgl. Quasi-Poisson-Modelle),
# man muss die Varianzfunktion explizit spezifizieren (es sind verschiedene
# weitere Möglichkeiten implementiert, z.B. mu^2 oder mu^3):
quasiglm1b <- glm(blotch ~ sitef + varietyf, data = leafblotch,
	family = quasi(link = "logit", var = "mu(1-mu)"))
summary(quasiglm1b)

# Man erkennt: Beide Verfahren liefern denselben Output. Der geschätzte
# Dispersionsparameter \phi ist mit 0.0887 weit unter dem Wert 1 (wie er bei
# Annahme einer Binomialverteilung festgelegt wäre). Bei Betrachtung der
# geschätzten Parameter erweisen sich die Sorten 1 und 2 als am stärksten
# resistent gegen die Blattkrankheit, während die Sorten 8-10 am stärksten
# befallen sind (in dieser Reihenfolge). Ähnliches gilt für die Felder.

# Beachte, dass man beim Schätzen unter Annahme einer Binomialverteilung die
# selben Schätzer erhält (aber andere Werte für Tests etc.); plus Warnung,
# weil Werte zwischen 0 und 1 bei Binomialverteilung ja nicht vorkommen dürften.
binomialModell  <- glm(blotch ~ sitef + varietyf, data = leafblotch, 
	family = binomial)
summary(binomialModell)
all.equal(coef(binomialModell), coef(quasiglm1a))


# b)

# Ein geeigneter Diagnostik-Plot (vgl. auch später) wäre, die Pearson-Residuen
# r_i^P = r_i/v(\hat\mu_i) gegen den geschätzten linearen Prädiktor (im GLM oft
# besser als gefittete Werte; auch wg. Skalierung, z.B. Poisson- oder
# Binomial-Modell, hier aber weniger relevant)

# gegen linearen Prädiktor
plot(predict(quasiglm1b), residuals(quasiglm1b, type = "pearson"), 
xlab = "linearer Prädiktor", ylab = "Pearson-Residuen")

# gegen gefittete Werte
plot(predict(quasiglm1b, type = "response"), residuals(quasiglm1b, type = "pearson"), 
xlab = "gefittete Werte", ylab = "Pearson-Residuen")

# Im Idealfall sollte man (wie im linearen Modell) keine Struktur erkennen
# können. Ferner sollte der Plot eine konstante Varianz aufweisen. Das ist hier
# nicht der Fall (die Varianz ist für kleine Werte des linearen Prädiktors
# geringer; bis ca. -5). Das ist ein Hinweis, dass die Varianzfunktion nicht
# gut gewählt ist.



# c)

# D(\beta) als Funktion
Dbeta <- function(eta) {
  1/(exp(-eta) + 2 + exp(eta))
}

# Design-Matrix X zusammenbauen
# leafblotch
X <- cbind(matrix(0, 10, 8), rbind(0, diag(1, 9)))
for (si in 1:8) {
    A <- matrix(0, 10, 8)
    A[, si] <- 1
    X <- rbind(X, cbind(A, rbind(0, diag(1, 9))))
}
X <- cbind(1, X)

# V-Marix (siehe Skript)
V <- t(X) %*% diag(Dbeta(predict(quasiglm1b))^2 * 
	residuals(quasiglm1b,type = "pearson")^2 / 
	(summary(quasiglm1b)$dispersion^2 * quasiglm1b$fitted * 
		(1-quasiglm1b$fitted))) %*% X

# geschätzte Kovarianz-Matrix F^{-1}VF^{-1}
Covbeta <- summary(quasiglm1b)$cov.scaled %*% V %*% 
	summary(quasiglm1b)$cov.scaled

# Varianzen (Beachte aber: wegen der vergleichsweise geringen Zahl an
# Beobachtungen ist diese Schätzungen vermutlich nicht besonders gut)
diag(Covbeta)
# bzw. Standardfehler
sqrt(diag(Covbeta))

# Vergleich mit Standardfehlern falls Var(y)  =  \phi\mu(1-\mu) richtig spezifiziert
sqrt(diag(summary(quasiglm1b)$cov.scaled))

# berechenbar mit 
W <- diag(Dbeta(predict(quasiglm1b))^2 / 
          (quasiglm1b$fitted *(1-quasiglm1b$fitted))/summary(quasiglm1b)$dispersion)
F2 <- t(X) %*% W %*% X
sqrt(diag(solve(F2)))


# (d)

#install.packages("gnm")
library(gnm)

# Modell fitten
glmwedderburn <- glm(blotch ~ sitef + varietyf, data = leafblotch, 
	family = wedderburn)

# Output
summary(glmwedderburn)

# Vergleich
summary(quasiglm1b)

# Man erkennt: die Parameterschätzer verändern sich. Man beachte, dass sich die
# Reihenfolge bei den Parametern für die Sorten nun etwas ändert, auch wenn die
# grobe Tendenz die selbe ist wie zuvor. Der Skalenparameter nimmt
# einen Wert von 0.988 an.

# Der Plot der Pearson-Residuen gegen den geschätzten linear Prädiktor,
 plot(predict(glmwedderburn), residuals(glmwedderburn, type = "pearson"), 
xlab = "linearer Prädiktor", ylab = "Pearson-Residuen")

# weist deutlich stärker in Richtung von gleichen Varianzen als zuvor, was
# darauf hindeutet, dass es sich hierbei um das bessere Modell handelt.