# 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
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))
# alles OK
with(phr.m, table(Vpn, interaction(Kontext, L)))

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

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

# $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 da Interaktion
p = phoc(dbc, .(d), .(Vpn), .(Eink, Dialekt))
# 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). 

# 4. Fuer den Data-Frame
rating = read.table(file.path(pfadu, "r.m.txt"))
# Inwiefern wird Rating von der Grammatikalitaet (Gram) und Bekanntheit (Fam) beeinflusst?

# beide Faktoren sind within
with(rating, table(Vpn, interaction(Gram, Fam)))

# Bild
ggplot(rating) + 
  aes(y = Rating, x = Gram, col = Fam) + 
  geom_boxplot()

# 1. Gram: (High > Moderate)
# 2. Fam: Old > New aber eventuell nur in 'Moderate'
# 3. daher wahrscheinlich eine Interaktion, da der Abstand zwischen High
# und Moderate groesser ist fuer New

ezANOVA(rating, .(Rating), .(Vpn), .(Gram, Fam))
#Effect DFn DFd          F            p p<.05        ges
#2     Gram   1  25 241.352814 2.368171e-14     * 0.52778958
#3      Fam   1  25  13.504271 1.136189e-03     * 0.07594048
#4 Gram:Fam   1  25   8.850306 6.412278e-03     * 0.03578053

# Rating wurde signifikant von der Grammatikalitaet beeinflusst
# (F[1,25] = 241.4, p < 0.001), und von Fam (F[1, 25] = 13.5, p < 0.01) 
# und es gab eine signifikante Interaktion zwischen diesen Faktoren (F[1, 25] = 8.9, p < 0.001).
ph = phoc(rating, .(Rating), .(Vpn), .(Gram, Fam))
round(phsel(ph$res), 3)
round(phsel(ph$res, 2), 3)

#                             t   df prob-adj
#High:New-Moderate:New       10.844 25     0
#High:Old-Moderate:Old     13.090 25        0
#High:New-High:Old         -1.251 25    1.000
#Moderate:New-Moderate:Old -3.926 25    0.004

# Post-hoc Bonferroni t-tests zeigten signifikante Unterschiede 
# zwischen High und Moderate sowohl fuer New (p < 0.001) 
# als auch fuer Old (p < 0.001); 
# es gab auch signifikante Unterschiede zwichen New und Old fuer Moderate 
# (p < 0.01), jedoch nicht fuer High.


# 5. Diese Tabelle
# http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf/stable.pdf
# aus Sussman et al (1997) zeigt sogenannte Lokus-Steigungen fuer 5 Sprecher (M# bis M#5) und 5 Sprecherinnen (F#1 bis F#5). Die Lokus-Steigungen sind in der Spalte unter `k` und sie kommen vor in silbeninitialer, silbenmedialer, und silbenfinaler Position (daher 10 k-Eintraege  pro Position; 3 k-Eintrage pro Sprecher oder Sprecherin).Inwiefern wird `k` vom Geschlecht und/oder der Silbenposition beeinflusst?

# Data-Frame bauen
werteinit = c(.75, .74, .82, .75, .61, .71, .88, .78, .84, .77)
wertemed = c(.79, .81, .79, .68, .69, .74, .81, .77, .84, .73)
wertefin = c(.68, .74, .62, .52, .45, .26, .34, .49, .58, .24)
werte = c(werteinit, wertemed, wertefin)
posn = c(rep("i", 10), rep("m", 10), rep("f", 10))
g = c(rep("M", 5), rep("W", 5), rep("M", 5), rep("W", 5), rep("M", 5), rep("W", 5))
vpn = rep(c(paste("M", 1:5, sep=""), paste("F", 1:5, sep="")), 3)

lok.df = data.frame(werte, P = posn, G = g, Vpn = vpn)
head(lok.df); dim(lok.df)

with(lok.df, table(Vpn, interaction(P, G)))
# Position within, Geschlecht between

lok.df$P = factor(lok.df$P,levels=c("i","m","f"))

ggplot(lok.df) + 
  aes(y = werte, x = P, col = G) + 
  geom_boxplot()


# Position: 
ggplot(lok.df) + 
  aes(y = werte, x = P) + 
  geom_boxplot()
# kaum Unterschiede zwischen initial und medial;
# final ist tiefer als die anderen beiden
# Geschlecht: wahrscheinlich kaum Unterschiede, wie man hier sieht:
ggplot(lok.df) + 
  aes(y = werte, x = G) + 
  geom_boxplot()

ggplot(lok.df) + 
  aes(y = werte, x = G, col = P) + 
  geom_boxplot()
# es wird sicherlich eine Interaktion geben,
# da der Abstand final vs. initial/medial groesser
# ist fuer Frauen als fuer Maenner.

ezANOVA(lok.df, .(werte), .(Vpn), .(P), between = .(G))
# Hier muessen die Freiheitsgrade wegen Sphericity geaendert werden
# Greenhouse-Geisser korrigierte Freiheitsgrade fuer P
round (c(2, 16) * 0.6553049, 1)
# fuer G:P
round (c(2, 16) * 0.6553049, 1)

# Die Steigungen wurden signifikant von der Position 
# (F[1.3, 10.5] = 60.0, p < 0.001) jedoch nicht vom Geschlecht beeinflusst 
# und es gab eine signifikante Interaktion zwischen diesen
# Faktoren (F[1.3, 10.5] = 14.2, p < 0.001).

# Post-hoc tests
p = phoc(lok.df, .(werte), .(Vpn), .(P,G))
round(phsel(p$res), 3)
#            t df prob-adj
#i:M-m:M -0.616  4    1.000
#i:M-f:M  3.099  4    0.544
#i:W-m:W  1.050  4    1.000
#i:W-f:W  7.012  4    0.033
#m:M-f:M  5.207  4    0.097
#m:W-f:W  7.669  4    0.023

round(phsel(p$res, 2), 3)
#           t    df prob-adj
#i:M-i:W -1.375 7.829    1.000
#m:M-m:W -0.751 7.436    1.000
#f:M-f:W  2.602 7.609    0.493

# Post-hoc tests zeigten einen signifikanten Einfluss
# zwischen initialen und finalen Steigungen in Frauen (p < 0.05) 
# und zwischen medialen und finalen Steigungen in Frauen (p < 0.05).
# Die Position hatte keinen signifikanten Einfluss auf die
# Steigungen fuer Maenner.
# Die Tests bestaetigten, dass die Steigungen
# nicht vom Geschlecht beeinflusst wurden.