## logistic regression

pfad = "das Verzeichnis der gespeicherten Dateien"
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
neuerwert = data.frame(jahr = 12)
ergebnis = predict(g, neuerwert, se.fit=T)
ergebnis$fit
abline(h=ergebnis$fit, lty=2, col="blue")
abline(v=12, col="blue", lty=2)

####
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) 
m = coef(g)[2]
k = coef(g)[1]

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)


## Logistic Regression und kategoriale Grenzen
# Abbildung der Proportionen als Funktion von F2
p = ui[,1]/apply(ui, 1, sum)

f2werte = as.numeric(rownames(ui))
plot(f2werte, p, ylab="Proportion /u/ Urteile", xlab="F2 (Hz)")

logui = glm(ui ~ f2werte, family=binomial)
m = coef(logui)[2]
k = coef(logui)[1]

curve(exp(m*x + k)/(1+ exp(m*x+k)), xlim=c(1200, 2400), add=T, col=2)
# 50% Umkipp-punkt ueberlagern
-k/m
abline(v=-k/m)


####################
## Signifikanz-Test

summary(g)
anova(g, test="Chisq")


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

lost2 = as.matrix(read.table(paste(pfad, "lost2.txt", sep="/")))
J = c(jahr, jahr)
G = c(rep(0,6), rep(1, 6)) 

### 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)

###########################################
## Logistic Regression mit nur 2 Ebenen

gmf = glm(lost2 ~ G, "binomial")
anova(gmf, test="Chisq")

## dasselbe
m = apply(lost2[1:6,], 2, sum)
f = apply(lost2[7:12,], 2, sum)
mf = rbind(m, f)
rownames(mf) = c(0, 1)
colnames(mf) = c("high", "low")
mf

l.mf = c(0,1)
gmf2 = glm(mf ~ l.mf, "binomial")
anova(gmf2, test="Chisq")


################# und ein aehnliches Ergebnis
## mit einem direkten chi-quadrat

chisq.test(mf)

####################################
### Drei Gruppen, jeweils 2 Ebenen

A = c(0, 0, 1, 1)
G = c(0, 1, 0, 1)
prop = lost3[,1]/apply(lost3, 1, sum)
interaction.plot(A, G, prop)

g = glm(lost3 ~ A * G, binomial)
anova(g, test="Chisq")