13-Statistische-Tests.Rmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---
# Statistische Tests

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

Hier kümmern wir uns um die meisten gängigen statistischen Tests aus QM2.  
Es sollte dazugesagt werden, dass wir Regression und ANOVA ausgeklammert haben um ihnen eigene Kapitel zu spendieren, der Kram ist nämlich einen Tacken komplexer als ein simpler t-Test.  
Als nächstes solltet ihr wissen, dass das hier nur ein *Auszug* aus gängigen statistischen Verfahren ist.  Was wir hier machen fällt unter *frequentistische Statistik*, und die coolen Kinder drüben in der *Bayesian*-Ecke finden das ziemlich langweilig.  
Und irgendwann, wenn die Autoren dieses Werks verstanden haben wie *Bayesian* Kram funktioniert, dann ergänzen wir das hier vermutlich auch noch.

## Tests auf Voraussetzungen

Die meisten der Tests hier setzen *Normalverteilung* und/oder *Varianzhomogenität* (aka *Homoskedastizität*) voraus, und auch wenn man diese Annahmen entweder mit *Robustheit* wegargumentieren oder im Falle der Normalverteilungsannahme mit ausreichend großem $N$ und dem zentralen Grenwzertsatz dezent ignorieren kann, gibt es natürlich auch Tests dafür.  
Weil eine mehr oder weniger objektive Testentscheidung immer dankbar ist.  

### Tests auf Normalvertielung

Normalverteilungstests testen eure empirische Verteilung auf Ähnlickeit zur Normalverteilung. Wer hätte das gedacht. Es gibt einen großen Haufen dieser Tests, und einige davon sind sogar richtig blöd.  
Dazu gehört zum Beispiel der *Kolmogorov-Smirnoff-Test*, der zwar richtig richtig kacke ist, aber da es der einzige Normalverteilungstest ist, den SPSS ohne manuellen Aufwand anbietet, ist er lächerlich weit verbreitet.  
Ja. Wenn ihr aus diesem Kapitel eins mitnehmt: Benutzt nicht den Kolmogorov-Smirnoff.

Wenn es um Teststärke geht hat der **Shapiro-Wilk**-Test die Nase vorn, der sollte also so der Standardtest sein, den ihr erstmal probiert. R hat den auch von Haus aus dabei:

```{r shapiro}
shapiro.test(qmsurvey$partnerinnen)
```

Solltet ihr mal andere Tests brauchen, hat das package `nortest` noch einige andere gängige:

```{r normtests, eval=FALSE}
library(nortest)

# Anderson-Darling (auch sehr gut)
ad.test(qmsurvey$partnerinnen)
# Shapiro-Francia (auch ganz gut)
sf.test(qmsurvey$partnerinnen)
# Lilliefors (okay)
lillie.test(qmsurvey$partnerinnen)
# Pearson's Chi^2-Anpassungstest (meh)
pearson.test(qmsurvey$partnerinnen)
# Kolmogorov-Smirnoff (einfach nicht benutzen)
ks.test(scale(qmsurvey$partnerinnen), y = pnorm)
```

### Tests auf Varianzhomogenität

Der gängigste Test auf *Varianzhomogenität* (oder auch *Homoskedastizität* bzw. analog *Heteroskedastizität* für *Varianzheterogenität*[^sked]) ist der **Levene Test**, den wir aus dem `car` package holen:

```{r levene_1, message=FALSE, warning=FALSE}
library(car)

leveneTest(partnerinnen ~ ons, data = qmsurvey)
```

Hier haben wir auf *Varianzheterogenität* getestet in den Gruppen, die durch `ons` (`"Hattest du schon einmal einen One-Night-Stand?"`) entstehen (`Ja` und `Nein`), mit dem abhängigen Merkmal `partnerinnen` (`"Wie viele Sexualpartner*innen hast du bisher gehabt?"`).  
Beachtet die Tilde (`~`) in der Funktion. Das ist eine von diesen *formulas*, die in R-Modelldefinitionen gerne auftauchen. Erstmal nicht weiter beachten, nur merken: Links steht die abhängige (intervallskalierte) Variable, und rechts steht die *Gruppierungsvariable*, in der Regel nominalskaliert.

```{r levene_2}
library(tadaatoolbox)

tadaa_levene(qmsurvey, partnerinnen ~ ons, print = "markdown")
```

## $\chi^2$

Zu einem sauberen $\chi^2$-Test gehört auch immer eine saubere Kontingenztabelle, und da greifen wir wie gehabt auf `sjPlot` zurück:

```{r chisq_xtab, message=FALSE, warning=FALSE}
library(sjPlot)

sjt.xtab(qmsurvey$ons, qmsurvey$behandlung,
         show.exp = TRUE, show.legend = TRUE, show.cell.prc = TRUE)
```

Den $\chi^2$-Test finden wir dann in *base R*, no packages required:

```{r chisq_base_r, message=FALSE, warning=FALSE}
chisq.test(qmsurvey$ons, qmsurvey$behandlung)
```

## z- und t-Test

### Eine Stichprobe

Einstichprobentests sind… sagen wir mal, _unüblich_.  
Solltet ihr doch mal in die Verlegenheit kommen einen **z-Test** mit bekannter Streuung durchführen zu müssen, so seid ihr vermutlich aus Versehen in der Epidemiologie gelandet.  
In R gibt's dafür jedenfalls von Haus aus keine Funktion (von der ich weiß), deswegen haben wir den Anwendungsfall auch in der `tadaatoolbox` abgedeckt:

```{r z-Test}
tadaa_one_sample(qmsurvey, x = alter, mu = 20, sigma = 3, print = "markdown")
```

Hier testen wir die Variable `alter` auf die Mitte $\mu = 20$ mit einer angenommenen Streuung von $\sigma = 3$.  
Ist das sinnvoll?  
Keine Ahnung, es ist ein Beispiel.  

Eine etwas gängigere Version ist ein Einstichproben-t-Test.  
Wieso auch nicht.  

```{r t-test_onesample}
tadaa_one_sample(qmsurvey, x = alter, mu = 20, print = "markdown")
```

Wir hätten den t-Test auch mit der R-eigenen `t.test`-Funktion machen können, die gucken wir uns als nächstes an.

### Zwei Stichproben

R bringt die t-Test-Funktion von Haus aus mit, die ist zwar nicht besonders schön, aber sie macht alles wichtige:

```{r t-test_twosample}
t.test(partnerinnen ~ ons, data = qmsurvey)
```

**Wichtige Argumente**:

- `var.equal`: Varianzen gleich? `TRUE`/`FALSE`, für Varianzhomogenität
    - Im Zweifelsfall lieber auf `FALSE` (default) lassen, sicherheitshalber.
- `paired`: Wenn `TRUE` wird ein *abhängiger* t-Test gemacht, bei `FALSE` (default) ein *unabhängiger*.

Beachtet bei der Funktion die Verwendung der Tilde: `~`, die zur *Modelldefinition* verwendet wird.  
In diesem Fall haben wir ein intervallskaliertes Merkmal (links von der Tilde), und das unabhängige nominalskalierte Merkmal (rechts von der Tilde). Die Tilde raucht in der [Regression] und [ANOVA](#ANOVA) nochmal auf, da gehen wir stärker auf Modelldefinition ein.

Alternativ haben wir da was in der toolbox:

```{r t-test_twosample_tadaa}
tadaa_t.test(qmsurvey, response = partnerinnen, group = ons, paired = FALSE, print = "markdown")
```

Die Funktion macht im Vergleich zu `t.test()` noch folgendes:  

- Automatische Varianzhomogenitätserkennung via eingebautem Levene-Test
- Automatische Effektgrößenberechnung (*Cohen's d*)
- Automatische Teststärkenberechnung via `pwr`-package (Spalte `Power`)

## Nichtparametrische Alternative

Sollte euch mal ein Skalenniveau verloren gehen oder eure Annahmen kaputtgegangen sein, könnt ihr immer noch auf *nichtparametrische Tests* umsteigen. Die spielen zwar in der Regel nicht bei den coolen Kindern mit, aber naja, sie sind da.

### Wilcoxon & Mann-Whitney

Die nichtparametrischen Alternativen zu abhängigen und unabhängigen t-Tests sind von Haus aus dabei mit der Funktion `wilcox.test`:

```{r wilcox_test}
wilcox.test(partnerinnen ~ ons, data = qmsurvey, paired = FALSE)
```

Und schon haben wir einen Rangsummentest gemacht.  
Wow.  
Wollt ihr einen Wilcoxon-Vorzeichentest? Dann tauscht `paired = FALSE` durch `paired = TRUE` aus.  
Kabläm.  
Die restlichen Argumente sind mehr oder weniger analog zu `t.test` und können der Hilfe unter `?wilcox.test` entnommen werden.  

Wenn ihr auch hier die *Tadaa!*-Variante wollt:

```{r wilcox_tadaa}
library(tadaatoolbox)

tadaa_wilcoxon(qmsurvey, response = partnerinnen, group = ons, paired = FALSE, print = "markdown")
```

### Kruskal-Wallis

Falls euch die ANOVA zu parametrisch ist, ist der Kruskall-Wallis praktisch wie für euch gemacht.  
Ich meine, es ist ja nicht so als ob ein sauberes Regressionsmodell nicht auch meistens ausreichen würde, aber naja, hier so:

```{r kruskal}
kruskal.test(partnerinnen ~ rauchen, data = qmsurvey)
```

…und ja, natürlich gibt es eine Alternative in der toolbox:

```{r kruskal_tadaa}
tadaa_kruskal(partnerinnen ~ rauchen, data = qmsurvey, print = "markdown")
```


## Kritische Werte und Theoretisch Verteilungen

Aus der Vorlesung und den Tutorien seid ihr es gewohnt die kritischen Werte für eure Hypothesentests aus irgendwelchen aus Büchern eingescannten Tabellen zu ziehen.  
Da die 1980er Jahre mittlerweile vorbei sind, scheint es mir angemessen diesen Zustand zu ändern, namentlich durch die Berechnung beliebiger kritischer Werte über die in `R` mitgelieferten Verteilungsfunktionen. 

Zentral zum Verständnis ist hierbei, dass es sich bei den kritischen Werten letztlich nur um Quantile handelt. Der (linksseitige) kritische Wert der Normalverteilung für ein $\alpha = 0.05$ entspricht demnach dem 5%-Quantil. In `R` gibt es standardmäßig Funktionen für die Quantile vieler gängiger theoretischer Verteilungen, von der Normalverteilung über die t-Verteilung zur F-Verteilung und sogar für die diskreten Verteilungen der nichtparametrischen Tests.

In diesem Dokument findet ihr kurze Hilfestellung zu den jeweiligen Testverteilungen. Eigenrecherche von Vorteil, ich will euch auch nicht alles vorkauen.

### Verteilungsfunktionen in R

R liefert euch diverse theoretische Verteilungen mit ihren Dichtefunktionen, kumulativen Wahrscheinlichkeitsverteilungen, Quantilsfunktionen und sogar Zufallsgeneratoren auf dem Silbertablett. Eine Übersicht könnte ihr euch hier verschaffen: `?Distributions`

Viel Spaß damit.

### z-Test (Normalverteilung)

Die wohl wichtigste theoretische Vertielung der Statistik. Siehe auch `?Normal`.  

* `dnorm` gibt die Dichtefunktion (relevant um die Funktion zu plotten)
    * `dnorm(0)` = Funktionswert der Normalverteilung bei x = 0 (0.3989423)
* `pnorm` die Wahrscheinlichkeitsfunktion (relevant für die Berechnung von Signifikanzniveaus)
    * **Kumulative Wahrscheinlichkeit**, im Grunde das Integral der Dichtefunktion
    * `pnorm(1)` = Wahrscheinlichkeit von $x = -\infty$ bis $x = 1$
    * `pnorm(1, lower.tail = FALSE)` = Wahrscheinlichkeit von $x = \infty$ bis $x = 1$
* `qnorm` gibt Quantile (relevant für kritische Werte)
    * `qnorm(0.05)` = 5%-Quantil = x-Wert, der in 5% der Fälle erreicht wird (-1.6448536)
    * `qnorm(0.05, mean = 90, sd = 2)` = Kritischer x-Wert für 5% bei einer Normalverteilung mit $\mu = 90$ und $\sigma = 2$
* `rnorm` gibt normalverteilte Zufallszahlen aus

### $\chi^2$-Test ($\chi^2$-Verteilung)

Die $\chi^2$-Verteilung. Siehe auch `?Chisquare`.  
Auch hier gilt:

* `dchisq` gibt die Dichtefunktion (relevant um die Funktion zu plotten)
    * `dchisq(5, df = 1)` = Funktionswert der $\chi^2$-Verteilung bei x = 5 und $df = 1$ (`0.01464498`)
* `pchisq` die Wahrscheinlichkeitsfunktion (relevant für die Berechnung von Signifikanzniveaus)
    * **Kumulative Wahrscheinlichkeit**, im Grunde das Integral der Dichtefunktion
    * `pchisq(1, df = 1)` = Wahrscheinlichkeit von $x = 0$ bis $x = 1$ bei $df = 1$
* `qchisq` gibt Quantile (relevant für kritische Werte)
    - Beispiel für kritischen Wert bei 95% Sicherheit bzw. $\alpha = 0.05$:
    * `qchisq(0.95, df = 1)` = 95%-Quantil = x-Wert, der in 95% der Fälle erreicht wird (`3.841459`)
* `rchisq` gibt $\chi^2$-Verteilte Zufallszahlen aus

### t-Testfamilie (t-Verteilung)

Die t-Verteilung. Siehe auch `?TDist`.  
Auch hier gilt:

* `dt` gibt die Dichtefunktion (relevant um die Funktion zu plotten)
    * `dt(5, df = 10)` = Funktionswert der t-Verteilung bei x = 5 und `df = 10` $(3.9600106\times 10^{-4})$
* `pt` die Wahrscheinlichkeitsfunktion (relevant für die Berechnung von Signifikanzniveaus)
    * **Kumulative Wahrscheinlichkeit**, im Grunde das Integral der Dichtefunktion
    * `pt(1, df = 10)` = Wahrscheinlichkeit von $x = -\infty$ bis $x = 1$ bei `df = 1`
* `qt` gibt Quantile (relevant für kritische Werte)
    * `qt(0.05, df = 10)` = 5%-Quantil = x-Wert, der in 5% der Fälle erreicht wird
* `rt` gibt t-Verteilte Zufallszahlen aus

Das Prinzip ist das gleiche für alle Verteilungen, die Präfixe `d, p, q, r` haben jeweils die gleiche Bedeutung.

### Wilcoxon-Test (Vorzeichentest, Signrank Test)

Die theoretische Verteilung für den Wilcoxon-Test findet ihr mittels `?SignRank`.  
Auch hier gilt:

* `dsignrank` gibt die Dichtefunktion (relevant um die Funktion zu plotten)
* `psignrank` die Wahrscheinlichkeitsfunktion (relevant für die Berechnung von Signifikanzniveaus)
* `qsignrank` gibt Quantile (relevant für kritische Werte)
* `rsignrank` gibt Zufallszahlen basierend auf der Verteilung

Die praktische Berechnung der ktitischen Werte gestaltet sich bei den nichtparametrischen Tests etwas weniger intuitiv als bei den parametrischen Tests, da es sich um diskrete Verteilungen handelt.  
Eine Erklärung dazu findet ihr [hier auf Stackoverflow](http://stats.stackexchange.com/questions/32445/wilcoxon-mann-whitney-critical-values-in-r) und [hier ergänzend für diesen Fall](http://stats.stackexchange.com/questions/139317/critical-value-for-wilcoxon-one-sample-signed-rank-test-in-r). 

Die kurze Version:

* Kritischer Wert = `qsignrank(1 - alpha, n) - 1`
* Kritischer Wert für einen einseitigen Test, $\alpha = 0.05$ und $n = 8$:
    * `qsignrank(0.95, 8) - 1` = 29

### U-Test / Mann-Whitney-(Wilcoxon)-U-Test / Rangsummentest

Analog zum Signed Rank Test, nur mit anderer Verteilung. Grundlage ist hier `?Wilcoxon`.

* `dwilcox` gibt die Dichtefunktion (relevant um die Funktion zu plotten)
* `pwilcox` die Wahrscheinlichkeitsfunktion (relevant für die Berechnung von Signifikanzniveaus)
* `qwilcox` gibt Quantile (relevant für kritische Werte)
* `rwilcox` gibt Zufallszahlen basierend auf der Verteilung

* Kritischer Wert = `qwilcox(1 - alpha, m, n) - 1`
* `m` und `n` sind dabei der Größen der beiden Stichproben
* Kritischer Wert für einen einseitigen Test, $\alpha = 0.05$ und $n_1 = 8$, $n_2 = 7$:
    * `qwilcox(0.95, 8, 8) - 1` = 47
    
Bei den Signifikanzniveaus der nichtparametrischen Tests seid ihr natürlich ganz vorsichtig, da ihr natürlich wisst, dass ihr hier nicht ohne Weiteres gültige Signifikanzniveaus angeben könnt. (\*hust\*)

## Poweranalyse

Für die Power-Analyse gibt es Möglichkeiten das ganze mit R statt G-Power zu machen, was natürlich unter Umständen etwas schwieriger nachzuvollziehen ist.  
Allgemeine Grundlage ist das package `pwr`.  
Für den z-Test sieht das wie folgt aus:

```r
library(pwr)
pwr.norm.test(d = 0.2, sig.level = 0.05, power = 0.8, alternative = "greater")
```

```
## 
##      Mean power calculation for normal distribution with known variance 
## 
##               d = 0.2
##               n = 154.5639
##       sig.level = 0.05
##           power = 0.8
##     alternative = greater
```

wobei gilt: 

* `d` = Effektgröße $\Delta$
* `sig.level` = Testniveau $\alpha$
* `power` = Teststärke $1-\beta$
* `alternative` = Testrichtung ("`greater`", "`less`" oder "`two.sided`")

Erklärungen in der R-Hilfe unter `?pwr.norm.test`

Das Ergebnis stimmt mit meiner semi-händischen Berechnung überein.
Die Funktion nimmt neben den genannten Argumenten auch den Stichprobenumfang n. Setzt man n explizit fest, beispielsweise auf `n = 60`, und lässt dafür die Teststärke weg, erhält man folgendes:


```r
pwr.norm.test(d = 0.2, sig.level = 0.05, n = 60, alternative = "greater")
```


```
## 
##      Mean power calculation for normal distribution with known variance 
## 
##               d = 0.2
##               n = 60
##       sig.level = 0.05
##           power = 0.4618952
##     alternative = greater
```

Ergo kann man diese Funktion für den vierseitigen Zusammenhang zwischen $\alpha$, $\beta$, `n` und Effektgröße $\Delta$ nutzen:  
Das Argument, das fehlt, wird auf Basis der anderen Argumente berechnet.

### t-Test

Bei den t-Tests ist das ganze etwas komplexer, allein wegen der Vielfalt an verschiedenen t-Tests und deren Umgang mit der Gesamtvarianz.  
Anbei die verschiedenen Funktionen:

* `pwr.t.test` für
    * Einstichprobentest: `pwr.t.test(…, type = "one.sample")`
    * 2 unabhängige Stichproben mit gleichem $n$: `pwr.t.test(…, type = "two.sample")`
    * Abhängige Stichproeben: `pwr.t.test(…, type = "paired")`
* `pwr.t2n.test` für
    * Unabhängige Stichprobene mit ungleichem $n$
    * Beispiel: `pwr.t2n.test(n1 = 12, n2 = 15, …)`
    
Und natürlich ist in allen Fällen die Testrichtung zu beachten via `alternative = <Richtung>`, wobei `<Richtung>` eben entweder `"less"`, `"greater"` oder `"two.sided"` ist.

### ANOVA

Auch für die ANOVA ist die Poweranalyse in R möglich, dafür wird die Funktion `pwr.anova.test` benutzt.  
Siehe auch: `?pwr.anova.test`  

Wichtig ist jedoch: Die Teststärkeberechnung ist im Zweifelsfall nur für *balancierte Designs* korrekt, das heißt wenn ihr ein Versuchsdesign mit unterschiedlich großen Gruppengrößen habt oder die Vorraussetzungen der ANOVA stark verletzt sind, sind die Ergebnisse der Teststärkenberechnung ggf. schlicht falsch. Leider gibt es da kein einfaches Allheilmittel soweit ich weiß, und euch bleibt nur die Option einfach von vornherein so sauber wie möglich zu arbeiten.

Eine Sache zur **Effektgröße** bei der ANOVA:  
Ihr kennt vermutlich $\eta^2$, aber die Funktionen zur Effektgrößenberechnung setzen **Cohen's f** voraus. Glücklicherweise gibt euch unsere ANOVA-Funktion aus der *tadaatoolbox* `tadaa_aov` beide Effektgrößen aus, aber wenn ihr mal $f$ selber berechnen müsst:

$$\text{Cohen's f} = \sqrt{\frac{\eta^2}{1- \eta^2}}$$

Für ein $\eta^2 = 0.2$ gibt das also ein $f = \sqrt{\frac{0.2}{1 - 0.2}} = 0.5$

Aber zurück zu `pwr.anova.test`:

* Beispielaufruf: `pwr.anova.test(k = 3, n = 3, f = 2, sig.level = 0.05)`
* `f` ist hierbei allerdings nicht der F-Wert aus dem F-Test, sondern die *Effektgröße* **Cohen's f**
* `n` ist die Größe jeder einzelnen Gruppe, nicht das Gesamt-N
    - Merke: Die Funktion nimmt **gleichgroße Gruppen** an!

Eine weitere Möglichkeit zur Berechnugn der Teststärke ist die Funktion `power.anova.test` aus den Standardpaketen (sprich nicht aus dem `pwr`-package):

* `power.anova.test(groups = 3, n = 3, between.var = 34.1, within.var = 2.1, sig.level = 0.05, power = NULL)`
* `n` ist die Größe jeder einzelnen Gruppe, nicht das Gesamt-N

Hier werden Binnen- und Treatment*varianzen* direkt angegeben, ergo ist keine dedizierte Effektgrößenbestimmung notwendig.

<!-- Footnotes -->
[^sked]: Die Namensgebung ist etwas irritierend, ja. Einfach merken: "homo" = "gleich", "hetero" = "unterschiedlich", und anscheinend für diesen Kontext: "Varianz" $\approx$ "Skedastizität"