12-Deskriptive-Statistik.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---
# (PART) Statistik {-}

# Deskriptive Statistik

```{r include=FALSE}
knitr::opts_knit$set(unnamed.chunk.label = "descriptives_")
knitr::opts_chunk$set(error = FALSE, warning = FALSE, message = FALSE)
library(tadaatoolbox)
library(irr)
library(sjstats)
```

In diesem Abschnitt widmen wir uns univariaten und bivariaten *deskriptiven Statistiken*, vom Mittelwert und Standardabweichung über $\chi^2$, *Somers' D* und Korrelationskoeffizienten.  
Wir halten uns (so grob) an die Struktur aus QM1 und QM2, das heißt: Wir kümmern uns beispielsweise zuerst um $\chi^2$ als deskriptive Statistik um dann auf *Cramer's V* zu sprechen zu kommen, und widmen uns erst im folgenden Kapitel dem $\chi^2$-Test als inferenzstatistisches Werkzeug.

## Univariate Statistiken

*Univariat* ist *fancy sprech* für "nur eine Variable", das heißt es geht erstmal noch nicht um Zusammenhänge irgendeiner Art, sondern primär darum wie ein Merkmal **verteilt** ist.  
Die dazugehörigen Kennwerte und Formeln sind sozusagen das kleine 1x1 der Statistik, dementsprechend sind sie auch in R schnell und einfach erreichbar.  
Als Beispieldatensatz nehmen wir wieder den *Game of Thrones Deaths*-Datensatz und den `qmsurvey`-Datensatz, beide aus dem Kapitel [Datenimport]:

```{r ch12_Daten_einlesen, message=FALSE, warning=FALSE}
library(readr)
qmsurvey  <- readRDS("data/qm_survey_ss2017.rds")
gotdeaths <- read_csv("data/got_deaths.csv", col_types = cols())
```

### Mittelwert, Median, Modus

Mittelwert und Median sind einfach: 

```{r ch_11_mean_median}
mean(qmsurvey$zufrieden)
median(qmsurvey$zufrieden)
```

Die Funktionen `mean` und `median` sind da wenig überraschend.  
Aber was machen wir mit dem Modus? In R an sich gibt es keine `modus`-Funktion, deswegen haben wir an dieser Stelle zwei Optionen:

1. Häufigkeitstabelle erstellen und ablesen.
2. Ein package benutzen, dass eine Funktion für den Modus hat.

Eine Häufigkeitstabelle bekommen wir sehr einfach mit `table`:

```{r freq_table}
table(qmsurvey$zufrieden)
```

Das ist zwar optisch wenig beeindruckend, aber wir können immerhin einfach erkennen, dass der Wert `4` die größte absolute Höufigkeit hat. Passt.

Alternativ können wir die `tadaatoolbox` benutzen:

```{r modus_toolbox}
library(tadaatoolbox)

modus(qmsurvey$zufrieden)
```

Vollkommen unerwartet kommen wir zum gleichen Ergebnis.  
Who'da thunk.  
Für die EnthusiastInnen:  
Die `modus`-Funktion macht übrigens auch nichts anderes als eine Häufigkeitstabelle via `table`, aus der dann das Maximum extrahiert wird.  
Praktisch nur `names(table(x)[table(x) == max(table(x))])` mit ein bisschen Glitzer und so.

Wenn wir das jetzt noch in schön und übersichtlich haben wollen, können wir auf `sjmisc` zurückgreifen:

```{r freq_sjt, warning=FALSE, message=FALSE}
library(sjmisc)

frq(qmsurvey$zufrieden, out = "viewer")
```

Da fällt eine Tabelle raus, die uns Median, die Häufigkeiten, fehlende Werte und in der Fußzeile gleich auch gesamt-N, Mittelwert und Standardabweichung anzeigt.

### Variabilität

Hier ist die Situation ähnlich wie im vorigen Abschnitt: Varianz und Standardabweichung sind geschenkt, Variation ist seltsam.

```{r var_sd}
var(qmsurvey$zufrieden)
sd(qmsurvey$zufrieden)
```

Wenig überraschend. Solange ihr euch merken könnt, dass *Standardabweichung* auf englisch ***s**tandard **d**eviation* heißt, sollte das kein Problem sein.  
Die Variation hingegen wird so selten als Maß benutzt, dass es dafür keine dedizierte Funktion gibt – analog der Modus-Sache von weiter oben. Da das allerdings _wirklich_ kein Mensch braucht, bleibt euch nichts anderes übrig als die Varianz zu berechnen und mit $n-1$ zu multiplizieren.  

¯\\\_(ツ)_/¯

Okay, ihr könnt natürlich auch *ganz cool* sein und euch eine **eigene Funktion** schreiben! Wie die großen!

So etwa:

```{r variation}
variation <- function(x) {
  # NA rauskicken
  x <- x[!is.na(x)]
  # Variation berechnen
  variation <- sum((x - mean(x))^2)
  # Und zurückgeben
  return(variation)
}
```

Das war... Definitiv etwas, das wir gerade getan haben!  

```{r variation_2}
variation(qmsurvey$zufrieden)
```

Funktioniert.  
Braucht nur wie gesagt kein Mensch, aber jetzt habt ihr immerhin mal so grob gesehen, wie Funktionen schreiben in R funktioniert.  
'Tis not witchcraft.  
Aber okay, ich schweife ab.

### Variabilität auf Ordinalskala

Auf Ordinalskala haben wir den Mittelwert nicht zur Verfügung. Das macht uns alle sehr traurig.  
Aus diesem Grund müssen wir auf den Median ausweichen, den wir dann noch dezent mit Quantilen garnieren können.  
Wenn wir nochmal die Variable `Zufriedenheit` aus `qmsurvey` nehmen:

```{r median}
median(qmsurvey$zufrieden)
```

Wenig beeindruckend.  
Für Quantile:

```{r quantiles}
quantile(qmsurvey$zufrieden)
```

Wir können auch beliebige andere Quantile berechnen, die `quantile`-Funktion gibt uns nur per default Quartile, weil die am gängigsten sind:

```{r quantiles_2}
quantile(qmsurvey$alter, probs = c(.333, .666))
```

Über `probs` geben wir der Funktion einen *Vektor aus Wahrscheinlichkeiten* für die entsprechenden Quantile. Hier haben wir also das 33%- und das 66%-Quantil berechnet.  
Um beispielsweise Dezile zu erhalten brauchen wir also einen Vektor der Wahrscheinlichkeiten von 0 bis 10 in 10% oder 0.1er Schritten, den wir einfach via `seq` erstellen können:

```{r quantiles_3}
quantile(qmsurvey$alter, probs = seq(from = 0, to = 1, by = 0.1))
```

Ein weiteres gängiges Maß für die Variabilität ist die **Spannweite**, die auf Ordinalskala und aufwärts sinnvoll sein kann.  
Gemeint ist damit nichts anderes als der Abstand vom kleinsten zum größten Wert, also nichts besonders komplexes. In R können wir da die `min` und `max`-Funktionen für die beiden Werte benutzen, oder `range` für beide zusammen:

```{r range}
min(qmsurvey$alter)
max(qmsurvey$alter)
range(qmsurvey$alter)

max(qmsurvey$alter) - min(qmsurvey$alter)
```

Die Altersspanne in QM (und damit so grob im Psychologiestudiengang) ist also 19 Jahre.  
Bahnbrechende Erkenntnisse die wir hier einfach so raushauen.

## Bivariate Statistiken

### Nominalskala

Auf Nominalskala ist $\chi^2$ unser bester Freund. Leider ist es auch ein etwas übelriechender, buckliger Freund, den wir nur selten zu uns nach Hause einladen – zumindest in Deskriptivstatistikland.

Als Beispiel nehemn wir mal *Geschlecht* und *Cannabiskonsum*:

```{r nomstats_1}
library(sjPlot)

sjt.xtab(qmsurvey$cannabis, qmsurvey$behandlung,
         show.exp = TRUE, show.legend = TRUE)
```

Die Berechnung in R ist einfach:

```{r nomstats_2}
chisq.test(qmsurvey$cannabis, qmsurvey$behandlung)
# Reihenfolge egal
chisq.test(qmsurvey$behandlung, qmsurvey$cannabis)
```

Das Output ist nicht so hübsch und etwas clunky und beinhaltet auch p-Wert und Freiheitsgrade, die euch eigentlich nur in einem inferenzstatistischen Kontext interessieren.  
Darüber hinaus kann es sein, dass der berechnete Wert nicht mit händischer Berechnung übereinstimmt.  
Der Grund dafür ist [Yates' Kontinuitätskorrektor](https://en.wikipedia.org/wiki/Chi-squared_test#Yates.27s_correction_for_continuity), die nur bei 2x2-Tabellen greift. Solltet ihr die Korrektur explizit deaktivieren wollen, könnt ihr das Argument `correct = FALSE` setzen.

Für kompakteres Output haben wir eine entsprechende Funktion in der *tadaatoolbox*:

```{r nomstats_3, message=FALSE}
library(tadaatoolbox)

nom_chisqu(qmsurvey$behandlung, qmsurvey$cannabis)
```

Als nächstes haben wir die auf $\chi^2$-basierenden Statistiken **Phi** ($\phi$) und **Cramer's V**, die entsprechenden Funktionen kommen auch aus der *tadaatoolbox*:

```{r nomstats_4}
nom_phi(qmsurvey$behandlung, qmsurvey$cannabis)
nom_v(qmsurvey$behandlung, qmsurvey$cannabis)
```

Oh nein, Phi funktioniert nicht. Wie unerwartet.  
Ihr mögt euch erinnern, dass Phi nur für 2x2-Tabellen definiert ist, wobei Cramer's V bei beliebigen Tabellendimensionen funktioniert – darüberhinaus sind die beiden Statistiken im Falle von 2x2-Tabellen auch noch identisch, weshalb sich die Frage aufwirft, wieso überhaupt jemand Phi benutzen sollte – aber nun gut, die Funktionen sind da.

Als letztes gibt es noch den **Kontingenzkoeffizient C** – keine Ahnung was der für eine Daseinsberechtigung hat, aber wenn ihr ihn mal braucht:

```{r nomstats_5}
nom_c(qmsurvey$behandlung, qmsurvey$cannabis)
```

Sagte ich *als letztes*? Das war gelogen. Es gibt noch **Lambda** ($\lambda$) in drei Geschmacksrichtungen:

1. "Sorum" (Variablen wie definiert)
2. "Andersrum" (x- und y-Variable vertauschen)
3. Symmetrisch

```{r nomstats_6}
nom_lambda(qmsurvey$behandlung, qmsurvey$cannabis)
nom_lambda(qmsurvey$behandlung, qmsurvey$cannabis, reverse = TRUE)
nom_lambda(qmsurvey$behandlung, qmsurvey$cannabis, symmetric = TRUE)
```

Es gibt zusätzlich die Funktion `assocstats` aus dem `vcd`-Package. Da gebt ihr allerdings nicht beide Variablen einzeln an, sondern einen `table` der beiden Variablen:

```{r nomstats_7}
library(vcd)

freq_table <- table(qmsurvey$behandlung, qmsurvey$cannabis)
freq_table

assocstats(freq_table)
```
So könnt ihr die gängigen Statistiken auf einen Blick sehen.  
Alternativ haben wir dafür natürlich auch was in der *tadaatoolbox*:

```{r nomstats_8}
tadaa_nom(qmsurvey$behandlung, qmsurvey$cannabis, print = "markdown")
```

Das Argument `print = "markdown"` ist nur für das Output in diesem Buch bzw. in RMarkdown-Dokumenten relevant, für die Nutzung in der Konsole könnt ihr das weglassen.  
Aber ihr benutzt ja jetzt eh alle brav RMarkdown, seit ihr das entsprechende Kapitel dazu hier überflogen habt, ja?  
Ja.  
Brav.

### Ordinalskala

Auf Ordinalskala haben wir primär **Gamma**, **Somers' D** und eine Handvoll **Tau** abzufrühstücken.  
Um's einach zu halten bedienen wir uns auch hier wieder der *tadaatoolbox*:

```{r ordstats_1}
ord_gamma(qmsurvey$lernen, qmsurvey$zufrieden)

ord_somers_d(qmsurvey$lernen, qmsurvey$zufrieden)
ord_somers_d(qmsurvey$lernen, qmsurvey$zufrieden, reverse = TRUE)
ord_somers_d(qmsurvey$lernen, qmsurvey$zufrieden, symmetric = TRUE)
```

Wobei **Somers' D** analog Lambda drei Geschmacksrichtungen hat.  

Oder auch hier, alle zusammen:

```{r ordstats_2}
tadaa_ord(qmsurvey$lernen, qmsurvey$zufrieden, print = "markdown")
```

Da sind dann auch gleich die Taus mit drin. Die könnt ihr auch einzeln bekommen:

```{r ordstats_3}
ord_tau(qmsurvey$lernen, qmsurvey$zufrieden, tau = "a")
ord_tau(qmsurvey$lernen, qmsurvey$zufrieden, tau = "b")
ord_tau(qmsurvey$lernen, qmsurvey$zufrieden, tau = "c")
```

Da $\tau_b$ auch als Korrelationskoeffizient durchgeht, gibts's den auch via `cor`:

```{r ordstats_4}
cor(qmsurvey$lernen, qmsurvey$zufrieden, method = "kendall")
```

### Intervallskala

Hier fangen die spaßigen Dinge an. Korrelationen und so.  
Und was ist Korrelation?  
Die standardisierte Kovarianz, nämlich.  
Und was ist Kovarianz?

```{r intstats}
cov(qmsurvey$alter, qmsurvey$beziehungen, use = "complete.obs")
```

Nämlich.  
`cov` gibt uns die **Kovarianz**, wenig überraschend, und über `use = "complete.obs"` spezifizieren wir, dass Beobachtungspaare mit fehlenden Werten (`NA`) ausgelassen werden sollen. Eine Alternative wäre `"pairwise.complete.obs"`, was in diesem Fall keinen Unterschied macht, aber prinzipiell wollt ihr eins von beidem benutzen um zu verhindern, dass ihr als Ergebnis nur ein trauriges `NA` bekommt.

Als nächstes also die Korrelationskoeffizienten:

- **Pearson's** $r$
- **Spearman's** $\rho$

Beide bekommen wir mit `cor`, abhängig vom `method`-Argument. Der Default ist `"pearson"`

```{r intstats_3}
cor(qmsurvey$alter, qmsurvey$beziehungen, use = "complete.obs")
cor(qmsurvey$alter, qmsurvey$beziehungen, use = "complete.obs", method = "spearman")
```

Über `method = "kendall"` bekommen wir auch Kendall's $\tau_b$. Just in case.

### Weitere Statistiken

Es gibt noch viele andere Statistiken, die für euch mehr oder weniger spannend sind. Viele andere sind für euch sogar vollkommen irrelevant.  
Point being: Das hier 

#### $\eta^2$

Ihr lernt $\eta^2$ in QM1 als deskriptive Statistik kennen, wo es als nicht-lineare Alternative zu $r$ oder $\rho$ verwendet wird. Das ist auch okay so, aber generell findet $\eta^2$ eher als Effektgröße bei der [ANOVA](#ANOVA) Anwendung, weshalb auch R eher diesen Kontext erwartet.  

Wenn ihr $\eta^2$ aber "einfach so" haben wollt, braucht ihr etwas Umweg:

```{r eta_sjstats}
library(sjstats)

anova <- aov(alter ~ alkohol, data = qmsurvey)
eta_sq(anova)
```

Wir haben hier zuerst ein *ANOVA*-Modell mittels `aov` erstellt, wo wir die Variable `alter` aufgeteilt haben in Gruppen abhängig von der Variable `alkohol`. Dann haben wir `eta_sq` aus `sjstats` für die Berechnung von $\eta^2$ benutzt, und... sind damit fertig.

Als Bonus-Feature haben wir da seit der Version `0.16.0` auch was passendes in der `tadaatoolbox`:

```{r eta_toolbox}
library(tadaatoolbox)

etasq(alter ~ alkohol, data = qmsurvey)
```

Unter der Haube benutzen wir da übrigens `DescTools::EtaSq`. Ihr seht, es gibt da etliche Möglichkeiten, wie so oft in R.

#### Cohen's Kappa ($\kappa$)

```{r}
library(irr)
data(anxiety) # Beispieldatensatz aus irr package laden
head(anxiety)

# Nur zwei Spalten / Rater für Kappa
anxiety_2 <- anxiety[c(1, 2)]

kappa2(anxiety_2)
```

#### Kendall's $W$

```{r}
library(irr)

kendall(anxiety)
```