library(lattice)
library(ez)
library(lmerTest)
library(multcomp)
source(file.path(pfadu, "phoc.txt"))
source(file.path(pfadu, "sig.fn.R"))
source(file.path(pfadu, "ph.step.R"))


#AV=Numerisch, UV=Numerisch, keine Random Factors
# aus f2.R, Aufgabe 2

# 2. Diese Daten:
tv = read.table(file.path(pfadu, "tv.df.txt"))
dim(tv)
# zeigen eine Messung der Sprachkompetenz (cdi) für 80 Schüler (Vpn) 
# sowie die Anzahl der Stunden am Tag, die sie fernsehen (tv.hours). 
# Prüfen Sie durch eine Abbildung und statistischen Test, 
# ob eine Verbindung zwischen Sprachkompetenz und Fernsehen vorliegt.



plot(cdi ~ tv.hours, data = tv)
reg = lm(cdi ~ tv.hours, data = tv)
summary(reg)
# Es gibt eine lineare Beziehung zwischen cdi und tv.hours (R^2 = 0.23, F[1,78] = 23.1, p < 0.001).
#dem plot hinzufügen:
abline(reg)


########################## oder mit xyplot() in library(lattice)
#zunächst Funktion einlesen:
mypanel = function(x, y, ...) {
  panel.xyplot(x, y, ...)
  panel.lmline(x, y)
}
#und dann
xyplot(cdi ~ tv.hours,data=tv, panel = mypanel)


summary(reg)
# Es gibt eine lineare Beziehung zwischen cdi und tv.hours (R^2 = 0.23, F[1,78] = 23.1, p < 0.001).
shapiro.test(resid(reg))
# OK
plot(resid(reg))
abline(h=0)
# OK
acf(resid(reg))
# OK




#AV=NUMERISCH, UV=KATEGORIAL, keine Random Factors
#aus f2.R, Aufgabe 3
# 3. Für diese Daten:
zweit = read.table(file.path(pfadu, "zweit.df.txt"))
dim(zweit)
head(zweit)
# nahmen Versuchspersonen (Vpn) an einem Test in einer 
# zweiten Sprache teil (l2score). Prüfen Sie durch eine Abbildung und statistischen Test, 
# ob l2score durch Geschlecht (G) beeinflusst wird.

bwplot(l2score ~ G, data = zweit)
with(zweit, table(Vpn, G))#logischerweise keine Paarung möglich
table(zweit$G)

#Vorsicht: Obwohl Vpn erwähnt wird, ist das hier KEIN Random Factor, 
# weil es für jede Vpn nur eine Beobachtung gibt
table(zweit$Vpn)

#deshalb t.test:
t.test(l2score ~ G, data = zweit)
# Geschlecht hat einen signifikanten Einfluss auf l2score (t[117.6] = 2.1, p < 0.05).
#wir müssen noch die Voraussetzung prüfen, dass beide Datensätze normalverteilt sein sollen 
with(zweit, tapply(l2score, G, shapiro.test))


ezANOVA(zweit,.(l2score),.(Vpn),between=.(G))
ezANOVA(zweit,.(l2score),between=.(G))


#AV=NUMERISCH, UV=KATEGORIAL, ein Random Factor
# 5. Prüfen Sie durch eine Abbildung und statistischen Test für diese Daten:
votspan = read.table(file.path(pfadu, "votspan.txt"))
dim(votspan)
head(votspan)

# inwiefern VOT vom Alter und von der Stadt beeinflusst wird.

bwplot(VOT ~ Alter | Stadt, data = votspan)


table(votspan$Vpn)
#auch hier gibt es nur EINE Beobachtung pro Vpn, aber wir haben zwei Faktoren
#(logischerweise between-Vpn-Factoren, weil es nur EINE Beobachtung pro Vpn gibt)
# so dass wir keinen t.test() anwenden können; also ezANOVA() mit Vpn als random factor:
with(votspan, table(Vpn, interaction(Alter, Stadt)))
ezANOVA(votspan, .(VOT), .(Vpn), between = .(Alter, Stadt))
# VOT wird signifikant vom Alter (F[1,44] = 46.2, p < 0.001) und von 
# der Stadt (F[1,44] = 17.5, p < 0.001) signifikant beeinflusst, 
# und es gibt eine signifikante (F[1, 44] = 5.7, p < 0.05) Interaktion zwischen diesen Faktoren.
# wegen Interaktion post-hoc:
p = phoc(votspan, .(VOT), .(Vpn), .(Alter, Stadt))
round(phsel(p$res), 3)
round(phsel(p$res, 2), 3)
# Bonferroni-korrigierte post-hoc t-tests bestätigten die signifikanten Unterschiede 
# zwischen alt und jung in Granada (p < 0.05) und in Sevilla (p < 0.001). 
# Der Unterschied zwischen den Städten war für jung (p < 0.01) 
# jedoch nicht für alt signifikant.


round(ezANOVA(votspan, .(VOT), .(Vpn), between=.(Alter, Stadt))$ANOVA$p,3)


#AV=NUMERISCH, UV=KATEGORIAL, mehr als ein Random Factor
#z.B. Aufgabe 7 aus f2.R:
# 7. Diese Daten:
phr = read.table(file.path(pfadu, "phr.df.txt"))
dim(phr)
# zeigen eine Messung der Sprechgeschwindigkeit (tempo) in 
# verschiedenen prosodischen Kontexten (Faktor Kontext) und 
# in verschiedenen Äußerungslängen (Faktor L). 
# Die Sprechgeschwindigkeitsmessungen wurden von 30 Sprechern (Vpn) 
# in 10 verschiedenen Äußerungen (Ag) erhoben. 
# Prüfen Sie durch eine Abbildung und statistischen Test, 
# inwiefern die Sprechgeschwindigkeit vom 
# Kontext und von der Äußerungslänge beeinflusst wird.

bwplot(tempo ~ Kontext, data = phr)
#Wir haben 2 Random Factors
#Vpn und Ag
#und brauchen daher ein mixed model mit lmer()
# within: (Kontext|Vpn)
with(phr, table(Vpn, Kontext))
# within: (Kontext|Ag)
with(phr, table(Ag, Kontext))
# within (L|Vpn)
with(phr, table(Vpn, L))
# between: (1|Ag)
with(phr, table(Ag, L))
# zusammen: (L+Kontext|Vpn) + (Kontext|Ag)

kon.lmer = lmer(tempo ~ Kontext * L  + (L+Kontext|Vpn) + (Kontext|Ag), data = phr)
kon.lmer = lmer(tempo ~ Kontext + L + Kontext:L + (L+Kontext|Vpn) + (Kontext|Ag), data = phr)
(kon.step = step(kon.lmer))
# Fixed effects:
#  Sum Sq  Mean Sq NumDF DenDF  F.value elim.num Pr(>F)
#Kontext:L   3.1213   0.7803     4   855   0.7890        1 0.5324
# L           1.2843   0.6422     2     7   0.6500        2 0.5509
#Kontext   488.6540 244.3270     2   859 247.2958       kept <1e-07

# Sprechgeschindigkeit wurde vom Kontext signifikant beeiflusst (F[2, 859] = 247.3, p < 0.001), aber nicht von der
# Äußerungslänge.

phr$Kontext # hat drei Stufen, daher ph.step
ph.step(kon.step,"b")
# Bonferroni-korrigierte post-hoc t-Tests zeigten, dass sich sowohl ip von paren, ip von pp und paren von pp 
# unterschieden (jeweils p < 0.001)



#AV=NUMERISCH, UV1=KATEGORIAL, UV2=NUMERISCH, mehr als ein Random Factor
#z.B.: letzte Aufgabe aus mmuebung.R 
#  In diesen Daten
lex = read.table(file.path(pfadu, "lex.txt"))
# wurden Reaktionszeiten von verschiedenen Sprechern (Subj) auf verschiedene Wörter (Word) gemessen. 
#Die Sprecher waren entweder Muttersprache englisch oder nicht (Faktor Lang). 
# Zusätzlich gibt es eine Verschlüsselung für die Wortlänge (Length): je kleiner die Zahl, 
# umso kürzer das Wort.  I
# nwiefern werden die Reaktionszeiten von der Sprache und  der Wortlänge beeinflusst?

#Die Besonderheit hier:
# während
lex$Lang
# ein steistufiger Faktor ist, ist
lex$Length
# eigentlich numerisch. Trotzdem kann und muss das mit lmer() gelöst werden ...

bwplot(rt~Lang,data=lex)
xyplot(rt~Length,data=lex)


with(lex, table(Lang, Subj))#-->(1|Subj)
with(lex, table(Lang, Word))#-->(Lang,Word)

with(lex, table(Length, Subj))#-->(Length|Subj)
with(lex, table(Length, Word))#-->(1|Word)

# in zwei Tests aufdröseln:
sprache.lmer = lmer(rt ~ Lang + (1|Subj) + (1+Lang|Word), data = lex)
lmerTest::step(sprache.lmer) 
lmerTest::anova(sprache.lmer)
# Die Reaktionszeiten wurden von Sprache beeinflusst (F[1,19.2]=6.3, p < 0.05)
#und
length.lmer = lmer(rt ~ Length + (Length|Subj) + (1|Word), data = lex)
step(length.lmer) 
# Außerdem wurden die Reaktionszeiten von der Länge beeinflusst (F[1,53.1]=14.7, p < 0.001)

# mehr müssen wir hier nicht beantworten
# nur als Fleißaufgabe
# geht es aber durchaus auch in einem Modell:
lex.lmer = lmer(rt ~ Lang * Length + (Length|Subj) + (Lang|Word), data = lex)
step(lex.lmer)
# Die Reaktionszeiten wurden von Sprache beeinflusst (F[1,19.0]=6.3, p < 0.05),
# die auch von der Länge (F[1,57.2]=19.1, p < 0.001)
# Außerdem gab es eine signifikante Interaktion zwischen Länge und Sprache (F[1,18.5]=8.3, p< 0.01).

# Da Length numerisch ist, können wir hier natürlich keinen Post-hoc-t-Test rechnen.
# (es sei denn wir würden LEngth mit factor(Length) künstlich zu einem 8-Kategorien-Faktor machen,
# was wir uns an dieser Stelle ersparen)
# Dieser plot hier
xyplot(rt~Length|Lang,data=lex, type = c("p", "smooth"))
# zeigt uns jedoch, dass der Einfluss von Length bei Lang=="Other" offenbar stärker ausgeprägt ist
# als bei Lang=="English", da bei Other der Anstieg der rt-Werte mit zunemender "Length" 
# deutlicher ist als bei English.







#AV=KATEGORIAL, UV=KATEGORIAL, kein Random Factor
# z.B. Aufgabe 5 aus muster.R
kj = read.delim(file.path(pfadu, "kj.txt"))

# 5. Prüfen Sie für diese Daten mit einer Abbildung und einem statistischen Test:
head(kj)
# inwiefern die Wahl des Frikatives (Faktor fric) als 's' (alveolar) oder 'S' (retroflex) 
# von der Emphase (Faktor emphatic) beeinflusst wird.
#wieviele less emphatic s, more emphatic s, less emphatic S, und more emphatic S gibt es?
tab = with(kj, table(emphatic, fric))
# pro Frikativkategorie Proportionen berechnen
p = prop.table(tab, 1)
barchart(p, auto.key=T, horizontal=F)
e.glm = glm(fric ~ emphatic, binomial, data = kj)
anova(e.glm,test="Chisq")
# Die Wahl des Frikatives wird signifikant von der Emphase beeinflusst (X^2[1]=20.1, p < 0.001).






#AV=KATEGORIAL, UV=NUMERISCH, kein Random Factor
#z.B. Aufgabe 4 aus f2.R
# 4. Für diese Daten:
elwi = read.table(file.path(pfadu, "elwi.df.txt"))
dim(elwi)
# wurde ein 13-stufiges F2-Kontinuum erstellt (F2) und 
# jede Versuchsperson musste pro Stimulus entscheiden, 
# ob 'will' oder 'wool' wahrgenommen wurde (Urteil). 
# Prüfen Sie durch eine Abbildung und statistischen Test, 
# ob die will/wool Entscheidung durch F2 beeinflusst wird. 
# Berechnen Sie den F2-Umkipppunkt, zu dem die Entscheidung 
# von 'will' nach 'wool' kippt, und überlagern Sie eine 
# Sigmoidalkurve mit der sig() Funktion auf die proportionalen Werte.
elwi$F2
# ist numerisch, wenn auch eine Art von Stufigkeit zu beobachten ist
levels(elwi$Urteil)
P = elwi$Urteil == "wool"
Q = !P
e.m = aggregate(cbind(P, Q) ~ F2, sum, data = elwi)
p = with(e.m, P/(P+Q))
e.m = cbind(e.m, p)
plot(p ~ F2, data = e.m)
wool.glm = glm(Urteil ~ F2, binomial, data = elwi)
# oder
wool.glm = glm(cbind(P, Q) ~ F2, binomial, data = e.m)


anova(wool.glm, test="Chisq")
# Das Urteil wird von F2 beeinflusst (X^2[1] = 65.1, p < 0.001)

cf = coef(wool.glm)
k = cf[1]; m = cf[2]
#mit den Koeffizienten die Sigmoidalkurve berechnen
sig(k, m, add=T)
#Umkipppunkt berechnen und überlagern
u = -k/m
abline(v=u, lty=2)

abline(h=0.5,lty=2)