# Einer nachmittaeglichen Diskussion mit Kommilitonen verdankt sich dieses 
# Zusatzbeispiel (vgl. Beispiel 4 aus a10.R)
library(ggplot2)
library(dplyr)
library(ez)
library(lmerTest)
library(emmeans)
source(file.path(pfadu, "phoc.txt"))

# aus f5.R
phr = read.table(file.path(pfadu, "phr.df.txt"))

# zeigen eine Messung der Sprechgeschwindigkeit (tempo) 
# in verschiedenen Kontexten (Faktor Kontext) und in 
# Aeusserungen unterschiedlicher Laenge (Faktor L). 
# Pruefen Sie durch eine Abbildung und statistischen Test, 
# inwiefern die Sprechgeschwindigkeit vom Kontext und 
# von der Aeusserungslaenge beeinflusst wird.

dim(phr)
levels(phr$Kontext)
levels(phr$L)
# Wir haben sehr viele Beobachtungen (naemlich 900), und koennten also 
# ohne Probleme ein MM rechnen. Das Beispiel stammt aber aus f5.R, das zu einer Zeit
# veroeffentlich wurde, als wir MMs noch nicht kannten. Daher mittelten wir
# und berechneten eine ANOVA (siehe unten in Wiederholung aus a5.R). 
# Das ist wohlgemerkt NICHT falsch im eigentlichen Sinne;
# es ist im Text von keinem Random Factor die Rede (d.h. man nimmt nur den Random Facor Vpn an;
# in einer moeglichen Aufgabenstellung in der Klausur waere dieser auch explizit erwaehnt). 
# Die Mittelung war notwendig, da es Wiederholungen gab 
# (und nicht etwa wegen eines nichtbeachteten zusaetzlichen Einflussfaktors)

# Eine Anova hat den "Nachteil", dass man Information wegmittelt; ob diese Information
# relevant ist oder nicht, sei mal dahingestellt; jedenfalls wird die Zufallsvariation durch
# mitteln reduziert; daher haben wir hier pro Fixed Factor NICHT die gleichen F-Werte, 
# wenn wir ANOVA und MM miteinander vergleichen.

# zunaechst die Lösung aus a5.R:
with(phr, table(Vpn, interaction(Kontext, L)))
#in "kurz" gibt es 7 Wdh., in "mittel" 2, in "lang" jedoch nur eine
# Ausserdem: L ist within, Kontext ist within. 
with(phr, table(Vpn, Kontext))
with(phr, table(Vpn, L))
# Aber es gibt, wie gesagt, mehrere Werte pro Stufen-Kombination. Daher mitteln.
phr.m = phr %>%
  group_by(Vpn,Kontext,L) %>%
  summarise(tempo = mean(tempo))
# alles OK
with(phr.m, table(Vpn, interaction(Kontext, L)))
# OK, weil ueberall nur 1 steht
ggplot(phr.m) + 
  aes(y = tempo, x = Kontext, col = L) + 
  geom_boxplot()
#oder gleichwertig
ggplot(phr.m) + 
  aes(y = tempo, x = L, col = Kontext) + 
  geom_boxplot()

# Eventuell kurz < (lang, mittel); ip < paren < pp. 
# Vielleicht keine Interaktion?

# Test
ezANOVA(phr.m, .(tempo), .(Vpn), within = .(Kontext, L))
# Keine signifikante Interaktion. 

# Aber: Sphericity-Korrektur muss beachtet werden! 
# Fuer Kontext nehmen wir HFe da GGe > 0.75. 
# Aber HFe kann ignoriert werden, da HFe > 1. 
# Daher aendern wir nichts an den Freiheitsgraden 
# noch an den Wahrscheinlichkeiten fuer Kontext. 
# (In `Mauchly's Test for Sphericity` war ja auch 
# Kontext derjenige Faktor, die die Sphaeriziaetsannahme
# nicht verletzte).

# Fuer L nehmen wir den GGe-Wert von 0.7434277 und 
# korrigieren dessen Freiheitsgrade damit:
c(2, 58) * .7434277
# [1]  1.486855 43.118807
# Und nehmen die unter p[GG] eingetragene Wahrscheinlichkeit von 
# 2.587909e-22 (die aber auch p < 0.001 ist)

# Schlussfolgerung:
# RT wurde signifikant vom Kontext (F[2,58] = 237.5, p < 0.001) 
# und signifikant von der Sprache (F[1.5, 43.1] = 249.0) beeinflusst, 
# und es gab keine signifikante Interaktion zwischen diesen Faktoren.
# Da bezueglich des Faktors Sprache eine Verletzung der Voraussetzung 
# der Sphaerizitaet vorlag, wurde eine Greenhouse–Geisser Korrektur 
# der Freiheitsgrade vorgenommen.

# Okay, nun das ganze ungemittelt als MM:

# Random Factor Vpn; wie oben gezeigt, sind Kontext und L hierzu "within", daher Formel:
# tempo~Kontext*L + (Kontext+L|Vpn)
phr.lmer = lmer(tempo~Kontext*L + (Kontext+L|Vpn),data=phr)
# Vorsicht: Warnmeldung!
#Warnmeldungen:
#  1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
#    unable to evaluate scaled gradient
#  2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
#    Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
#  3: Model failed to converge with 1 negative eigenvalue: -2.3e-03 
# Seltsamerweise konvergiert allerdings dies hier:
phr.lmer2 = lmer(tempo~Kontext*L + (Kontext|Vpn) + (L|Vpn),data=phr)
# aber zweimal den gleichen Random Factor angeben "soll"darf" man ja nicht, und beide Modelle
# ergeben vollkommen andere Freiheitsgrade, und daher Ergebnisse (die einander zumindest aehneln):
anova(phr.lmer) #z.B. DenDF fuer Kontext = 20.69
anova(phr.lmer2) #z.B. DenDF fuer Kontext = 862.0 (!!!)

# Dennoch mit dem richtigen, aber nicht konvergierenden Modell weitermachen!
phr.lmer = lmer(tempo~Kontext*L + (Kontext+L|Vpn),data=phr)
(phr.step = step(phr.lmer)) # u.a. das Probleme machende Kontext+L 
# faellt aus (Kontext+L|Vpn) heraus;
# nochmal zur Verdeutlichung:
# (Kontext+L|Vpn) geht davon aus, dass Kontext und L jeweils auf eine
# versuchspersonspezifische Weise variieren, weshalb für jede Vpn:Kontext:L-Kombination
# spezifische Intercepts und Steigungen ausgerechnet werden muessen; die Intercepts stehen 
# dann fuer die Unterschiede zwischen den Versuchsperson, die Steigungen dafür, dass diese Vpn
# die Unterschiede zwischen Kontexten und L-Stufen auf UNTERSCHIEDLICHE WEISE realisieren;
# ist letzteres nicht der Fall, koennen die individuellen Steigungen durch eine gemeinsame 
# Steigung ersetzt werden; das ist dann der Fall, wenn alle Versuchspersonen die Unterschiede
# zwischen Kontexten und L-Stufen auf eine SEHR AEHNLICHE WEISE realisieren; 
# dass das Modell hier nicht "konvergiert" hängt vermutlich mit der oben genannten 
# unterschiedlichen Anzahl an Wiederholungen zusammen, siehe:
with(phr, table(Vpn,interaction(Kontext,L))) # mal 7, mal eine, mal 2 Wiederholungen...

#weiter im R-Code:
anova(get_model(phr.step))
#           Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#  Kontext 488.65  244.33     2   866  38.647 < 2.2e-16 ***
#  L       859.08  429.54     2   866  67.944 < 2.2e-16 ***

#vgl. mit dem ANOVA-Ergebnis:
ezANOVA(phr.m, .(tempo), .(Vpn), within = .(Kontext, L))  
#     Effect DFn DFd           F            p p<.05         ges
#2   Kontext   2  58 237.5315852 1.155660e-28     * 0.224680732
#3         L   2  58 248.9863508 3.410839e-29     * 0.331212775
#4 Kontext:L   4 116   0.7008831 5.928824e-01       0.002432093

#wir ignorieren kurz die Sphaerizitaets-Korrektur an dieser Stelle!

# Die Hauptaussage bleibt gleich: Effekte fuer Kontext und L, keine Kontext:L-Interaktion.
# Man beachte aber die voellig unterschiedlichen F-Werte, z.B. fuer Kontext 38.647 vs. 237.5315852
# F = Zufallsvariation / Systematische Variation
# bei MM viel mehr Daten "im Blick", naemlich
nrow(phr) #900
# bei ANOVA "nur"
nrow(phr.m) #270
#--> voellig andere Freiheitsgrade (die von der Anzahl der Beobachtungen abhaengen);
# die Grundaussage bleibt aber erhalten!!!

# -> vgl. im Kontrast dazu das Beispiel 4 in a10.R: gleiche Anzahl an Beobachtungen fuehrt zu identischen
# F-Werten bei ANOVAs vs. MMs

# Wir haben jeweils drei Stufen in Kontext und in L; obschon wir keine Interaktion haben,
# waere also ein zwei post-hoc-Tests wünschenswert (da sowohl Kontext als auch 
# L signifikant waren, wir aber nicht wissen, in jeweils welchen Paarungen)
# -Post-hoc Tests fuer ~L und fuer ~Kontext

#Kein Problem nach dem Mixed Model:

ph.Kontext = pairs(emmeans(get_model(phr.step),~Kontext))
phsel(ph.Kontext)
# ip-paren: p<0.001; ip-pp: p < 0.001, paren-pp: p < 0.001
ph.L = pairs(emmeans(get_model(phr.step),~L))
phsel(ph.L)
#kurz-lang: p < 0.001; kurz-mittel: p < 0.001, lang-mittel: n.s.




#Okay, das geht aber nicht (mit dem dataframe phr.m) mit phoc() nach ezANOVA()

#Wir koennen zwar rechnen
phoc(phr.m, .(tempo), .(Vpn), .(Kontext,L))
# das ergibt aber keinen Sinn, da wir keine Kontext:L-Interaktion festgestellt hatten

# entweder Kontext oder L wegzulassen geht aber nicht, da phoc sonst abstuerzt:
#phoc(phr.m, .(tempo), .(Vpn), .(Kontext)) #stuerzt ab
#phoc(phr.m, .(tempo), .(Vpn), .(L)) #stuerzt ab

# Wir bleiben also dabei (fuer die Klausur):
# phoc() NUR AUSFUEHREN, WENN EINE INTERAKTION ZWISCHEN ZWEI "FIXED" FAKTOREN VORLAG!!!!!
# (entschuldigen Sie bitte mein SCHREIEN!)

# technisch geht es aber eigentlich schon, ist aber sehr umstaendlich, da man jeweils vorher
# erneut mitteln muesste, z.B.:
phr.m.Kontext = phr.m%>%
  group_by(Vpn,Kontext)%>%
  summarise(tempo=mean(tempo))
phoc(phr.m.Kontext, .(tempo), .(Vpn), .(Kontext))
#                   t df     prob-adj
#ip-paren -13.861472 29 7.579419e-14
#ip-pp    -21.361306 29 8.206619e-19
#paren-pp  -8.716232 29 4.057795e-09
# vgl.
phsel(ph.Kontext)


########################
#und fuer L:
phr.m.L = phr.m%>%
  group_by(Vpn,L)%>%
  summarise(tempo=mean(tempo))
phoc(phr.m.L, .(tempo), .(Vpn), .(L))
#                       t df     prob-adj
#kurz-lang   -20.4629379 29 2.646703e-18
#kurz-mittel -25.9397817 29 3.861079e-21
#lang-mittel   0.7776804 29 1.000000e+00
#vgl.
phsel(ph.L)

#Auch hier gilt also: andere Freiheitsgrade, andere t- und p-Werte,
# aber grundsaetzlich bleibt die Grundaussage erhalten!!!

#Ich wiederhole mich sicherheitshalber lieber nochmal:
# Wir bleiben also dabei, dass fuer die Klausur gilt:
# phoc() nach ezANOVA() NUR AUSFUEHREN, WENN EINE INTERAKTION ZWISCHEN ZWEI "FIXED" FAKTOREN VORLAG!!!!!

#pairs(emmeans(get_model(MODEL.lmer))) kann und soll aber auch auf einzelne Faktoren, 
# die mehr als zwei Stufen aufweisen, angewendet werden!