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



################################################################################
#
# -------- Aufgabe 21 ----------------------------------------------------------
#
################################################################################


# b) 

bike <- read.table("bike.dat", header=T)

# Poisson-Modell
model1a <- glm(y ~ lecture + offset(log(time)), family=poisson, data=bike)
summary(model1a)

# äquivalente Schreibweise
model1b <- glm(y ~ lecture, offset=log(time), family=poisson(link=log), 
               data=bike)
summary(model1b)

# Quasi-Poisson-Modell
model2 <- glm(y ~ lecture, offset=log(time), family=quasipoisson(link=log), 
              data=bike)
summary(model2)


################################################################################
#
# -------- Aufgabe 22 ----------------------------------------------------------
#
################################################################################

library(pscl)
data(bioChemists)
help(bioChemists)
attach(bioChemists)

# a)

# Verteilung der Anzahl an Artikeln
plot(as.numeric(names(table(art))),table(art),type="h", xlab="articles",
ylab="freq")

# nach Geschlecht
boxplot(art ~ fem, varwidth=T)

# nach Familienstand
boxplot(art ~ mar, varwidth=T)

# nach Kinderzahl
boxplot(art ~ kid5, varwidth=T)

# nach Prestige
plot(phd,art)
lines(lowess(phd,art))

# nach Betreuer
plot(ment,art)
lines(lowess(ment,art))

# c)

# Poisson-Modell
poissonModell <- glm(art ~ fem + mar + kid5 + phd + ment, family=poisson, 
data=bioChemists)
summary(poissonModell)

# Quasi-Poisson-Modell 
qpoissonModell <- glm(art ~ fem + mar + kid5 + phd + ment, family=quasipoisson, 
                     data=bioChemists)
summary(qpoissonModell)


# d)

# Negativbinomial-Modell
negbinomial <- glm.nb(art ~ fem + mar + kid5 + phd + ment, link=log, 
data=bioChemists)
summary(negbinomial)

# f)

help(zeroinfl)

# Poisson - Binomial(Logit-Link)
zeroinflPBl <- zeroinfl(art ~ fem + mar + kid5 + phd + ment, data=bioChemists, 
                        dist="poisson", link="logit")
summary(zeroinflPBl)

# Poisson - Binomial(Probit-Link)
zeroinflPBp <- zeroinfl(art ~ fem + mar + kid5 + phd + ment, data=bioChemists, 
                        dist="poisson", link="probit")
summary(zeroinflPBp)

# Negativbinomial - Binomial(Logit-Link)
zeroinflNBl <- zeroinfl(art ~  fem + mar + kid5 + phd + ment, data=bioChemists, 
                        dist="negbin", link="logit")
summary(zeroinflNBl)

# Negativbinomial - Binomial(Probit-Link)
zeroinflNBp <- zeroinfl(art ~  fem + mar + kid5 + phd + ment, data=bioChemists, 
                        dist="negbin", link="probit")
summary(zeroinflNBp)

zeroinflPBl2 <- zeroinfl(art ~ fem + mar + kid5 + phd + ment|fem + mar, 
                         data=bioChemists, dist="poisson", link="logit")
summary(zeroinflPBl2)


# h)

help(hurdle)

# Poisson - Poisson
hurdlePP <- hurdle(art ~ fem + mar + kid5 + phd + ment, data=bioChemists, 
dist="poisson", zero.dist="poisson")
summary(hurdlePP)
# Test, ob Koeffizientenvektoren unterhalb und oberhalb der Hürde gleich sind:
# H_0: \beta_1 = \beta_2 wird abgelehnt 
hurdletest(hurdlePP)

# Negativbinomial - Poisson
hurdleNP <- hurdle(art ~ fem + mar + kid5 + phd + ment, data=bioChemists, 
dist="negbin", zero.dist="poisson")
summary(hurdleNP)

# Poisson - Binomial
hurdlePB <- hurdle(art ~ fem + mar + kid5 + phd + ment, data=bioChemists, 
dist="poisson", zero.dist="binomial", link="logit")
summary(hurdlePB)

# Negativbinomial - Binomial
hurdleNB <- hurdle(art ~ fem + mar + kid5 + phd + ment, data=bioChemists, 
dist="negbin", zero.dist="binomial", link="logit")
summary(hurdleNB)

HurdlePP2 <- hurdle(art ~ fem + mar + kid5 + phd + ment|fem + mar, 
data=bioChemists, dist="poisson", zero.dist="poisson")
summary(HurdlePP2)


# i)

# hier am leichtesten zu bestimmen: AIC
AIC(poissonModell,negbinomial,hurdlePP,hurdleNP,hurdlePB,hurdleNB,
zeroinflPBl,zeroinflPBp,zeroinflNBp,zeroinflPBl2)