# Zwei-Faktoren-Varianzanalyse

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

# 1.
auf = read.table(file.path(pfadu, "auf.txt"))
# Die Daten zeigen Reaktionszeiten auf schwedische Woerter 
# von franzoesischen und englischen Versuchspersonen (Faktor Lang) 
# nach einem 0 oder 6 monatigen Aufenthalt (Faktor Monat) in Schweden. 
# Werden die Reaktionszeiten von der Sprache und/oder Aufenthaltsdauer beeinflusst?

# Welche Faktoren sind within oder between
with(auf, table(Vpn, interaction(Lang, Monat)))
# Lang ist between, Monat ist within
# Bild
ggplot(auf) + 
  aes(x = RT, col = Monat) + 
  geom_density() + 
  xlab("Reaktionszeit") + 
  ylab ("Wahrscheinlichkeitsdichte") + 
  facet_wrap(~Lang)
# am besten:
ggplot(auf) + 
  aes(y = RT, x = Monat, col = Lang) + 
  geom_boxplot() + 
  xlab("Monat") + 
  ylab ("Reaktionszeit (ms)")


# 1. wird RT von Lang beeinflusst?
# 2. wird RT von Sprache beeinflusst?
# 3. Gibt es zwischen Lang und Sprache eine Interaktion?
# Ist der Unterschied zwischen 0 und 6 Monaten anders fuer die 
# Englaender als fuer die Franzosen? 
# Wenn ja, dann gibt es wahrscheinlich eine Interaktion. Hier wahrscheinlich ja...
# Test
ezANOVA(auf, .(RT), .(Vpn), within = .(Monat), between = .(Lang))
#Effect DFn DFd        F           p p<.05
# 2       Lang   1   8 19.75207 0.002155491     *
# 3     Monat    1   8 11.05323 0.010467542     *
# 4 Lang:Monat   1   8 12.87092 0.007109191     *


# post-hoc t-test
p = phoc(auf, .(RT), .(Vpn), .(Monat, Lang))
# Die Ergebnisse fuer Faktor 1, gerundet auf drei Nachkommastellen:
round(phsel(p$res, 1), 3)
# Die Ergebnisse fuer Faktor 2:
round(phsel(p$res, 2), 3)

# Eine ANOVA mit abhaengiger Variable Reaktionszeit und 
# mit between-Faktor Sprache (2 Stufen: Englisch/Franzoesisch) 
# und within-Faktor Aufenthaltsdauer (2 Stufen: 0 oder 6 Monate) 
# zeigte einen signifikanten Einfluss von der Sprache 
# (F[1, 8] = 19.8, p < 0.01) und von der Aufenthaltsdauer 
# (F[1,8] = 11.1, p < 0.05) und auch eine signifikante Interaktion 
# zwischen diesen Faktoren (F[1, 8] = 12.9, p < 0.01). 
# Post-hoc Bonferroni-korrigierte t-Tests zeigten signifikante 
# Unterschiede zwischen Englaedern und Franzosen in der Reaktionszeit 
# bei 0 Monaten (p < 0.05); sonst gab es keine anderen signifikanten Unterschiede.

# Dies ist uebrigens ein interessanter Fall, wie die Bonferroni-Korrektur
# dazu fuehrt, dass ziemlich deutliche Unterschiede ploetzlich nicht mehr als 
# statistisch signifikante Unterschiede berichtet werden duerfen. Betrachten wir
# den Boxplot, so stellen wir fest, dass der Unterschied bei den Franzosen
# zwischen 0m und 6m sehr deutlich zu sein scheint.
# Der post-hoc-durchgefuehrte gepaarte t-Test hat aber ein Bonferroni-adjustiertes 
# p von 0.148. Ohne die Bonferroni-Korrektur waere das
0.148/6
# = 0.02466667
# D.h. haetten wir nur Franzosen, haetten wir einen gepaarten t-Test durchfuehren koennen,
# der ein statistisch signifikantes Ergebnis erbracht haette 
# (wie die Abbildung ja auch zeigt). Die Bonferroni-Korrektur
# ist aber auch die "konservativste" Korrektur-Methode ("konservativ" im Sinne von
# "vorsichtig damit, Signifikanzen zu verkuenden").

# 2.
phr = read.table(file.path(pfadu, "phr.df.txt"))
dim(phr)
levels(phr$Kontext)
# 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.

with(phr, table(Vpn, interaction(Kontext, L)))
# L ist within, Kontext ist within. 
# Aber es gibt mehrere Werte pro Stufen-Kombination. Daher mitteln.
phr.m = phr %>%
  group_by(Vpn,Kontext,L) %>%
  summarise(tempo = mean(tempo))
# oder 
phr.m2 = aggregate(tempo~Vpn+Kontext+L,FUN = mean, data = phr)
# alles OK
with(phr.m, table(Vpn, interaction(Kontext, L)))
with(phr.m2, table(Vpn, interaction(Kontext, L)))

ggplot(phr.m) + 
  aes(y = tempo, x = Kontext, col = L) + 
  geom_boxplot()

# oder

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.

# 3. Fuer den Data-Frame

dbc = read.table(file.path(pfadu, "dbc.txt"))
# inwiefern wird wird die Dauer (d) vom Dialekt und/oder Einkommen beeinflusst?
# Beide Faktoren sind between
with(dbc, table(Vpn, interaction(Dialekt, Eink)))
# 0der zwei tables durch
with(dbc,table(Vpn, Dialekt))
with(dbc,table(Vpn, Eink))
# Bild
ggplot(dbc) + 
  aes(y = d, x = Eink, col = Dialekt) + 
  geom_boxplot()
# Frage 1: wird d von Einkommen beeinflusst? Ja:
ggplot(dbc) + 
  aes(y = d, x = Eink) + 
  geom_boxplot()
# Frage 2: wird d vom Dialekt beinflusst? Koennte sein...
ggplot(dbc) + 
  aes(y = d, x = Dialekt) + 
  geom_boxplot()
# Frage 3: Interaktion? Eindeutig ja.
ggplot(dbc) + 
  aes(y = d, x = Eink, col = Dialekt) + 
  geom_boxplot()
# Der Unterschied zwischen N und S ist nicht derselbe in High vs Low

ezANOVA(dbc, 
        .(d), 
        .(Vpn), 
        between = .(Eink, Dialekt))
#zunächst den LEvene-Test beachten:
# $`Levene's Test for Homogeneity of Variance`
# DFn DFd   SSn     SSd         F         p   p<.05
# 1   3  36 37.7 2766.9 0.1635043 0.9202332      

#ist aber vollkommen in Ordnung (da p = 0.9202332)!

#also:

# $ANOVA
# Effect DFn DFd         F           p p<.05       ges
# 1         Eink   1  36 12.028921 0.001375342     * 0.2504516
# 2      Dialekt   1  36  4.208065 0.047562566     * 0.1046572
# 3 Eink:Dialekt   1  36  9.273933 0.004329839     * 0.2048405

# post-hoc-t-Tests, da Interaktion
p = phoc(dbc, .(d), .(Vpn), .(Eink, Dialekt)) #OHNE "between =" !!!
# Faktor 1:
round(phsel(p$res, 1), 3)
# Faktor 2:
round(phsel(p$res, 2), 3)

# Schlussfolgerung
# Die Dauer wurde signifikant vom Einkommen (F[1, 36] = 12.0, p < 0.01) 
# und vom Dialekt (F[1, 36] = 4.2, p < 0.05) beeinflusst 
# und es gab eine signifikante Interaktion zwischen dieses Faktoren 
# (F[1,36] = 9.3, p < 0.001).  
# Post-hoc Bonferroni korrigierte t-tests zeigten einen signifikanten 
# Unterschied zwischen Nord und Sueddialekten fuer das hohe Einkommen 
# (p < 0.05), sowie einen signifikanten Unterschied zwischen 
# Einkommensgruppen im Nord-Dialekt (p < 0.01).