Kapitel 13 Statistische Tests

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.

13.1 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.

13.1.1 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:

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  qmsurvey$partnerinnen
#> W = 0.76247, p-value = 1.506e-11

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

13.1.2 Tests auf Varianzhomogenität

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

#> Levene's Test for Homogeneity of Variance (center = median)
#>       Df F value    Pr(>F)    
#> group  1  28.747 5.983e-07 ***
#>       93                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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.

Table 1: Levene's Test for Homogeneity of Variance (Brown-Forsythe Adaption)

Term df F p
ons 1 28.75 < 0.001
Residuals 93



13.2 \(\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:

Hattest du schon
einmal einen One
Night Stand 		Warst du schon
einmal in
psychotherapeutischer
Behandlung 		Total
Ja Nein
Ja 12
12
12.6 % 		26
26
27.4 % 		38
38
40 %
Nein 19
19
20 % 		38
38
40 % 		57
57
60 %
Total 31
31
32.6 % 		64
64
67.4 % 		95
95
100 %
χ2=0.000 · df=1 · φ=0.018 · p=1.000

observed values
expected values
% of total

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

#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  qmsurvey$ons and qmsurvey$behandlung
#> X-squared = 2.5447e-31, df = 1, p-value = 1

13.3 z- und t-Test

13.3.1 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:

Table 2: z-Test with alternative hypothesis: \(\mu_1 \neq\) 20

\(\mu_1\) alter SE z \(CI_{95\%}\) p Cohen's d Power
21.7 0.3 5.71 (13.22 - 30.18) < .001 0.57 1



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.

Table 3: One Sample t-test with alternative hypothesis: \(\mu_1 \neq\) 20

\(\mu_1\) alter df SE t \(CI_{95\%}\) p Cohen's d Power
21.7 101 0.38 4.49 (20.95 - 22.45) < .001 0.44 0.99



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.

13.3.2 Zwei Stichproben

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

#> 
#>  Welch Two Sample t-test
#> 
#> data:  partnerinnen by ons
#> t = 6.4053, df = 40.321, p-value = 1.231e-07
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#>   5.23016 10.05054
#> sample estimates:
#>   mean in group Ja mean in group Nein 
#>           9.605263           1.964912

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 nochmal auf, da gehen wir stärker auf Modelldefinition ein.

Alternativ haben wir da was in der toolbox:

Table 4: Welch Two Sample t-test with alternative hypothesis: \(\mu_1 \neq \mu_2\)

Diff \(\mu_1\) Ja \(\mu_2\) Nein t SE df \(CI_{95\%}\) p Cohen's d Power
7.64 9.61 1.96 6.41 1.19 40.32 (5.23 - 10.05) < .001 1.58 1



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)

13.4 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.

13.4.1 Wilcoxon & Mann-Whitney

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

#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  partnerinnen by ons
#> W = 2002, p-value = 2.043e-12
#> alternative hypothesis: true location shift is not equal to 0

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:

Table 5: Wilcoxon rank sum test with continuity correction with alternative hypothesis: \(M_1 \neq M_2\)

Difference \(M_1\) Ja \(M_2\) Nein W p
7 8 1 2002 < .001



13.4.2 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:

#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  partnerinnen by rauchen
#> Kruskal-Wallis chi-squared = 8.8281, df = 2, p-value = 0.01211

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

Table 6: Kruskal-Wallis Rank Sum Test

\(\chi^2\) df p
8.83 2 < .05



13.5 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.

13.5.1 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.

13.5.2 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

13.5.3 \(\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

13.5.4 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.

13.5.5 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 und hier ergänzend für diesen Fall.

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

13.5.6 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*)

13.6 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:

## 
##      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:

## 
##      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.

13.6.1 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.

13.6.2 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 Treatmentvarianzen direkt angegeben, ergo ist keine dedizierte Effektgrößenbestimmung notwendig.

  1. Die Namensgebung ist etwas irritierend, ja. Einfach merken: “homo” = “gleich”, “hetero” = “unterschiedlich”, und anscheinend für diesen Kontext: “Varianz” \(\approx\) “Skedastizität”↩