library(lme4)
library(multcomp)
###################################################
alt = read.table(url("http://www.phonetik.uni-muenchen.de/~jmh/lehre/sem/ss11/statfort/alter.txt"))

###################################################

tab = with(alt, table(Alter, Correct))
barchart(tab, auto.key=T)
###################################################

a = lmer(Correct ~ Alter + (1|Hoerer) + (1|Sprecher), family=binomial, data = alt)
fixef(a)
ranef(a)
contrasts(alt$Alter)
# Für Reihe 2
L = 0.7488927 + 1.4288235 + 0.13556251 -0.9330030
exp(L)/(1 + exp(L))
# das gleiche
fitted(a)[2]
###################################################

# Signifikanz. Bei 2 Stufen und einen Faktor:
summary(a)
# das gleiche
summary(glht(a, linfct=mcp(Alter="Tukey")))
###################################################

mes = read.table(url("http://www.phonetik.uni-muenchen.de/~jmh/lehre/sem/ss11/statfort/mes.txt"))

tab = with(mes, table(vtype, Pre))
barchart(tab, auto.key=T)
o = lmer(Pre ~ vtype + (1|spk), family=binomial, data=mes)
###################################################

summary(o)
neuv = relevel(mes$vtype, "o")
o2 = update(o, ~ . -vtype + neuv)
summary(o2)
###################################################

summary(glht(o, linfct=mcp(vtype="Tukey")))
###################################################
tab = with(mes, table(vtype, ptonic, Pre))
barchart(tab, auto.key=T)
###################################################

o = lmer(Pre ~ vtype * ptonic + (1|spk), family=binomial, data=mes)
ohne = lmer(Pre ~ vtype + ptonic + (1|spk), family=binomial, data=mes)
anova(o, ohne)

plabs = factor(with(mes, paste(mes$vtype, mes$ptonic)))
beide = lmer(Pre ~ plabs + (1|spk), family=binomial, data=mes)
summary(glht(beide, linfct=mcp(plabs = "Tukey")))
###################################################
summary(o)
lab = relevel(mes$vtype, "o")
o2 = lmer(pre ~ lab * ptonic + (1|spk), family=binomial, data=mes)
summary(o2)
###################################################
# 2 Faktoren einer numerisch
bwplot(Pre ~ clodur | ptonic, data = mes)
###################################################

o = lmer(Pre ~ clodur * ptonic + (1|spk), family=binomial, data=mes)
ohne = lmer(Pre ~ clodur + ptonic + (1|spk), family=binomial, data=mes)
anova(o, ohne)
o2 = update(ohne, ~ . - ptonic)
anova(o, o2)

o3 = update(o, ~ . - clodur)
anova(o, o3)
###################################################
bwplot(Pre ~ clodur | vtype, data = mes)

###################################################

o = lmer(Pre ~ clodur * vtype + (1|spk), family=binomial, data = mes)
o2 = lmer(Pre ~ clodur + vtype + (1|spk), family=binomial, data = mes)
anova(o, o2)

###################################################

temp = mes$vtype == "e"
mes.e = mes[temp,]
e = lmer(Pre ~ clodur + (1|spk), family=binomial, data = mes.e) 
summary(e)
###################################################

histogram(Pre ~ clodur | vtype, data = mes)
###################################################
summary(o)
lab = relevel(mes$vtype, "o")
summary(o3)
###################################################

yuki = read.table(url("http://www.phonetik.uni-muenchen.de/~jmh/lehre/sem/ss11/statfort/yuki.txt"))
###################################################

prop = with(yuki, Aus/(Aus+Ers))
yuki = cbind(yuki, prop)
bwplot(prop ~ Stim | Typ, data = yuki, horizontal=F)
###################################################

o = lmer(cbind(Aus, Ers) ~ Stim * Typ + (1|Vpn), family=binomial, data = yuki)
o2 = lmer(cbind(Aus, Ers) ~ Stim + Typ + (1|Vpn), family=binomial, data = yuki)
anova(o, o2)
###################################################

temp = yuki$Typ == "L"
yukiL = yuki[temp,]
yukiH = yuki[!temp,]
###################################################

H = lmer(cbind(Aus, Ers) ~ Stim + (1 + Stim | Vpn), family=binomial, data = yukiH)
w.H = coef(H)[[1]]
k.H = w.H[,1]
m.H = w.H[,2]
u.H = -k.H/m.H
vpn.H = rownames(w.H)

L = lmer(cbind(Aus, Ers) ~ Stim + (1 + Stim | Vpn), family=binomial, data = yukiL)
w.L = coef(L)[[1]]
k.L = w.L[,1]
m.L = w.L[,2]
u.L = -k.L/m.L
vpn.L = rownames(w.L)
###################################################

all(vpn.L == vpn.H)
par(mfrow=c(1,2))
boxplot(u.L - u.H, main = "Umkipppunkt (Unterschied)")
boxplot(m.L - m.H, main = "Neigung (Unterschied)")
t.test(u.L, u.H, paired=T)
t.test(m.L, m.H, paired=T)
###################################################

K.L = mean(k.L)
M.L = mean(m.L)
U.L = mean(-K.L/M.L)
Prop.L = with(yukiL, tapply(prop, Stim, mean))

K.H = mean(k.H)
M.H = mean(m.H)
U.H = mean(-K.H/M.H)
Prop.H = with(yukiH, tapply(prop, Stim, mean))
###################################################

ylim = c(0,1)
xlim = c(1,11)
curve(exp(M.L * x + K.L) / (1 + exp(M.L * x + K.L)), xlim=xlim, ylim=ylim, ylab="Aussage (Proportion)", xlab = "Stimulus Nummer")
points(1:11, Prop.L)
abline(v=U.L, lty=2)

par(new=T)
curve(exp(M.H * x + K.H) / (1 + exp(M.H * x + K.H)), add=T, col=2)
points(1:11, Prop.H, col=2)
abline(v=U.H, lty=2, col=2)