#===========================================================================
# Zusatzmaterial zu Übungsblatt 6, Aufgabe 13 + 14
# Statistik 4 für Nebenfachstudierende SoSe 2017
#
# Datum: 22.06.2017
# Autor: Micha Schneider, Margret Oelker, Nora Fenske
#=============================================================================

# data
library("AER")
data("BankWages")
?BankWages

###############################################################################
# Multinomiales Modell
###############################################################################

# fit multinomial model for all employees
# nnet:
library("nnet")
m1 <- multinom(job ~ education + minority, data = BankWages)
summary(m1)

# VGAM
library("VGAM")
m2 <- vglm(factor(job, ordered=FALSE) ~ education + minority, data = BankWages,
           family=multinomial(refLevel=1))
summary(m2)

# Signifikanz-Test 
z_value <- coef(summary(m2))[4,1]/coef(summary(m2))[4,2]
z_value
# Vergleich mit kritischem Wert:
abs(z_value) > qnorm(0.975)
# Berechnung des p-Werts:
2*pnorm(-abs(z_value))

###############################################################################
# Kumulatives Modell
###############################################################################
# cumulative logit model for male employees 

# VGAM
library("VGAM")
m3 <- vglm(job ~ education + minority, family=cumulative(parallel=T), data = BankWages,
           subset = gender == "male")
summary(m3)

# MASS
library("MASS")
m4 <- polr(job ~ education + minority, data = BankWages,
           subset = gender == "male", Hess = "TRUE")
round(summary(m4)$coef,4)

all.equal(AIC(m3),AIC(m4)) # AICs gleich

# Achtung: Vorzeichen der Koeffizienten "vertauscht", 
#          Aber selbe inhaltliche Interpretation, da andere Strukturannahme (-)

# Achtung: m2 und m3 nicht mit AIC vergleichbar, da m3 nur ein Modell für die 
#          Männer fittet und m2 für Männer und Frauen (-> n sehr unterschiedlich)

###############################################################################
# Sequentielles Modell
###############################################################################
# Example: sequential logit model for male employees (inhaltlich nicht sinnvoll)
m5 <- vglm(job ~ education + minority, family=sratio(parallel=T), data = BankWages,
           subset = gender == "male")
summary(m5)