# Klausurvorbereitende Übung

# Vergessen Sie nicht, diese Textdatei regelmaessig zu speichern!

# Viel Erfolg!

# Hier die benoetigten Packages und Funktionen (bitte in R einlesen):
pfadu = "http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf"
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
}



# Hier die benoetigten data.frames:
hoch = read.table(file.path(pfadu, "hoch.df.txt"))
ctm = read.table(file.path(pfadu, "ctm.df.txt"))
vdata = read.table(file.path(pfadu, "vdata.k.4.txt"))


######################################################################
# FRAGEN (1 - 5) #####################################################
######################################################################

# 1. Fuer diese Daten:
dim(hoch)
# beantworten Sie durch eine Abbildung und statistischen Test die Frage,
# ob die Wahl zwischen /i/ und /e/ (Faktor V) von F1 (Faktor F1) beeinflusst wird. 
# Zu welchem F1-Wert kommt der Umkipppunkt zwischen /i/ und /e/ vor?
# Ueberlagern Sie auch eine Sigmoidalkurve auf ihre Abbildung! 
# Erlaeutern Sie schriftlich, was das Problem mit dem Umkipppunkt ist
# (wo liegt dieser Umkipppunkt?)!

levels(hoch$V)
P = hoch$V == "i"
hoch$P = hoch$V == levels(hoch$V)[2]
hoch$Q = !hoch$P

hoch.sum = hoch %>%
  group_by(F1) %>%
  summarise(P = sum(P), Q = sum(Q)) %>%
  mutate(p = P/(P+Q))

plot(p ~ F1, data = hoch.sum)

hoch.glm = glm(V ~ F1,family = binomial,data = hoch)
anova(hoch.glm,test="Chisq")

# F1 hat einen sig. Einfluss auf die Vokalkategorie
# (X^2[1]=4.7, p < 0.05 )
k=coef(hoch.glm)[1]
m=coef(hoch.glm)[2]
umkipp = -k/m
umkipp

sig(k,m,add = TRUE)

#Der Umkipppunkt ist bei 363 Hz (und damit außerhalb des 
#Variationsbereichs von F1 im df "hoch").
#####################################################################

# 2. Fuer diesen Data-Frame
dim(ctm)
# pruefen Sie durch eine Abbildung und einen statistischen Test, 
# inwiefern die F2-Werte von einem ersten Vokal (V1) von den F2-Werten 
# des danach kommenden Vokals (V2) beeinflusst werden. 
# Was ist der vorhergesagte V1-Wert fuer einen V2-Wert von 2600?

plot(V1~V2,data=ctm)
reg = lm(V1~V2,data=ctm)
summary(reg)
abline(reg)

shapiro.test(resid(reg))#Okay
plot(resid(reg))#fragwürdig
acf(resid(reg))#Okay

#V2 hatte einen sig. Einfluss auf V1 
# (R^2=0.44, F[1,37]=29.6, p < 0.001)

predict(reg,data.frame(V2 = 2600))
#Für einen V2-Wert von 2600 wird ein V1-Wert von 1964.3 vorhergesagt.

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

# 3.
# Wenn ich 50 Lose aus einem Hut mit den Zahlen -100 bis +100 ziehe, 
# a) wie hoch ist dann die Wahrscheinlichkeit, dass der Mittelwert der 50 Zahlen
# zwischen -5 und +5 liegt?
# b) fuehren Sie den obigen Vorgang 200 Mal mit der proben()-Funktion durch; 
# erstellen Sie dann noch ein Histogramm, 
# das die Mittelwerte ihrer 200 Versuche gezaehlt darstellt.

mu = mean(-100:100)
SE = sd(-100:100)*sqrt(200/201)/sqrt(50)

pnorm(5,mu,SE) - pnorm(-5,mu,SE)
#Die Wahrscheinlichkeit liegt bei 45.8 %

hist(proben(unten = -100,oben = 100,k = 50,N = 200))

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

# 4. In diesem Datensatz
head(vdata)
# sind Vokaldauern (dur) von gespannten und ungespannten (Faktor Tense) Vokalen
# aufgefuehrt, die von verschiedenen Sprechern (Faktor Vpn) in Pseudowoertern
# mit unterschiedlichen konsonantischen Kontexten (Faktor Cons) in 
# zwei Sprechgeschwindigkeiten (Faktor Rate) gesprochen wurden.
# Pruefen Sie durch eine Abbildung und einen statistischen Test,
# ob die Gespanntheit (Tense) der Vokale und die Sprechgeschwindigkeit (Rate)
# einen Einfluss auf die Vokaldauern (dur) hatten!
# (Vorsicht: die Rechendauer kann leicht erhoeht sein.)

ggplot(vdata) +
  aes(y = dur, x = Tense, col = Rate) +
  geom_boxplot()

dim(vdata)

with(vdata,table(Vpn,Tense))
# Tense ist within in Bezug zu Vpn
with(vdata,table(Vpn,Rate))
# Rate ist within in Bezug zu Vpn
with(vdata,table(Cons,Tense))
# Tense ist within in Bezug zu Cons
with(vdata,table(Cons,Rate))
# Rate ist within in Bezug zu Cons

vdata.lmer = lmer(dur ~ Tense*Rate + 
                    (Tense+Rate|Vpn) +
                    (Tense + Rate|Cons),
                  data = vdata)
vdata.step = step(vdata.lmer)
anova(get_model(vdata.step))

# Sowohl Tense (F[1,3.4]=29.8, p<0.01) 
# als auch Rate (F[1,3.0]=40.1,p<0.01) 
# hatten einen sig. Einfluss
# auf die Dauern; außerdem gab es eine
# sig. Interaktion zwischen Tense und Rate
# (F[1,1707.9] = 188.9, p < 0.001)


vdata.ph = pairs(emmeans(get_model(vdata.step),~Tense:Rate))

phsel(vdata.ph,1)
phsel(vdata.ph,2)
# Nur in der gespannten Bedingung unterscheiden sich langsam und schnell
# signifikant voneinender (p <0.01), nicht aber in der ungespannten;
# Nur in der langsamen Bedingung gibt es sig. Unterschiede zwischen
# gespannt und ungespannt (p < 0.05), nicht aber in der schnellen Bedingung.
######################################################################

# 5. Sind die Daten in
dat = c(-13,18,20,-6,17,28,19,23,21,2)

# signifikant > 0?
shapiro.test(dat)
wilcox.test(dat)

# Ja, die Daten sind signifikant größer als 0!
# (V = 50, p < 0.05)

# - ENDE - ###########################################################