library(lmerTest)
library(emmeans)
library(ggplot2)
library(dplyr)
library(ez)
source(file.path(pfadu, "phoc.txt"))

# 1.: Viele Daten (dauern laenger):
pfric = read.table(file.path(pfadu, "pfric2.df.txt"))

# 1. Fuer diese Daten (ein Beispiel mit sehr vielen Beobachtungen):
dim(pfric)
# Neun Versuchspersonen haben Woerter bzw. Logatome (Pseudowoerter)
# (=Faktor LW), die verschiedene Frikativkonsonanten enthielten
# (=Faktor K), produziert. Es wurden mehrere akustische Parameter
# gemessen, darunter der Wert "k0".
# Pruefen Sie durch eine Abbildung und einen 
# statistischen Test - ohne zu mitteln -, inwiefern 
# k0 von K (C,s, oder S), von LW (Logatom oder Wort), 
# und von Rate abhaengt.

ggplot(pfric) +geom_boxplot() + aes(y=k0,x=K,col=Rate) + facet_wrap(~LW)

with(pfric,table(K,Vpn))
with(pfric,table(LW,Vpn))
with(pfric,table(Rate,Vpn))
#das Folgende wird relativ lange dauern:
pfric.lmer = lmer(k0~K*LW*Rate + (K+LW+Rate|Vpn),data=pfric)
pfric.step = step(pfric.lmer)
pfric.step
anova(get_model(pfric.step))
# k0 war nicht davon abhaengig, ob es in einem Wort oder Logatom gemessen wurde,
# aber von der Frikativklasse K (F[2,8]=41.7, p<0.001) und von der Rate
# (F[1,8]=18.5, p< 0.01); ausserdem gab eine signifikante Interaktion zwischen
# K und Rate (F[2,3842.7]=3.7, p < 0.05).
#daher kann man einen neuen Plot machen, in dem man nicht zwischen 
#Logatom und Wort unterscheidet:
ggplot(pfric) +geom_boxplot() + aes(y=k0,x=K,col=Rate)
library(emmeans)
pfric.ph = pairs(emmeans(get_model(pfric.step),~K:Rate))
#aus der hier ausgegebenen Warnmeldung entnehmen wir bitte die beiden Befehle!
emm_options(pbkrtest.limit = 3880)
emm_options(lmerTest.limit = 3880)
#das hier dauert dann wirklich lange:
pfric.ph = pairs(emmeans(get_model(pfric.step),~K:Rate))
pfric.ph
phsel(pfric.ph,1)
phsel(pfric.ph,2)


# Post-hoc-t-Tests (Bonferroni-korrigiert) zeigten
# zwischen schnellem und langsamen C (p<0.05) und s (p<0.05), und S (p<0.05). 
# In der schnellen Bedingung
# unterschieden sich C und s (p < 0.05) sowie s und S (p < 0.001), jedoch nicht
# C und S.
# In der langsamen Bedingung unterschieden sich alle drei Paarungen, 
# also C vs. s, C vs. S (jeweils mit p<0.05), und s vs. S (p<0.001)

#2.: pairs(emmeans()) auch bei nur einem Faktor sinnvoll
# (wenn dieser mehr als zwei Stufen hat)
# 
x24 = read.table(file.path(pfadu, "x24.txt"))
# wurden fuer 3 Vokale (V) Entfernung zum Vokalmittelpunkt (ent)
# fuer verschiedene Sprecher (Vpn), die verschiedene Woerter (Wort) produziert haben.
# Inwiefern werden die Entfernungen vom Vokal beeinflusst?

summary(x24)
# Also
#AV: ent
#UVen: V (Vokal), Vpn, Wort
# hiervon interessant fuer uns: nur V -->fixed
# auszuklammern sind Vpn, Wort -->random

# Bild: ent in Abhaengigkeit von V
ggplot(x24) + aes(y=ent,x=V) + geom_boxplot()

# V ist also ein fixed factor
# Vpn, Wort sind random factors
# haben wir genuegend Daten fuer ein MM?
dim(x24)
# Ja! Wir brauchen zunaechst die library
library(lmerTest)
# wir muessen fuer beide random factors bestimmen
# ob V in Bezug zu ihnen within oder between ist:

with(x24, table(Vpn, V)) # --> within (V|Vpn)

with(x24, table(Wort, V)) # --> # within (V|Wort)
# unsere Formel lautet also: ent ~ V + (V|Vpn) + (V|Wort)
# MM
x24.lmer = lmer(ent ~ V + (V|Vpn) + (V|Wort), data = x24)
# Fuer die Pruefstatistik brauchen wir die step()-Funktion:
(x24.step = step(x24.lmer))

anova(get_model(x24.step))
# Die Entfernungen sind signifikant vom Vokal beeinflusst worden
# (F[2, 803] = 23.6, p < 0.001).

# Das Model wurde zwar auch vereinfacht 
# (V wurde bezueglich Vpn und Wort jeweils von within- zu between hin geaendert,
# d.h es wurde keine Vpn- bzw. Wort-spezifischen Steigungen fuer V berechnet), 
# das muessen Sie aber nicht berichten.



# Da es nur eine unabhaengige Variable gab, deren Einfluss uns interessierte,
# muessen wir hier eigentlich keinen Post-hoc-Test durchfuehren 
# (da es keinen zweiten Faktor gibt, mit dem V interagieren koennte).

# Aber: es gab drei Stufen!! Es koennte sein, dass sich nur a von 
# den beiden anderen Vokalen unterscheidet, aber nicht u vs. i (siehe Abbildung!)
# Um zu pruefen,
# wie die einzelnen Stufen von V sich
# unterscheiden, fuehren Sie pairs(emmeans()) mit nur einem
# Faktor nach der Tilde durch (wie gesagt: nur, sofern der Faktor 
# drei oder mehr Stufen hat):
pairs(emmeans(get_model(x24.step),~V))
# contrast   estimate         SE  df t.ratio p.value
# a - i    -0.4073304 0.05960399 803  -6.834  <.0001
# a - u    -0.2369580 0.05960399 803  -3.976  0.0002
# i - u     0.1703724 0.05960399 803   2.858  0.0121
# 
# P value adjustment: tukey method for comparing a family of 3 estimates 

# Die Entfernungen sind signifikant vom Vokal beeinflusst worden
# (F[2, 803] = 23.6, p < 0.001). Im Einzelnen unterschieden sich sowohl
# a von i (p<0.001), a von u (p<0.001), als auch im geringeren Masse 
# i von u (p<0.05)