# Episode VI: ExperimentieRen

# Wie immer beginnen wir mit dem Laden der benötigten Pakete.
library(tidyverse)
library(tidycomm)



# Zunächst benötigen wir Daten.
# Dafür laden wir Teile der Daten einer Experimentalbefragung zu automatisiertem Journalismus (AJ).
# Probanden sahen einen entweder von einem Journalisten (J) oder einem Algorithmus (A) generierten Text.
# Im Beitrag ging es um Finanzentwicklungen des DAX.
# Sie gaben dann auf 5-stufen Skalen an, ob der Text objektiv, lebhaft und gut geschrieben ist.
aj <- read_csv('haim-graefe_automated-news_abs-partly.csv')



# Im ersten Schritt wollen wir einen t-Test rechnen.
# Wir vergleichen, ob Personen, die einen J-Text sahen, den Text anders bewerten, als Personen mit A-Text.
# Wir filtern und berechnen Mittelwerte.
aj %>% 
  group_by(text, item) %>% 
  summarise(mittelwert = mean(einschaetzung))

# NA?
?mean

# Okay, also NA values vorher strippen, wie die Dokumentation sagt ...
aj %>% 
  group_by(text, item) %>% 
  summarise(mittelwert = mean(einschaetzung,
                              na.rm = TRUE))


# Spannend!
# Und das kriegen wir auch locker visuell hin.
aj %>% 
  group_by(text, item) %>% 
  summarise(mittelwert = mean(einschaetzung,
                              na.rm = TRUE)) %>%  
  ggplot(aes(x = text,
             y = mittelwert)) +
  geom_bar(stat = 'identity') +
  facet_wrap(vars(item))



# Okay, jetzt der t-Test.
# Der ist dank tidycomm sehr einfach umsetzbar.
?t_test


# Allerdings braucht der an der Stelle ein "breites" Format.
# Kein Problem, umformen können wir ja längst.
?pivot_wider
aj %>% 
  pivot_wider(names_from = item,
              values_from = einschaetzung)


# Jetzt aber: t-Test
# Wichtig ist dabei der paired-Parameter, der angibt, ob es sich um gepaarte Stichproben handelt.
aj %>% 
  pivot_wider(names_from = item,
              values_from = einschaetzung) %>% 
  t_test(group_var = text,
         objektiv, lebendig, gut_geschrieben,
         paired = FALSE)

# Ergebnisse:
# -  Bewertungen für "lebendig" und "gut geschrieben" unterscheiden sich signifikant voneinander
# - insbesondere bei "lebendig" liegt fast ein Skalenpunkt auf der 5-stufigen Skala zwischen A und J





# Für die einfaktorielle Varianzanalyse schauen wir uns an, welchen Einfluss das Geschlecht hat.
# Das Geschlecht ist auf drei Stufen codiert, m, w und d.
# Wir wollen die AnOVa nur für "lebendig" als abhängige Variable berechnen.
# Damit entfällt das Verbreitern der Daten, wenn wir stattdessen nur auf "lebendig" filtern.
?unianova
aj %>% 
  filter(item == 'lebendig') %>% 
  unianova(group_var = geschlecht,
           einschaetzung,
           descriptives = TRUE)


# Nicht signifikant. 

# Wir schauen uns nichtsdestotrotz noch den Post-hoc-Test an.
# Geschulte StatistikerInnen mögen jetzt anmerken, dass das VERBOTEN ist, weil ja nicht signifikant.
# Trotzdem ... tidycomm kennt Tukeys HSD, der sich an dieser Stelle zuschalten lässt.
aj %>% 
  filter(item == 'lebendig') %>% 
  unianova(group_var = geschlecht,
           einschaetzung,
           descriptives = TRUE,
           post_hoc = TRUE)


# Fieserweise verschachtelt tidycomm den Post-hoc-Test in ein weiteres tibble.
# Das rauzufieseln, ist aber kein Problem.
# Und "rausfieseln" heißt im tidyverse pull.
?pull

aj %>% 
  filter(item == 'lebendig') %>% 
  unianova(group_var = geschlecht,
           einschaetzung,
           descriptives = TRUE,
           post_hoc = TRUE) %>% 
  pull(post_hoc)

# Wieder nicht signifikant.
# (Geschulte StatistikerInnen: Habichdochgesagt!)

  
  

# Können wir jetzt A/J und das Geschlecht kombinieren? Also, in so einer mehrfaktoriellen AnOVa?
# Jup, klar. Das heißt dann übrigens MAnOVa.
# Erstmal visuell.
aj %>% 
  filter(item == 'lebendig') %>% 
  group_by(text, geschlecht) %>% 
  summarise(mittelwert = mean(einschaetzung, na.rm=T)) %>% 
  ggplot(aes(x = text,
             y = mittelwert,
             color = geschlecht)) +
  geom_line(stat = 'identity')


# Huch, was ist denn da falsch?
# *liest, recherchiert, stellt fest, dass ggplot nicht weiß, welche Punkte verbunden gehören*
aj %>% 
  filter(item == 'lebendig') %>% 
  group_by(text, geschlecht) %>% 
  summarise(mittelwert = mean(einschaetzung, na.rm=T)) %>% 
  ggplot(aes(x = text,
             y = mittelwert,
             color = geschlecht,
             group = geschlecht)) +
  geom_line(stat = 'identity')



# Das sieht vertraut aus!

# Die statistische MAnOVa-Auswertung gibt es die derzeit nicht in so schön.
# Wir müssen dafür zunächst die Daten in ein sauberes separates tibble speichern.
aj_manova_daten <-
  aj %>% 
  filter(item == 'lebendig')

# Zur Berechnung arbeiten wir mit der R-Funktion lm und dem sjstats-Paket (--> installieren!).
# - lm (Linear Model) berechnet ein lineares Modell, wie bei einer Regression
# - sjstats > anova_stats extrahiert daraus die nötigen Parameter für eine Varianzanalyse
?lm

library(sjstats)
vignette('sjstats')
browseVignettes('sjstats')
?anova_stats


# Für lm müssen wir die MAnOVa als Formel im R-Stil beschreiben.
# Das klingt kompliziert, ist aber eigentlich sehr intuitiv:
# - abhängige Variable (AV) kommt auf die eine Seite einer Tilde (~)
# - unabhängige Variablen (UVs) auf die andere Seite der Tilde (~)
# - mehrere UVs lassen sich mit + verbinden
# - wenn sie mit * verbunden sind, wird zusätzlich der Interaktionseffekt berechnet
# That's it. 

# AV ~ UV1 + UV2
# Sprich: die AV wird vorhergesagt durch UV1 und UV2 (separat)

# AV ~ UV1 * UV2
# Sprich: die AV wird vorhergesagt durch UV1, UV2 und den Interaktionseffekt aus UV1/UV2

# Übertragen auf unseren Anwendungsfall:
# einschaetzung ~ text * geschlecht
aj_lm <-
  lm(einschaetzung ~ text*geschlecht,
     data = aj_manova_daten)

# Das war es schon. 
# Wir können uns die (M)AnOVa-Informationen anzeigen lassen.
# In den Zusatzinformationen des lm (gaaanz unten) steckt außerdem das (korrigierte) R².
anova_stats(aj_lm)
summary(aj_lm)




# In einige Fällen eines Experimentaldesigns kann auch eine Regression gefordert sein.
# Es gibt verschiedene Arten der Regression (zB logistische), wir beschränken uns hier auf die lineare.
# (Infos zu anderen Formen gibt es zB hier: https://bookdown.org/ndphillips/YaRrr/regression-on-non-normal-data-with-glm.html)
# Dafür nutzen wir noch einmal die Worlds-of-Journalism-Befragungsdaten aus tidycomm.
# Wir versuchen, Autonomie (AV) durch Erfahrung (UV1) und prekäre Arbeitsverhältnisse (UV2) vorhersagen.
woj_befragung <- 
  WoJ %>% 
  mutate(arbeitet_prekaer = if_else(temp_contract == 'Temporary', 1, 0, missing = 0))


# Wie bereits bei der MAnOVa definieren wir wieder ein "lineares Modell" (lm).
# Wie bereits bei der MAnOVa brauchen wir dafür die Formel-Beschreibung.
woj_lm <- lm(autonomy_selection ~ work_experience + arbeitet_prekaer,
             data = woj_befragung)
summary(woj_lm)


# Hier erhalten wir bereits zahlreiche zentrale Indikatoren:
# - "(Intercept)" ist das, was SPSS die "(Konstante)" nennt
# - je Variable erhalten wir das B ("Estimate"), den Standardfehler (Std. Error), T und die Signifikanz
# - außerdem erhalten wir Teststatistiken weiter unten (F, df und p)
# - das (korrigierte) R² finden wir ebenfalls fast ganz unten


# Es fehlen uns 
#  (1) ein Vergleich der beiden separat eingeschlossenen UVs
#  (2) Beta-Werte.


# Zunächst zum (1) Vergleich der verschiedenen Modelle.
# Das passiert hier nicht automatisch. Stattdessen müssen wir die verschiedenen Modelle selbst berechnen.
woj_lm_erfahrung <- lm(autonomy_selection ~ work_experience,
                       data = woj_befragung)
summary(woj_lm_erfahrung)

woj_lm_erfahrung_vertrag <- lm(autonomy_selection ~ work_experience + arbeitet_prekaer,
                               data = woj_befragung)
summary(woj_lm_erfahrung_vertrag)


# Anschließend können wir sie mithilfe einer AnOVa gegenüberstellen.
# (Hintergrundinfo: SPSS macht das genauso, bloß klammheimlich ...)
anova(woj_lm_erfahrung, woj_lm_erfahrung_vertrag)


# Was lernen wir?
# Model 2 ist nicht signifikant (p = .337) besser, wir bleiben also beim ersten Modell.


# Und was ist mit den (2) Beta-Werten?
# Die sind in vielen anderen Disziplinen für die Regression unüblich.
# Deshalb brauchen wir weiteres Paket, das schlicht "lm.beta" heißt (--> installieren).
library(lm.beta)
?lm.beta

# Das klingt einfach ...
lm.beta(woj_lm_erfahrung)




# Hier endet unser gemeinsamer R-Auftakt.
# Wer es bis hierhin geschafft hat, ist solide vorbereitet auf den weiteren Weg.
# Viel Erfolg :)