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

################################################################################
#
# -------- Aufgabe 26 ----------------------------------------------------------
#
################################################################################

# a)

# Daten
library(geepack)
data(ohio)
ohio[1:30,]

glm1 <- glm(resp ~ age*smoke, data=ohio, family=binomial(link=probit))
summary(glm1)


# b)

library(lme4)
glmer1 <- glmer(resp ~ age*smoke + (1| id), family=binomial(link=probit),
data=ohio, nAGQ=20)
summary(glmer1)

sigmab2 <- VarCorr(glmer1)$id[1,1]
betavector <- fixef(glmer1)
betavector/sqrt(1 + sigmab2)

coef(glm1)

# Stimmt hier nur einigermaßen; Grund: Approximation des Integrals
# bei glmer: Gauß-Hermite-Quadratur (adaptive Wahl der Stützstellen)


# c)

glm2 <- glm(resp ~ age*smoke, data=ohio, family=binomial(link=logit))
summary(glm2)


# d)


# gemischte Modelle

library(glmmML)

?glmmML
# nur Random Intercept Modelle möglich

# ML mit Gauß-Hermite-Quadratur

ghq20 <- glmmML(resp ~ age*smoke, method="ghq", cluster=id, data=ohio, n.points=20)
summary(ghq20)

ghq15 <- glmmML(resp ~ age*smoke, method="ghq", cluster=id, data=ohio, n.points=15)
summary(ghq15)

ghq10 <- glmmML(resp ~ age*smoke, method="ghq", cluster=id, data=ohio, n.points=10)
summary(ghq10)

ghq5 <- glmmML(resp ~ age*smoke, method="ghq", cluster=id, data=ohio, n.points=5)
summary(ghq5)


# ML mit Laplace-Approximation ( = Gauß-Hermite-Quadratur mit n.points=1)

laplace1 <- glmmML(resp ~ age*smoke, method="Laplace", cluster=id, data=ohio)
summary(laplace1)

laplace2 <- glmmML(resp ~ age*smoke, method="ghq", cluster=id, data=ohio, n.points=1)
summary(laplace2)

# Ergebnisse werden von der Anzahl der Stützstellen beeinflusst

# e)

library(lme4)

# ML mit Gauß-Hermite-Quadratur (mit adaptiver Wahl der Stütztstellen)

glmer2 <- glmer(resp ~ age * smoke + (age | id), family="binomial", data=ohio, nAGQ=20)


# ML mit Laplace-Approximation ( = Gauß-Hermite-Quadratur mit nAGQ=1)

glmer3 <- glmer(resp ~ age * smoke + (age | id), family="binomial", data=ohio, nAGQ=1)
summary(glmer3)


# f)

library(MASS)

# mit Penalized Quasi-Likelihood  (als Laplace-Approximation motivierbar)
pql1 <- glmmPQL(resp ~ age * smoke , random = ~age | id, data=ohio,
 family=binomial(link=logit))
summary(pql1)