################################################################################
#
# Beispiel: Risiko für CBR (Chronische Bronchitis) anhand Staubbelastung
#
# Vorlesung Lineare Modelle (Helmut Küchenhoff) 1.7 - 2.7.2014
# Thema: Das logistische Regressionsmodell 
#
# data: Staubdaten von 1246 Patienten mit Exposition und Krankheit cbr
#       file: dust.asc (http://www.statistik.lmu.de/~helmut/limo09/dust.asc) 
#
# contents: (1) Einlesen der Daten
#  	        (2) Explorative Analyse
#           (3) logistische Regression
#
# last.modified: 6.7.2015
#
################################################################################

#### (1) Daten Einlesen ########################################################

staub <- read.table("dust.asc", 
  header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
#staub <- read.table("dust.asc", header=TRUE, sep="", na.strings="NA", dec=".", 
#  strip.white=TRUE)

attach(staub)
head(staub)  ## die ersten werte ##
table(cbr)
table(smoking)
mean(cbr)
mean(smoking)
plot(dust,cbr)
abline(lm(cbr~dust))

### Transformation ###
staub$ldust<- log(1+staub$dust)
attach(staub)
plot(ldust,cbr)
abline(lm(cbr~ldust))
plot(lm(cbr~ldust))
#### Explorative Analyse mit interaktiver grafik (iplots)
#library(iplots)
#ibar(cbr)
#ibar( smoking)
#iplot(dust,expo)
#ihist(ldust)
#ihist(expo)

#### Logistische regression
library(Epi)

st1<-glm(cbr ~ smoking+ expo+ ldust, family=binomial(logit),data=staub)
summary(st1)
confint(st1)
### odds ratios   ###
exp(confint(st1))
exp(st1$coefficients)
summary(st1)
anova(st1)
library(car)
Anova(st1)

library(effects)
plot(allEffects(st1))
### Mit Rauchen als Kategoriale Variable) 
staub$csmoking<-as.factor(staub$smoking)
st1<-glm(cbr ~ csmoking+ expo+ ldust, family=binomial(logit),data=staub)
ylim<-c(0,1)
plot(allEffects(st1),ylim=rep(c(-2,0),3))

###
st2<- glm(cbr~csmoking+expo+ldust + I(ldust^2),family=binomial,data=staub)
summary(st2) 

### Staubeffekt
x<-seq(0,5,by=0.1) 
y<- st2$coefficients[4]*x+ st2$coefficients[5]*x^2
plot(x,y,type="l")  

###Effekt am Anfang eher gering 
### Schwellenwertmodell 
staub$ldustg<- (staub$ldust-1.27)*(staub$ldust>1.27)   
st3<- glm(cbr~smoking+expo+ldustg, family=binomial, data=staub)   
summary(st3)
y1<-st3$coefficients[4]* (x-1.27)*(x>1.27) 
points(x,y1,type="l",col="red")

anova(st3)
library(car)
Anova(st3)
###  ROC Kurve #####
library(Epi)
ROC(form = cbr ~ smoking+expo+ldust, data=staub)

ROC(form = cbr ~ smoking, data=staub)