# Aehnlich wie in der Klausur: gemischte Fragen:
# pfadu = "http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf"
# Anmerkung: es werden immer alle packages und selbstgeschriebene
# Funktionen bereitgestellt werden.
# Das heisst aber nicht notwendigerweise, dass Sie die alle brauchen!
library(ggplot2)
library(dplyr)
library(ez)
library(lmerTest)
library(emmeans)
source(file.path(pfadu, "phoc.txt"))
source(file.path(pfadu, "sig.fn.R"))

proben <- function(unten = 1,
                   oben = 6,
                   k = 10,
                   N = 50)
{
  # default: wir werfen 10 Wuerfel 50 Mal
  alle <- NULL
  for (j in 1:N) {
    ergebnis = mean(sample(unten:oben, k, replace = T))
    alle = c(alle, ergebnis)
  }
  alle
}


#Dann die benoetigten Data.frames:
tv = read.table(file.path(pfadu, "tv.df.txt"))
fremd = read.table(file.path(pfadu, "fremd.df.txt"))
glottal = read.table(file.path(pfadu, "glottal.txt"))
vcv = read.table(file.path(pfadu, "vcvC.txt"))
form10 = read.table(file.path(pfadu, "form10.df.txt"))
svot = read.table(file.path(pfadu, "svot.txt"))
stefan = read.table(file.path(pfadu, "stefan.txt"))

# 1. Diese Daten:
dim(tv)
# zeigen eine Messung der Sprachkompetenz (cdi) fuer
# 80 Schueler (Vpn) sowie die Anzahl
# der Stunden am Tag, die sie fernsehen (tv.hours).
# Pruefen Sie durch eine Abbildung und statistischen Test,
# ob eine Verbindung zwischen Sprachkompetenz und
# der Anzahl der Fernsehstunden vorliegt.
plot(cdi ~ tv.hours, data = tv)
# Trotz der Nennung von "80 Schuelern" ist hier kein Mixed Model noetig,
# da fuer jede Versuchsperson nur eine einzige Beobachtung (also insgesamt 80) ...
dim(tv)
# ...vorhanden ist
# (d.h., die Daten sind unabhaengig voneinander und koennen deshalb durch
# Regression bearbeitet werden)
reg = lm(cdi ~ tv.hours, data = tv)
abline(reg)
#oder
ggplot(tv) +
  aes(y = cdi, x = tv.hours) +
  geom_point() +
  geom_abline(intercept = coef(reg)[1], slope = coef(reg)[2])
summary(reg)
# Es gibt eine lineare Beziehung zwischen cdi und tv.hours (R^2 = 0.22, F[1,78] = 23.1, p < 0.001).
shapiro.test(resid(reg))
# OK
plot(resid(reg))
# OK
acf(resid(reg))
# OK



# 2. Die Daten in diesem Data-Frame:
dim(glottal)

# zeigen ob ein aspirierter (t) oder glottalisierter (tQ)
# Verschluss von einem Erkennungssystem
# richtig (ja) erkannt wurde oder nicht (nein).
# Erstellen Sie eine Abbildung, und fuehren Sie
# einen statistischen Test durch, um zu pruefen,
# ob die Erkennung von dem Verschluss beeinflusst wurde.
head(glottal)
p1 = ggplot(glottal)
p2 = aes(fill = erkannt, x = verschluss)
p3 = geom_bar(position = "fill")
p1 + p2 + p3
#beide Variablen (die abhaengige und die unabhaengige) sind kategorial, daher glm()
g = glm(erkannt ~ verschluss, binomial, data = glottal)
anova(g, test = "Chisq")
# Die Erkennung wird vom Verschluss nicht signifikant beeinflusst,
# obwohl es eine ganz leichte Tendenz dazu gibt.




# 3. Diese Daten:
dim(fremd)
# zeigen fuer 30 Sprecher (Vpn) eine Messung der Sprechgeschwindigkeit (tempo),
# wenn sie in ihrer Muttersprache oder ihrer Zweitsprache (Sprache) reden.
# Pruefen Sie durch eine Abbildung und statistischen Test,
# inwiefern die Sprechgeschwindigkeit von der Sprache beeinflusst wird.

ggplot(fremd) +
  aes(y = tempo, x = Sprache) +
  geom_boxplot()

#Fuer den Test theoretisch drei Moeglichkeiten:
#t.test(), ezANOVA(), oder lmer(); fuer letzteres ist die Anzahl der Beobachtungen
#sehr gering (dim(fremd) zeigt 60 )

with(fremd, table(Vpn, Sprache))
# Sprache ist within; es gibt keine fehlenden Werte und keine Wiederholungen;
# eigentlich ein Fall fuer einen gepaarten t-test (den man auf zwei Weisen realisieren kann);
# der Vollstaendigkeit halber: man kann das auch als ezANOVA() mit einem within-Faktor berechnen:
ezANOVA(fremd, .(tempo), .(Vpn), .(Sprache))
# tempo wurde signifikant von der Sprache beeinflusst (F[1,29] = 28.2, p < 0.001).

# oder eben auf zweierlei Weise als gepaarter t-test:
fremddiff = fremd %>%
  group_by(Vpn) %>%
  summarise(tempo = diff(tempo))
#dann auch die bessere Abbildung:
boxplot(fremddiff$tempo)
abline(h = 0)
# oder komplizierter und trotzdem unschoener
ggplot(fremddiff) + aes(y = tempo, x = factor(0)) + geom_boxplot() + geom_hline(yintercept =
                                                                                  0)
shapiro.test(fremddiff$tempo)
# fremddiff$tempo ist normalverteilt, und daher duerfen wir einen t-Test anwenden:


t.test(fremddiff$tempo)
# tempo wurde signifikant von der Sprache beeinflusst (t[29] = 5.3, p < 0.001).
# oder 3. auf diese Weise (mit dem urspruenglichen data.frame):
t.test(tempo ~ Sprache, paired = TRUE, data = fremd)
# Es reicht natuerlich vollkommen EINER dieser drei moeglichen Loesungswege aus...
#(und die Hilfslinie mit abline() duerfen Sie auch getrost weglassen)

# ABER: auch wenn Sie die ezANOVA() benutzen, sollten Sie bitte die Differenzen berechnen
# und davon die Abbildung erstellen sowie die Normalverteilung der Differenzen testen
# (mit shapiro.test(fremddiff$tempo))

# 4. Diese Daten:
dim(vcv)

# zeigen Reaktionszeiten (RT) von verschiedenen Versuchspersonen (Subject)
# auf 3 Vokale (Faktor Vowel) in verschiedenen davor und danach
# kommenden Konsonanten-Kontexten (Context). Pruefen Sie durch eine Abbildung und
# statistischen Test, ob die Reaktionszeiten vom Vokal (Vowel) und/oder
# Sprachgruppe (Faktor Lang) beeinflusst wurden.

# RT ist numerisch; Vowel kategorial; mehrere Random Factors (Subject und Context);
# relativ viele Beobachtungen (dim(vcv) = 810;7), daher MM mit lmer() moeglich:
ggplot(vcv) +
  aes(y = RT, x = Vowel) +
  geom_boxplot() +
  facet_wrap( ~ Lang)

#oder:
ggplot(vcv) +
  aes(y = RT, x = Vowel, col = Lang) +
  geom_boxplot()
#oder:
ggplot(vcv) +
  aes(y = RT, x = Vowel, linetype = Lang) +
  geom_boxplot()
# oder beliebige Kombinationen davon

with(vcv, table(Subject, interaction(Vowel, Lang)))
with(vcv, table(Context, interaction(Vowel, Lang)))#oder einfacher aufgeteilt:

with(vcv, table(Subject, Vowel))#-->Vowel within in Bezug zu Subject
with(vcv, table(Subject, Lang))#--> Lang between in Bezug zu Subject

with(vcv, table(Context, Lang))#-->Lang within in Bezug zu Context
with(vcv, table(Context, Vowel))#-->Vowel within in Bezug zu Context


#Formel:
# RT~Vowel+Lang+Vowel:Lang-->RT~Vowel*Lang + (Vowel|Subject) + (Lang+Vowel|Context)


vcv.lmer = lmer(RT ~ Vowel * Lang + (Vowel |
                                       Subject) + (Vowel + Lang | Context), data = vcv)
vcv.step = step(vcv.lmer)

anova(lmer(RT ~ Vowel + (Vowel |
                           Subject) + (Vowel | Context), data = vcv))

get_model(vcv.step)
anova(get_model(vcv.step))
# Die Daten zeigen einen signifikanten Einfluss vom Vokal (F[2,18.6] =  5.75, p < 0.05)
# aber keine Einfluesse von der Sprachgruppe und es gab auch keine Interaktion zwischen diesen Faktoren.

# da Vowel dreistufig ist, kann man post-hoc testen,
# inwiefern sich a von i, i von u, und a von u unterscheiden
pairs(emmeans(get_model(vcv.step),  ~ Vowel))
#auch
phsel(pairs(emmeans(get_model(vcv.step),  ~ Vowel)))
#Ein tukey-korrigierter post-hoc Test zeigte, dass fuer die Unterschiede nur die Paarungen
# a vs. i (p<0.05) und a vs. u (p<0.05), aber nicht die Paarung i vs. u verantwortlich sind.


# 5. Fuer diese Daten:
dim(form10)
# pruefen Sie durch eine Abbildung und statistischen Test,
# ob F1 aus der Dauer vorhergesagt werden kann.
# Was waere der vorhergesagte F1-Wert fuer eine Dauer von 50 ms?
# Ist diese Vorhersage sehr verlaesslich?
plot(F1 ~ Dauer, data = form10)
form10.lm = lm(F1 ~ Dauer, data = form10)
abline(form10.lm)
#oder sehr umstaendlich
ggplot(form10) +
  aes(y = F1, x = Dauer) +
  geom_point() +
  geom_abline(intercept = coef(form10.lm)[1],
              slope = coef(form10.lm)[2])
summary(form10.lm)
predict(form10.lm, data.frame(Dauer = 50))
# Es gibt eine signifikante lineare Beziehung zwischen
# F1 und Dauer (R^2 = 0.05, F[1,198] = 12.0, p < 0.001)
#Das Modell sagt einen F1-Wert von 791.7138 fuer einen
#Dauerwert von 50 vorher; wegen des geringen R^2-Wertes macht diese
#Vorhersage aber so gut wie keinen Sinn.

shapiro.test(resid(form10.lm))
# OK
plot(resid(form10.lm))
# OK
acf(resid(form10.lm))
# OK (auch wenn bei Lag=2 der blaue Bereich minimal ueberschritten wird).