## logistic regression
## Bitte lost.txt, lost2.txt, lost3.txt, ui.txt
## aus der Webseite in einem Verzeichnis speichern


pfad = "/Volumes/Data_1/d/"
lost = as.matrix(read.table(paste(pfad, "lost.txt", sep="/")))
lost2 = as.matrix(read.table(paste(pfad, "lost2.txt", sep="/")))
lost3 = as.matrix(read.table(paste(pfad, "lost3.txt", sep="/")))
ui = as.matrix(read.table(paste(pfad, "ui.txt", sep="/")))

jahr = as.numeric(rownames(lost))
jahr = jahr - 1950

g = glm(lost ~ jahr, binomial)


p = lost[,1]/apply(lost, 1, sum)
# log-odds
lodd = log(p/(1-p))
plot(jahr, lodd, type="b")
# Regressionslinie ueberlagern
abline(g, col=2)
text(jahr, predict(g), "x", col=3)


# Vorhersage
m = coef(g)[2]
k = coef(g)[1]
abline(h=m * 12 + k, lty=2, col="blue")
abline(v= 12)

####
p = 0.8
L = log(p/(1-p))
# in die andere Richtung
exp(L)/(1+exp(L))

###### Abbildung Jahr x Proportionen
p = lost[,1]/apply(lost, 1, sum)
plot(jahr,p) 

curve(exp(m*x + k)/(1+ exp(m*x+k)), xlim=c(0, 60), add=T, col=2)
vorher = predict(g)
text(jahr, exp(vorher)/(1+exp(vorher)), "x", col=3)
# Vorhersage fŸr 1962
x = 12
exp(m*x + k)/(1+ exp(m*x+k))


###########################
### Zwei unabhaengige Variablen

lost = as.matrix(read.table(paste(pfad, "lost.txt", sep="/")))
lost2 = as.matrix(read.table(paste(pfad, "lost2.txt", sep="/")))
r = rownames(lost2)
J = as.numeric(substring(r, 1, 2))
G = factor(substring(r, 4, 4))

### Interaction-plot
p = lost2[,1]/apply(lost2, 1, sum)
interaction.plot(J, G, p)


# gibt eine Warnung bez. G - ignorieren
mehrg = glm(lost2 ~ J + G, binomial)

# pruefen ob wir dasselbe Ergebnis bekommen
# wie anova(glm(lost ~ jahr, binomial), test="Chisq")
#
g2 = glm(lost2 ~ J, binomial) 
anova(g2,test="Chisq")
 
## zurueck zur mehrfachen Regression
mehrg
anova(mehrg, test="Chisq")


########################################
### Interaktion


g = glm(lost2 ~ J + G + J:G, binomial)
# dasselbe: g = glm(lost2 ~ J * G, binomial)
anova(g, test="Chisq")


library(MASS)
stepAIC(g)