# 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:

glottal = read.table(file.path(pfadu, "glottal.txt"))
fremd = read.table(file.path(pfadu, "fremd.df.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"))


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


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

# zeigen ob ein aspirierter (t) oder ein 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 + ylab("Proportionen")
#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, ganz leichte Tendenz dazu gibt ).




###################################################################################
# 2. 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
#aber zu gering (dim(fremd) zeigt 60, with(fremd, table(Vpn, Sprache)) jeweils nur eine Beobachtung)

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) + 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))

###################################################################################
# 3. 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)


get_model(vcv.step) #RT ~ Vowel + (Vowel | Subject) + (Vowel | Context)
anova(get_model(vcv.step))
# Die Daten zeigen einen signifikanten Einfluss vom Vokal (F[2,18.5] =  5.8, 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
vcv.ph = pairs(emmeans(get_model(vcv.step),  ~ Vowel))
#dann
phsel(vcv.ph,1)
#Ein tukey-korrigierter post-hoc Test zeigte, dass nur die Unterschiede in den Paarungen
# a vs. i (p<0.05) und a vs. u (p<0.05), aber nicht in der Paarung i vs. u signifikant sind.


###################################################################################
# 4. 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 ergibt 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).




###################################################################################
# 5. Fuer die Daten in diesem Data-Frame:
dim(svot)
# pruefen Sie mit einer Abbildung und einem 
# statistischen Text, inwiefern VOT von 
# Gender und Region beeinflusst wird. 
# (Vpn enthaelt das Sprecherkuerzel).


...

###################################################################################
# 6. Die voice-onset-time-Werte für /t/, 
# gesprochen von 16 Sprechern mit L1-Franzoesisch, 
# waren wie folgt:

vot = c(-4,28,23,1,-13,20,21,-15,4,29,16,22,35,3,3,-4)

# (a) Wie wahrscheinlich es ist, dass /t/ 
# voellig unaspiriert ist 
# ("voellig unaspiriert" = mit einem VOT-Wert von 0 ms)?
# Berichten Sie hierbei auch das 95%-Konfidenzintervall 
# für die Daten in vot!

# Folgend den Angaben in https://www.phonetik.uni-muenchen.de/~jmh/lehre/sem/ss19/statistikR_files/t-Test.html
# könnte man diese Frage nach dem 95%-Konfidenzintervall auf die folgende,
# äußerst komplizierte Weise lösen: verwende qt() anstelle von qnorm(),
# mit den Parametern 
# mu (= mean(vot) ), 
# SE (= sd(vot)/sqrt(length(vot)) ), 
# df (= Anzahl der Beobachtungen minus 1,
# also (length(vot)-1) ),
# und den Wahrscheinlichkeitswerten 2.5% (0.025) und 97.5% (0.975):
a = mean(vot) + sd(vot)/sqrt(length(vot)) * qt(0.025, df=length(vot)-1)
b = mean(vot) + sd(vot)/sqrt(length(vot)) * qt(0.975, df=length(vot)-1)
a
b
# Puh, kompliziert. Viel, viel einfacher ist:

#Wie wahrscheinlich es ist, dass /t/ 
# voellig unaspiriert ist? lässt sich übersetzen in 
# "die Wahrscheinlichkeit, dass die Werte in vot gleich 0 sind"
# --> klassischer one-sample t.test mit mu=0, also
t.test(vot)
# One Sample t-test
# 
# data:  vot
# t = 2.7137, df = 15, p-value = 0.01601
# alternative hypothesis: true mean is not equal to 0
# 95 percent confidence interval:
#   2.26633 18.85867
# sample estimates:
#   mean of x 
# 10.5625 

# t[15]=2.7, p < 0.05 (diesmal ausnahmsweise ganz genau: p = 0.01601);
# Die Wahrscheinlichkeit, dass /t/ völlig unaspiriert ist, liegt bei 1.6 Prozent.
# Das 95 %-Konfidenzintervall liegt zwischen 2.26633 und 18.85867 

# (diese Werte werden beim t-Test-Ergebnis angegeben!!!); vgl. die händisch errechneten
a
b

#
# (b) Die VOT-Werte derselben 16 franzoesischen Sprecher 
# nach einem 6 monatigen Aufenthalt in Deutschland waren wie folgt:

vot2 = c(2,28,26,4,-9,23,19,-6,11,33,24,31,41,10,7,8)



#
# Pruefen Sie anhand einer Abbildung und einem statistischen Test, 
# ob und falls ja, wie die voice-onset-time von dem Aufenthalt 
# in Deutschland beeinflusst wurde!

boxplot(vot2-vot)
shapiro.test(vot2-vot)#Okay!
t.test(vot2-vot)

# Ein gepaarter t-Test ergab, dass nach einem sechsmonatigem Aufenthalt in Deutschland
# französische Muttersprachler eine höhere VOT hatten (t[15] = 5.9, p < 0.001).