# @Author: andreas.bender@stat.uni-muenchen.de
# @Date:   2015-07-09 13:05:29
# @Last Modified by:   andreas.bender@stat.uni-muenchen.de
# @Last Modified time: 2015-07-09 13:09:09

library(ggplot2)
theme_set(theme_bw())
library(reshape2)


## -------------------------------- Aufgabe 1   ----------------------------  ##
bauteile <- read.csv("../../daten/bauteile.csv", sep=";", dec=".", header=TRUE)
bauteile$observed <- with(bauteile, y/n)

# Anlegen neuer Variablen, die Einzelbeobachtungen beinhalten
belastung <- c(
  rep(bauteile$belastung, times=bauteile$y), 
  rep(bauteile$belastung, times=(bauteile$n - bauteile$y)))
versagt <- c(
  rep(1, times=sum(bauteile$y)), 
  rep(0, times=sum(bauteile$n - bauteile$y)))

bauteile.long <- data.frame(belastung, versagt)

# Ueberpruefung
# table(versagt, belastung)
mosaicplot(~  belastung + versagt, data=bauteile.long, col=c("grey90", "grey30"), 
  las = 1)

# Logit-Modell
mod.bauteile<- glm(versagt ~ belastung, data=bauteile.long, family=binomial())
summary(mod.bauteile)


# graphische Darstellung
# geschaetzte Wahrscheinlichkeiten
bauteile$predicted <- predict(mod.bauteile, bauteile, type="response")
# datensatz in long format
mbauteile <- melt(bauteile[, c("belastung", "observed", "predicted")], id.vars="belastung")
ggplot(mbauteile, aes(x=belastung, y=value)) + 
  geom_point(aes(col=variable)) + 
  geom_line(aes(col=variable)) + 
  scale_color_discrete("Beobachete und prädiktierte\n Versagenswahrscheinlichkeit") + 
  ylab("Versagenswahrscheinlichkeit") + 
  xlab("Belastung")

## ----versagenbauteileBeispiele, echo=TRUE--------------------------------
# Veranschaulichung der geschaetzten Funktion
bauteile.long.new <- data.frame(belastung = seq(min(bauteile$belastung), 
  max(bauteile$belastung), length=1000))
bauteile.long.new$predicted <- predict(mod.bauteile, newdata=bauteile.long.new, 
  type="response")

ggplot(bauteile.long.new, aes(x=belastung, y = predicted)) + 
  geom_line() + 
  # x_50
  geom_hline(yintercept=0.5, lty=3) + 
  geom_vline(xintercept=-coef(mod.bauteile)[1]/coef(mod.bauteile)[2], lty=3) + 
# Wahrscheinlichkeit des Versagens bei Belastung von 5000
  geom_vline(xintercept=5000, lty=1, col=2) + 
  geom_hline(yintercept=1/(1+exp(sum(-coef(mod.bauteile)*c(1, 5000)))), col=2)


## --------------------------- Aufgabe 2 ------------------------------------ ##
library(faraway)
data(teengamb)
teengamb$gambbin <- (teengamb$gamble > 3) * 1
teengamb$sex <- factor(teengamb$sex, levels=0:1, labels=c("male", "female"))
# anpassen eines logistischen Modells 
glm.gamble <- glm(gambbin ~ income + sex, data = teengamb, family = binomial())
summary(glm.gamble)