# Episode IV: Nun geht's wirklich los


# Für die Episoden IV-VI nutzen wir, zusätzlich zum tidyverse, auch das tidycomm-Package.
# Also installieren und dann hier entsprechend laden.
library(tidyverse)
library(tidycomm)


# Was ist das?
vignette('tidycomm')


# Hmm, keine Vignette? Und wenn wir online einfach mal "tidycomm" suchen ...?
# Aha, mehrere Vignetten!
browseVignettes('tidycomm')


# Und: tidycomm hat zwei Datensätze mit an Bord, die wir nutzen können.
# Hier nutzen wir "WoJ", was eine Stichprobe aus der Worlds-of-Journalism-Befragungen enthält.
# Darin wurden in mehreren Wellen 2012-2016 JournalistInnen weltweit zu Ihrer Arbeit befragt.
# Wir speichern die Daten zunächst in eine eigene Variable, um damit weiter zu arbeiten.
woj_befragung <- WoJ
str(woj_befragung)
woj_befragung



# 1200 Befragungen also, mit JournalistInnen aus aus 5 Ländern.
# Um hier einen deskriptiven Überblick zu bekommen, bietet tidycomm quasi Soforthilfe.
woj_befragung %>% 
  describe()



# Was ist passiert?
# Für jede numerische (!) Variable erhalten wir Informationen, über ...
# - die Anzahl der Fälle
# - fehlende Werte
# - Mittelwert (M) und Standardabweichung (SD)
# - Mindestwert (Min) und Maximalwert (Max)
# - die Spannweite (Range), also max-min
# - das 25er- (Q25), 50er- (Median, Mdn) und 75er- (Q75) Perzentil (auch: Quantil)
# - der Median (Mdn)
# - die Schiefe der Verteilung (Skew bzw. Skewness)
# - die Wölbung der Verteilung (Kurtosis)



# Natürlich wissen wir längst, wie wir solche Daten auch visualisieren können.
# Zum Beispiel als Histogramm.
# (Zur Erinnerung: Die Daten sind nur ein [verzerrter] Ausschnitt der gesamten Befragung!)
# Lt. Beschreibung ist "autonomy_selection" übrigens die Antwort auf folgende Frage:
#   "Thinking of your work overall, how much freedom do you personally have in selecting news stories you work on?"
#   (1 = "no freedom at all" bis 5 = "complete freedom")
woj_befragung %>% 
  ggplot(aes(x = autonomy_selection)) +
  geom_histogram()



# Was aber ist mit den nicht-numerischen Variablen?
woj_befragung %>% 
  ggplot(aes(x = employment)) +
  geom_histogram()

woj_befragung %>% 
  ggplot(aes(x = employment)) +
  geom_histogram(stat = 'count')

woj_befragung %>% 
  ggplot(aes(x = country)) +
  geom_histogram(stat = 'count')



# Und wenn wir das kombinieren wollen, also autonomy_selection je Land?
woj_befragung %>% 
  group_by(country) %>% 
  describe()

woj_befragung %>% 
  ggplot(aes(x = country,
             y = autonomy_selection)) +
  geom_bar(stat = 'identity')


# Übrigens kann ggplot auch fancy "violin plots", die mehr Streuungsmaße darstellen.
# Solche Visualisierungen lesen sich wie aufgestellte Histogramme.
woj_befragung %>% 
  ggplot(aes(x = country,
             y = autonomy_selection)) +
  geom_violin()

# Der Einfachheit halber können wir sie auch "hinlegen".
woj_befragung %>% 
  ggplot(aes(x = country,
             y = autonomy_selection)) +
  geom_violin() +
  coord_flip()

# Es zeigt sich (wie schon in der describe()-Funktion vorhin), dass ...
# - die meisten Nennungen in allen Ländern auf die "4" entfallen
# - die Daten linksschief (rechtssteil) sind
# - es kaum (in Österreich: keine) 1-Nennungen gibt


# Und noch eine fancy Darstellung der Häufigkeit je Gruppe: Boxplots
# Sie funktionieren ähnlich wie "violin plots", zeigen aber ...
# - als Box den Interquartilsabstand (IQR) vom 25er- zum 75-Quantil (die "mittleren" 50% der Stichprobe)
# - als waagrechte dicke Linie in der Box den Median (Q50)
# - als senkrechte Box-Verlängerungs-Linien ("Whisker") das etwa 1.5-fache des IQR
# - als Punkte die sonstigen Ausreißer
woj_befragung %>% 
  ggplot(aes(x = country,
             y = autonomy_selection)) +
  geom_boxplot()


# Histogramme, Boxplots und "violin plots" helfen dabei, seine Daten kennenzulernen.
# Für die finale Version einer wissenschaftlichen Arbeit eignen sie sich aber nur in Ausnahmefällen.



# Gerade bei Befragungsdaten interessiert zunächst häufiger die reine Verteilung der Antwortoptionen.
# Man könnte auch sagen: Wie zwei kategoriale Variablen zusammenspielen.
# Dafür nutzen wir Kreuztabellen, die wir, das kennen wir ja schon, mit group_by und count anlegen können.
# Da die Daten dann im Long-Format vorliegen, müssen wir sie mit pivot_wider verbreitern.
# Als Beispiel schauen wir uns die Anstellungsverhältnisse über die Ländern an.
woj_befragung %>% 
  group_by(country) %>% 
  count(employment)

?pivot_wider

woj_befragung %>% 
  group_by(country) %>% 
  count(employment) %>% 
  pivot_wider(names_from = 'employment',
              values_from = 'n')



# Mit tidycomm geht es auch etwas "schöner".
?crosstab
woj_befragung %>% 
  crosstab(employment, country)


# Was fehlt, ist der statistische Test dazu, ein Chi²-Test.
# Den gibt es standardmäßig bei R mit an Bord (chisq.test).
# Aber aufmerksame LeserInnen haben natürlich auch den Parameter "chi_square" in crosstab entdeckt.
?crosstab
woj_befragung %>% 
  crosstab(employment, country,
           chi_square = TRUE)

# Here we go:
# Chi² = 158.16 bei 8 Freiheitsgraden und einem p < .001 (***)



# Zur Visualisierung eignen sich für solche kategorialen Zählungen übrigens Balkendiagramme. 
# Die lassen sich (für den Vergleich) nebeneinander oder (für Verteilungen) gestapelt anordnen.
woj_befragung %>% 
  group_by(country) %>% 
  count(employment) %>% 
  ggplot(aes(x = employment,
             y = n)) +
  geom_bar(stat = 'identity')


# Anstelle der gesamten Stichprobe können wir das auf Länder aufteilen.
# "facet_wrap" haben wir dafür ja schon kennengelernt.
# Eine Alternative dazu sind gestapelte Balken, die mit unterschiedlichen Farben gefüllt sind.
# Alles, was ggplot dafür braucht, ist ein Hinweis auf die Füllfarbe (fill) als gruppierende Variable.
woj_befragung %>% 
  group_by(country) %>% 
  count(employment) %>% 
  ggplot(aes(x = country,
             y = n,
             fill = employment)) +
  geom_bar(stat = 'identity')


# Schöner im Liegen?
# Dann kippen wir das Ganze eben.
# Und ergänzen in dem Zuge auch besser lesbare Farben, indem wir die Farbe wie eine Achse ansprechen.
woj_befragung %>% 
  group_by(country) %>% 
  count(employment) %>% 
  ggplot(aes(x = country,
             y = n,
             fill = employment)) +
  geom_bar(stat = 'identity') +
  coord_flip() +
  scale_fill_brewer(type = 'qual')



# Wollten wir, der Vollständigkeit halber, keine gestapelten Balken, dann ...
# Nun ja, Hilfefunktion und Suchmaschinen und Cheatsheets sind ja mittlerweile mehr als bekannt.
woj_befragung %>% 
  group_by(country) %>% 
  count(employment) %>% 
  ggplot(aes(x = country,
             y = n,
             fill = employment)) +
  geom_bar(stat = 'identity',
           position = 'dodge') +
  coord_flip() +
  scale_fill_brewer(type = 'qual')



# Damit sind alle Fragen in woj_befragung beschreibbar.
# Alle Fragen?
# Nicht ganz. Eine kleine Gruppe an Variablen sticht aus der Masse heraus.
# Denn die Variablen ethics_1 bis ethics_4 sind eigentlich vier Items eines Faktors.
# Sie sagen was darüber aus, wie stark Journalismus ethischen Prinzipien verpflichtet sein sollte.
# Dafür bietet tidycomm eine sehr einfache Funktion:
?add_index


# Für einen einfachen Mittelwertindex lautet das Skript also lediglich wiefolgt:
woj_befragung %>% 
  add_index(ethik_mittelwertindex, 
            ethics_1, 
            ethics_2, 
            ehtics_3, 
            ethics_4)

# Da ist nicht allzuviel sichtbar, aber wir erinnern uns, dass wir Tibbles auch anzeigen können:
woj_befragung %>% 
  add_index(ethik_mittelwertindex, 
            ethics_1, 
            ethics_2, 
            ethics_3, 
            ethics_4) %>% 
  View()


# Hervorragend!
# Doch halt: Laden die vier Items überhaupt auf einen Faktor, den wir hier als Index berechnet haben?
# Dafür nutzen wir in der Regel Cronbachs Alpha. 
# Und auch dafür bietet tidycomm eine sehr einfache Funktion:
?get_reliability

woj_befragung %>% 
  get_reliability()

woj_befragung %>% 
  add_index(ethik_mittelwertindex, ethics_1, ethics_2, ethics_3, ethics_4) %>% 
  get_reliability()

# "Geht so" könnte man sagen. Also das Ergebnis.
# Die methodische Umsetzung hingegen in R an dieser Stelle ist grandios einfach!



# Dieser (metrische) Index lässt sich vielleicht mit der (quasi-metrischen) Arbeitserfahrung erklären?
# Dafür korrelieren wir die beiden.
# Übrigens lassen sich auch Korrelationsbeziehungen über Visualisierungen "kennenlernen".
woj_befragung %>% 
  add_index(ethik_mittelwertindex, ethics_1, ethics_2, ethics_3, ethics_4) %>% 
  ggplot(aes(x = ethik_mittelwertindex,
             y = work_experience)) +
  geom_point()


# Unübersichtlich?
# Mithilfe von geom_smooth lässt sich eine angenäherte Trendgerade einbauen:
# - zeigt die abwärts, können wir eine negative Korrelation erwarten
# - zeigt die aufwärts, wird die Korrelation wohl positiv ausfallen
# - je horizontaler, desto schwächer die Korrelation
# - je steiler, desto stärker die Korrelation
woj_befragung %>% 
  add_index(ethik_mittelwertindex, ethics_1, ethics_2, ethics_3, ethics_4) %>% 
  ggplot(aes(x = ethik_mittelwertindex,
             y = work_experience)) +
  geom_point() +
  geom_smooth(method = 'lm',
              se = FALSE)


# Für eine Veröffentlichung eignet sich diese Form der Visualisierung in der Regel nicht.
# Hier setzen wir vielmehr auf die bekannten Korrelationskoeffizienten.
# Die berechnet R zwar auch von Haus aus (?cor), aber das ist sehr unhandlich.
# Viel praktischer an der Stelle ist mal wieder tidycomm.
?correlate
woj_befragung %>% 
  add_index(ethik_mittelwertindex, ethics_1, ethics_2, ethics_3, ethics_4) %>% 
  correlate(ethik_mittelwertindex, work_experience)

# Negativ, signifikant, aber ein seeehr schwacher Zusammenhang.