library(car)
library(languageR)
library(multcomp)
pfad = "/Volumes/Data_1/d/"
attach(paste(pfad, "anova1", sep=""))
soa = read.table(paste(pfad, "soa.txt", sep="/"))
vot = read.table(paste(pfad, "vot.txt", sep="/"))
alc = read.table(paste(pfad, "alcdata.txt", sep="/"))

attach(soa)
boxplot(F1 ~ Pfinal * W, ylab="F1 (Hz)")

### lmer() Durchfuehrung I

# RM-Manova
soa.t = Anova.prepare(soa, c("s", "w", "w", "d"))
soa.lm = lm(soa.t$d ~ 1)
Anova(soa.lm, idata=soa.t$w, idesign = ~ W * Pfinal)

# Mixed model
F1.lmer = lmer(F1 ~ Pfinal + (1 | Vpn) + (1 | W), data=soa)

# F1.lmer = lmer(F1 ~ Pfinal + (1 + Pfinal | Vpn) + (1 +Pfinal | W))

print(F1.lmer, corr=F)
ranef(F1.lmer)
fitted(F1.lmer)

levels(Pfinal)
Pfinal2 = relevel(Pfinal, "short")
F1b.lmer = lmer(F1 ~ Pfinal2 + (1 | Vpn) + (1 | W))
print(F1b.lmer, corr=F)

2 *  ( 1 - pt(3.419, 16))

# Markov-Chain Monte-Carlo sampling
F1.fnc = pvals.fnc(F1.lmer)
# Wegen dem Bug
F1.lmer = lmer(F1 ~ Pfinal + (1 | Vpn) + (1 | W))

# Man soll *nicht* den Pr(>|t|) Wert nehmen, 
# weil die Freiheitsgrade nicht genau eingeschaetzt
# werden koennen sondern
# den pMCMC Wert.
F1.fnc$fixed

detach(soa)
### lmer() Durchfuehrung II
# d: VOT Werte von /b, d, g/. Inwiefern haben
# Artikulationsstelle und Gender einen Einfluss auf VOT?
attach(vot)
table(Wort, Vpn)

# Formel
d.lmer = lmer(d ~ Cons * G + (1|Wort) + (1|Vpn), data=vot)

print(d.lmer, corr=F)

# in p-Werte umwandeln
d.fnc = pvals.fnc(d.lmer, withMCMC=T)

d.fnc$fixed

# Wir brauchen jetzt einen F-table mit MCMC
# Werten, um den Signikanz vom *Faktor* zu pruefen

# Wegen dem Bug, noch einmal
d.lmer = lmer(d ~ Cons * G + (1|Wort) + (1|Vpn), data=vot)

# Wo kommen Consd und Consg vor?
names(d.fnc$mcmc)
d.aov = aovlmer.fnc(d.lmer, d.fnc$mcmc, 2:3)

# auf den $MCMC$p Wert schauen...
d.aov$MCMC$p
# Die Artikulationsstelle hatten einen signifikanten
# Einfluss auf VOT  (F = 11.25, MCMC  p < 0.05)

# Problem: wie bekomme ich einen p-Wert fuer die Interaktion?
# Ein Moeglichkeit. lmer() noch einmal ohne Interaktion
# ausfuehren, und dann anova() einsetzen, um
# die Wahrscheinlichkeit festzulegen, dass
# sich die beiden Mixed-Models unterscheiden

d.lmerohne = lmer(d ~ Cons + G + (1|Wort) + (1|Vpn), data=vot)
anova(d.lmer, d.lmerohne)

# Naechstes Problem. Muesste ich nicht den selben
# MCMC-p-Wert bekommen, wenn ich zB den Faktor neu kodiere
# und den lmer() mit den Faktoren in die andere Reihenfolge
# durchfuehre? zB

Cons2 = relevel(Cons, "d")
d.lmer2 = lmer(d ~  G * Cons2  + (1 | Wort) + (1 | Vpn))
d.fnc2 = pvals.fnc(d.lmer2, withMCMC=T)
names(d.fnc2$mcmc)
d.lmer2 = lmer(d ~  G * Cons2  + (1 | Wort) + (1 | Vpn))
d.aov2 = aovlmer.fnc(d.lmer2, d.fnc2$mcmc, 3:4)
# $MCMC$p jetzt [1] 0.0289 - statt  0.025


####### 3. Mixed Models und post-hoc tests

# post-hoc t-tests Es gibt eine Interaktion zwischen
# Gender und Artikulationsstelle. Wie pruefe ich,
# ob Interaktionen vorliegen? 
# Vielleicht lmer() nochmal durchfuehren mit
# G und Cons2 (oder Cons) zusammengelegt?
both = factor(paste(G, Cons, sep="."))
table(both)
post.lmer = lmer(d ~  both  + (1 | Wort) + (1 | Vpn))
# NB, F.b ist Basis
levels(both)
print(post.lmer, corr=F)

# Matrix der moeglichen Kontraste erstellen
K = diag(length(fixef(post.lmer)))
n = names(fixef(post.lmer))
rownames(K)= colnames(K) = n

# prueft ob both.F.d und both.F.g sich
# von der Basis-Stufe, both.F.b unterscheiden
plot(confint(glht(post.lmer,linfct= K[2:3,])))

# NB Bonferroni Korrektur muesste dann angewandt werden
summary(glht(post.lmer, linfct = K[2:3,]), test = adjusted("none"))

# Wie koennten wir pruefen, ob F.d sich von M.d unterscheiden?
# In dem wir die entsprechenden Reihen von K *subtrahieren*
K = rbind(K, K[2,]-K[5,])
rownames(K)[7] = "F.d - M.d"
plot(confint(glht(post.lmer,linfct= K[c(2:3,7),])))
summary(glht(post.lmer, linfct = K[c(2:3,7),]), test = adjusted("none"))

## 4. lmer() und logistische Regression #########################
# Haben Mode und Alkoholisierung
#  einen Einfluss auf die Hesitationen 
# 4.1 Ohne den Sprecher als Random Factor zu erklaeren (falsch)
hes = as.matrix(alc[,1:2])
alc.g = glm(hes ~ Mode * Alc, "binomial", data=alc)
anova(alc.g, test="Chisq")

# 4.2 mit lmer() und Sprecher als Random Faktor
alc.lmer = lmer(hes ~ Mode * Alc + (1 | Vpn), family="binomial", data=alc)

# Funktioniert: Bedeutung?
print(alc.lmer, corr=F)

# ich bekomme zum Glueck das gleiche Ergebnis in der anderen Reihenfolge
print(lmer(hes ~ Mode * Alc + (1 | Vpn), family="binomial", data=alc), corr=F)

# Funktioniert nicht -- kann die Signifikanzen von Faktoren nicht berechnen!!
anova(alc.lmer, test="Chisq")

# Funktioniert...auch nicht
alc.fnc = pvals.fnc(alc.lmer)


#######
# zusaetzliche Daten von Katalin
vr = read.table(paste(pfad, "vr.txt", sep="/"))
vr[,3] = factor(vr[,3])
attach(vr)
boxplot(ratio ~ stress * vowel)

# repeated measures Manova
table(subj, vowel, stress)
vr.t = Anova.prepare(vr, c("s", "w", "w", "d"))
vr.lm = lm(vr.t$d ~ 1)
vr.rmaov = Anova(vr.lm, idata=vr.t$w, idesign= ~vowel * stress)


# lmer

vr.lmer = lmer(ratio ~ vowel * stress + (1|subj))
print(vr.lmer, corr=F)
vr.pnc = pvals.fnc(vr.lmer, withMCMC=T)
vr.lmer = lmer(ratio ~ vowel * stress + (1|subj))
names(vr.pnc$mcmc)
vr.aov = aovlmer.fnc(vr.lmer, vr.pnc$mcmc, 2:3)

detach(vr)
######################
vrs = read.table(paste(pfad, "vrs.txt", sep="/"))

vrs[,3] = factor(vrs[,3])
vrs[,4] = factor(vrs[,4])

vrs.m = anova.mean(vrs[,6], vrs[,1], vrs[,2], vrs[,3], vrs[,4], vrs[,5])
colnames(vrs.m) = c("ratio", "Subj", "V", "S", "N", "Q")

attach(vrs.m)
# wird nicht gehen: einige Zellen sind nicht vollstaendig
table(V, S, N, Q)
detach(vrs.m)

# Nochmal, Faktor S weglassen
vrs.m = anova.mean(vrs[,6], vrs[,1], vrs[,2],  vrs[,4], vrs[,5])
colnames(vrs.m) = c("ratio", "Subj", "V", "N", "Q")

# RM-Manova
vrs.t = Anova.prepare(vrs.m, c("d", "s", "w", "w", "w"))
vrs.lm = lm(vrs.t$d ~ 1)
vrs.aov = Anova(vrs.lm, idata=vrs.t$w, idesign= ~ V * N * Q)

# lmer
colnames(vrs) = c("Subj", "V", "S", "N", "Q", "ratio")
vrs.lmer = lmer(ratio ~ V * N * Q + (1|Subj), data=vrs)
vrs.fnc = pvals.fnc(vrs.lmer, withMCMC=T)
vrs.lmer = lmer(ratio ~ V * N * Q + (1|Subj), data=vrs)
names(vrs.fnc$mcmc)
vrs.loav = aovlmer.fnc(vrs.lmer, vrs.fnc$mcmc, 2:3)
vrs.aov