library(MASS)
library(lme4)
library(lattice)
library(ez)
library(multcomp)
source(file.path(pfadu, "phoc.txt"))
source(file.path(pfadu, "sigmoid.txt"))

# 1.
# 10 Sprecher produzierten Silben zu einer normalen und
# langsameren Sprechgeschwindigkeit. Die Silbendauern 
# waren wie folgt:

# normal
norm = c(165, 178, 143, 187, 186, 127, 138, 155, 157, 171)
# langsam
lang = c(182, 189, 179, 196, 188, 153, 154, 178, 169, 191)
boxplot(norm - lang)
t.test(norm - lang)
# Hat die Sprechgeschwindigkeit einen Einfluss auf die Silbendauer?

# Die Silbendauer wurde von der Sprechgeschwindigkeit
# beeinflusst (t[9] = 5.6, p < 0.001)

# 2.
form = read.table(file.path(pfadu, "F2dj.txt"))
#  30 Sprecher, 15 deutsch und 15 japanisch (Faktor L) produzierten jeweils
# /i:/ Vokale in einem labialen, alveolaren, oder velar Kontext
# (Faktor Stelle). F2 wurde gemessen (F2). Inwiefern wird
# F2 von der Sprache und/oder Artikulationsstelle beeinflusst?
boxplot(F2 ~ L * Stelle, data = form)
# vel > alv = lab; D und J unterscheiden sich nur in lab
ezANOVA(form, .(F2), .(Vpn), .(Stelle), .(L))
# F2 wurde von der Sprache (F[1,28] = 32.7, p < 0.001)
# und von der Stelle F[1.2, 34.6] = 91.6, p < 0.001)
# signifikant beeinflusst und es gab eine signifikante
# Interaktion zwischen diesen Faktoren (F[1.2, 34.6] = 3.9, p < 0.05).

p = phoc(form, .(F2), .(Vpn), .(Stelle, L))
round(phsel(p$res), 3)
round(phsel(p$res, 2), 3)

# Post-hoc Bonferroni korrigierte
# t-tests zeigten signifikante Unterschiede zwischen labial und alveolar für
# deutsch (p < 0.001) aber nicht für japanisch; zwischen
# labial und velar für deutsch (p < 0.001) und für japanisch (p < 0.001);
# und zwischen alveolar und velar für deutsch (p < 0.01) und
# für japanisch (p < 0.001). Es gab signifikante Unterschiede
# zwischen deustch und japanisch für labial (p < 0.001) jedoch nicht
# für die anderen zwei Artikulationsstellen.

# 3.
# Die folgenden Daten zeigen F1 (Hz) und Dauerwerte (ms) für verschiedene
# Vokale. 

F1 = c(576, 370, 612, 1216, 409, 1502, 946, 998, 189, 787, 210, 737)

dauer = c(178, 138, 94, 278, 158, 258, 198, 188, 98, 179, 138, 98)
plot(dauer, F1)
o = lm(F1 ~ dauer)
summary(o)

# Inwiefern kann F1 aus der Dauer vorhergesagt werden?

# F1 kann aus der Dauer vorhergesagt werden
# (R^2 = 0.63, F[1, 10] = 17.2, p < 0.01)

# 4.
ital = read.table(file.path(pfadu, "ital.txt"))
# Prüfen Sie, inwiefern die Wahl
# ob präaspiriert wird oder nicht (Faktor Pre)
# von der Region (Faktor Reg) abhängig ist.
tab = with(ital, table(Reg, Pre))
p = prop.table(tab, 1)
barchart(p, auto.key=T, horizontal=F)
o = lmer(Pre ~ Reg + (1|Vpn) + (1|Wort), family=binomial, data = ital)
o2 = update(o, ~ . -Reg)
anova(o, o2)
# Region hat einen signifikanten Einfluss auf die 
# ±Präaspiration-Verteilung (c^2[2] = 47.3, p < 0.001).
summary(glht(o, linfct = mcp(Reg = "Tukey")))
# Post-hoc Tukey Tests zeigten signifikante Unterschiede
# zwischen S und N (p < 0.001) und zwischen S und C (p < 0.001)
# jedoch nicht zwischen N und C.

# 5.
walt = read.table(file.path(pfadu, "walt.txt"))
# Ein Kontinuum wurde synthetisiert zwischen englischem 'will' und 'wool'
# in 13 Schritten (Faktor Stim) und 15 Hörer mussten pro 
# Stimulus beurteilen (Faktor Urteil), ob sie 'will' oder 'wool' hörten.

# (i) Bilden Sie eine Sigmoid-Funktion der Urteile ab, und überlagern
# Sie den Bevölkerungsumkipppunkt.
tab = with(walt, table(Stim, Urteil))
p = prop.table(tab, 1)
o = lmer(Urteil ~ Stim + (1+Stim|Vpn), family=binomial, data = walt)
k = fixef(o)[1]
m = fixef(o)[2]
sigmoid(p, k, m)
abline(v= -k/m)

# (ii) Berechnen Sie den hörerspezifischen Umkipppunkt (einen 
# Umkipppunkt pro Hörer).
walt.vpn = coef(o)[[1]]
walt.u = -walt.vpn[,1]/walt.vpn[,2]

# 6.
wjung = read.table(file.path(pfadu, "wjung.txt"))
# (i) 18 junge Versuchspersonen nahmen an dem selben Test Teil wie in 5.
# Berechnen Sie den hörerspezifischen Umkipppunkt (einen 
# Umkipppunkt pro Hörer) für diese Daten.
o = lmer(Urteil ~ Stim + (1+Stim|Vpn), family=binomial, data = wjung)
wjung.vpn = coef(o)[[1]]
wjung.u = -wjung.vpn[,1]/wjung.vpn[,2]

# (ii) Vergleichen Sie die Daten aus 5(ii) und 6(i), 
# um zu prüfen, ob die Hörergruppe (alt vs. jung)
# einen Einfluss auf die Umkipppunkte hat.
alter = factor(c(rep("A", length(walt.u)), rep("J", length(wjung.u))))
u = c(walt.u, wjung.u)
boxplot(u ~ alter)
t.test(u ~ alter)
# Alter hat einen signifikanten Einfluss auf den Umkipppunkt
# (t[31] = 5.2, p < 0.001).