remove(list=ls())
#setwd("/Volumes/MichaelsHD/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt09/R")
setwd("/Volumes/Daten/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt09/R")
set.seed(23514)

#Daten einlesen
oring <- read.table("./oring.dat", header=T)
temp <- oring$temp
temp
def <- oring$defekt
def

#Modell mit logit-Funktion
defmod.logit <- glm(def ~temp,family=binomial(link=logit))
summary(defmod.logit)

#Modell mit probit-Funktion
defmod.probit <- glm(def ~temp,family=binomial(link=probit))
summary(defmod.probit)

#Vorhersagen
pred.logit <- predict(defmod.logit,data.frame(temp=c(31)),type="response")
pred.logit
pred.probit <- predict(defmod.probit,data.frame(temp=c(31)),type="response")
pred.probit

#Plot der Daten
plot(temp,def,xlim=c(30,85))
#Hinzufuegen der Regressionskurven
temp2 <- sort(c(temp,31))
ylogit <- c(pred.logit,defmod.logit$fit[order(temp)])
lines(temp2,ylogit)
yprobit <- c(pred.probit,defmod.probit$fit[order(temp)])
lines(temp2,yprobit,col="2")
#hinzufuegen eines Punktes bei 31F
points(31,1,col="3")
#Hinzufuegen einer Legende
legend(x=30,y=0.2,legend=c("logit","probit"),col=c(1,2),lty=c(1,1))


#andere Moeglichkeit (glattere Kurven):
plot(temp,def,xlim=c(30,85))
#logit-Linkfunktion aus dem boot-Paket
library(boot)
x <- seq(30,85,by=0.5)
logitcoef <- defmod.logit$coef
ylogit2 <- inv.logit(logitcoef[1] + logitcoef[2]*x )
lines(x,ylogit2)
probitcoef <- defmod.probit$coef
#Erinnerung: Probitfunktion ist N(0,1)-Verteilungsfunktion
yprobit2 <- pnorm(probitcoef[1] + probitcoef[2]*x)
lines(x,yprobit2,col="2")
#hinzufuegen eines Punktes bei 31F
points(31,1,col="3")
#Hinzufuegen einer Legende
legend(x=30,y=0.2,legend=c("logit","probit"),col=c(1,2),lty=c(1,1))