library(ggplot2)
library(dplyr)
library(ez)
source(file.path(pfadu, "phoc.txt"))
blang  = read.table(file.path(pfadu, "blang.txt"))
v.df = read.table(file.path(pfadu,"vokal.txt"))
dg = read.table(file.path(pfadu, "dg.txt"))
ssb = read.table(file.path(pfadu, "ssb.txt"))

##########################################################
# 1. Zwei Faktoren 
######################################################################
head(dg); dim(dg)
#
# Wenn nur ein 1 Wert pro Reihe belegt ist,
# dann sind beide Faktoren 'between'
with(dg, table(Vpn, interaction(Region, Gen)))

# 2 Faktoren, 3 Fragen:
# Wird F2 von Gender beeinflusst? Von der Region?
ggplot(dg) + 
  aes(y = F2, x = Gen, colour = Region) + 
  geom_boxplot()

##########################################################
# 2. Interaktionen 
#############################################################

# Liegt eine Interaktion zwischen Gender und Region
# vor? Bedeutung: sind die Unterschiede zwischen A, B, C
# ähnlich für männlich und weiblich? (Oder umgekehrt:
# sind die Unterschiede zwischen männlich und weiblich
# ähnlich in den Regionen A, B, C?). Wenn nein, dann
# liegt eine Interaktion vor. 
# Man kann teilweise eine Interaktion auch graphisch darstellen,
# indem die Medianwerte von dem Boxplot miteinander 
# verbunden werden:
# Hier sind die Medianwerte pro Stufe
dg.m = dg %>%
  group_by(Gen,Region) %>%
  summarise(F2 = median(F2))
# Hier werden die Medianwerte getrennt für männlich
# und weiblich verbunden
ggplot(dg.m) + 
  aes(y = F2, x = Region, group=Gen, colour = Gen) +  
  geom_line()


# Wenn die Linien mehr oder weniger parallel zueinander
# verlaufen, liegt keine Interaktion vor. Je mehr sie davon
# abweichen, umso wahrscheinlicher wird die
# Interaktion signifikant sein.
ezANOVA(dg, .(F2), .(Vpn), between =.(Region, Gen))
# F2 wurde von der Region (F[2,54] = 119.6, p < 0.001)
# und vom Geschlecht (F[1,54] = 106.2, p < 0.001) 
# signifikant beeinflusst und es gab eine signifikante
# Interaktion (F[2,54] = 12.1, p < 0.001) zwischen
# diesen beiden Faktoren. Post-hoc Bonferroni-korrigierte
# t-Tests zeigten...

##########################################################
# 3. post-hoc t-tests
#############################################################

## (Wenn eine signifikante Interaktion vorliegt).
## Post-hoc t-Tests vergleichen alle möglichen Stufen-Paare
## d.h. es gibt A, B, C sowie m, w daher Kombinationen von
## A-m mit A-w, A-m mit B-m, A-m mit B-w...
## Es gibt 15 mögliche Paare
## (errechnet sich as 6!/(4!2!) = 30/2 = 15)
## 6 weil es 3 * 2 = 6 Stufen gibt
## 6! bedeutet: 6 * 5 * 4 * 3 * 2 * 1. Es gibt also 15 Tests:
factorial(6)/(factorial(4) * factorial (2))

phoc(dg, .(F2), .(Vpn), .(Region, Gen))
## Bonferroni-Korrektur:
## Zusätzlich werden die Wahrscheinlichkeiten
## mit der Anzahl der Tests multipliziert
## (d.h. es ist 15 Mal so schwierig, Signifikanzen
## zu bekommen). z.B. ein Wert von 0.04
## ist nach dieser Berechnung 0.04 * 15 = 0.6
## daher nicht signifikant

## Auswahl
vok.ph = phoc(dg, .(F2), .(Vpn), .(Region, Gen))
## Wir benötigen aber nicht alle 15 Tests
## sondern nur diejenigen wenn ein Faktor variiert
## und die anderen konstant sind, daher:
# Faktor 1 variiert, Faktor 2 ist konstant
# Bedeutung: gibt es Unterschiede zwischen den
# Dialekten (Faktor 1 variiert) getrennt für Männer
# und Frauen?
p1 = phsel(vok.ph$res, 1)
round(p1, 3)
# Faktor 2 variiert, Faktor 1 ist konstant
# Bedeutung: gib es Unterschiede zwischen männlich/weiblich
# (Faktor 2 variiert) getrennt in den Dialekten?
p2 = phsel(vok.ph$res, 2)
round(p2, 3)
# Post-hoc Bonferroni-korrigierte
# t-Tests zeigten signifikante Unterschiede
# zwischen den Dialekten A und C
# (Für Männer: p < 0.001; Für Frauen p < 0.001)
# sowie zwischen
# den Dialekten B und C 
# (Für Männer: p < 0.001; Für Frauen p < 0.001)
# Der Unterschied
# zwischen den Dialekten A und B war nicht signifikant,
# weder in Männern noch in Frauen. 
# Schließlich gab es signifikante Unterschiede
# zwischen Männern und Frauen in Dialekt A (p < 0.001)
# und in Dialekt B (p < 0.001) aber nicht in Dialekt C.

##########################################################
# 5. Wiederholungen in within-Stufen
#############################################################

#
# Tabelle: wieviele Stichproben pro Vpn?
with(ssb, table(Vpn, interaction(Wort, Alter)))

# Über die Wiederholungen mitteln
ssbm = ssb %>%
  group_by(Wort,Alter,Vpn) %>%
  summarise(F2 = mean(F2))
head(ssbm)
dim(ssbm)
# Tabelle: wieviele Stichproben pro Vpn?
with(ssbm, table(Vpn, interaction(Wort, Alter)))

# Abbildung
ggplot(ssbm) + 
  aes(y = F2, x = Alter, colour = Wort) + 
  geom_boxplot()

# Statistik
ezANOVA(ssbm, .(F2), .(Vpn), .(Wort), between = .(Alter))

##########################################################
# 6. Greenhouse-Geisser-Korrektur
#############################################################

# Korrektur der Freiheitsgrade
2 * 0.6860511
20 * 0.6860511

# Aufgerundet
round(2 * 0.6860511, 1)
round(20 * 0.6860511, 1)


1 - pf(9.8908882, 2 * 0.6860511, 20 * 0.6860511)