Skip to Tutorial Content

Hinweise

Eine R-Nutzerin, die bereits vorher einen R-Kurs belegt hat, bewertete dieses Tutorial insgesamt mit einer Schwierigkeit von 6 (0=sehr leicht, 10=sehr schwer). Sie brauchte für dieses Tutorial ungefähr 1h10min. Klicke auf “Nächstes Kapitel” und es geht los.

Einführung

Leider musst Du Dich heute noch weiter mit statistischen Tests beschäftigen, genauer gesagt nonparametrischen Tests. Diese verwenden wir, wenn die untersuchten Variablen kein Intervallskalenniveau haben oder die Voraussetzungen für einen statistischen Test nicht erfüllt sind (z. B. keine Normalverteilung in der Population). Dazu gibt es zwei Dinge zu sagen: (1) in der Psychologie wird das Intervallskalenniveau so gut wie nie erreicht – die Psychologen tun aber gerne trotzdem so als ob, sie “spielen” Intervallskalenniveau. (2) Bootstrapping und Permutationstests wären eigentlich die bessere Alternative um das Problem zu lösen, werden aber erst im Masterstudium ausführlich behandelt.

Man könnte vermuten, dass es es keinen guten Grund gibt statistische Tests so zu lernen wie wir das gerade tun. Stattdessen könnte man direkt mit Bootstrapping anfangen und hätte die komplette Inferenzstatistik in 2-3 Vorlesungen abgehandelt und alles unter dem Hut eines einzigen Konzepts, statt einem Dutzend verschiedener Tests. Warum tun wir das (noch) nicht? Weil in der Literatur immer noch die klassischen Tests verwendet werden und Du für Dein empirisch-experimentelles Praktikum und Deine Abschlussarbeit die Literatur lesen und verstehen musst. Die Trägheit des Systems ist also Schuld daran, dass Du Deine kostbare Lebenszeit verschwenden musst. Nun ja, ganz verschwendet ist sie vielleicht nicht, denn zumindest lernst Du hier die Universalsprache R.

Da es also wieder etwas trocken wird, hab ich als kleinen Höhepunkt die Poweranalyse in dieser Sitzung integriert. Außerdem schauen wir uns auch kurz die Lowess-Kurve und Partialkorrelationen an. Diese Themen sind etwas willkürlich gewählt, aber für eine Erfrischung gut geeignet.

Wir verwenden übrigens wieder die Vorlesungsdaten:

data <- readRDS("data/students.rds")
data

Nonparametrische Tests

\(\chi^2\)-Test

Wir stellen uns am Anfang eine ganz wichtige Frage! Ist die Ernährungsweise unabhängig vom Geschlecht? Ja, diese Frage ist wichtig. Das wichtige an ihr für uns ist vor allem, dass sie sich wunderbar für einen \(\chi^2\)-Test eignet.

result <- chisq.test(data$geschlecht, data$vegetarier)
result
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data$geschlecht and data$vegetarier
## X-squared = 3.9077, df = 1, p-value = 0.04807
result$observed
##                data$vegetarier
## data$geschlecht ja nein
##        männlich  2   21
##        weiblich 27   57

Wie zu erwarten: Männer wollen immer noch ihr Steak. Der Effekt ist gar nicht klein und auf einem \(\alpha\)-Niveau von 5% haben wir ein signifikantes Ergebnis. Aufpassen müssen wir, wenn wir das Ergebnis per Hand nachrechnen wollen, denn man sieht im Output, dass hier noch eine Korrektur verwendet wird, die zu einem genaueren \(p\)-Wert führt. Das Problem ist, dass sich dadurch der \(\chi^2\)-Wert verändert und der Effekt würde falsch berechnet werden:

sqrt(3.9077 / 107) # für Effektgröße w und bei 2x2-Tabelle ist das auch Phi
## [1] 0.1911035

Tatsächlich ist der Effekt:

cor(as.numeric(as.factor(data$geschlecht)), as.numeric(as.factor(data$vegetarier)))
## [1] -0.2166947

Wobei hier die spezifische Kodierung zu einem negativen Effekt führt, was für uns allerdings keine Rolle spielt, denn wir sehen aus der Tabelle, dass Frauen eher Vegetarier sind als Männer. Die komischen as.factor und as.numeric sind hier notwendig, denn das Geschlecht und die Ernährungsweise sind als Characters gespeichert.

Wie kommen wir über den \(\chi^2\)-Wert auf den richtigen Effekt? Einfach die Korrektur abschalten:

chisq.test(data$geschlecht, data$vegetarier, correct = FALSE)
## 
##  Pearson's Chi-squared test
## 
## data:  data$geschlecht and data$vegetarier
## X-squared = 5.0244, df = 1, p-value = 0.02499

Und nun passt die Formel:

sqrt(5.0244/107)
## [1] 0.2166957

Die Detektivarbeit war erfolgreich und wir sehen auch, dass die Korrelation äquivalent zur Effektgröße \(\Phi\) und \(w\) ist – zumindest für 2×2-Tabellen. Jetzt bist Du dran. In der Übung zur Methodenlehre haben wir aus einem Artikel folgende Tabelle extrahiert:

a <- c(19, 27, 9, 13, 8, 27)
b <- c(25, 26, 24, 15, 9, 4)
gutachter <- matrix(a+b, nrow=2, byrow=T)
colnames(gutachter) <- c("Plant Relocation", "Tax evasion", "Independence")
gutachter <- as.data.frame(gutachter, row.names = c("Recognize", "Do not recognize"))
gutachter

Es ging um die Frage ob Gutachter ethische Probleme in Finanz-Anträgen identifizieren können. Es gab drei Szenarien (in den Spalten) und als AV wurde erfasst, ob das Problem erkannt wurde oder nicht. Man kann schon nach einem kurzen Blick erkennen, dass bei Steuerhinterziehung die Erkennungsrate am höchsten ist. Die Frage ist generell, ob es Unterschiede in der Erkennungsrate zwischen den Szenarien gibt. Als \(\chi^2\)-Wert kam im Paper 8.6 heraus. Prüfe diesen nach:

chisq.test(gutachter)

\(U\)-Test

Wir haben an Tag 4 einen \(t\)-Test gerechnet um zu prüfen, ob sich die Körpergröße zwischen Männern und Frauen unterscheidet. Die Voraussetzungen für diesen Test sind vielleicht nicht erfüllt. Eine Alternative ist hier der \(U\)-Test, der in R jedoch unter dem Wilcoxon-Test läuft (es gibt viele Namen für diesen Test, auch noch Mann-Whitney-Wilcoxon, o. Wilcoxon-Mann-Whitney o. eben Wilcoxon rank-sum test):

m <- data[data$geschlecht == "männlich", "groesse"]
w <- data[data$geschlecht == "weiblich", "groesse"]
wilcox.test(m, w, paired = FALSE) # paired = FALSE ist der default, kann man also auch weglassen
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  m and w
## W = 1793.5, p-value = 3.403e-10
## alternative hypothesis: true location shift is not equal to 0

Im Output sieht man mal wieder, dass eine Korrektur angewandt wird. Wenn Du also was per Hand rechnest, vergiss nicht die Korrektur auszustellen. Übrigens ist \(W\) hier die Teststatistik, die wir eigentlich als \(U\) kennen.

Ist gar nicht schwer, oder? Das kannst Du auch. In der Übung haben wir uns ein Paper angeschaut, bei dem es um die Frage ging ob Frauen und Männer in unterschiedlichen Umwelten (friedlich oder kompetitiv) besser oder schlechter performen. Im Folgenden findest Du die Anzahl der gelösten Aufgaben:

solved_female <- c(24, 20, 18, 17, 15, 15, 14, 14, 14, 13, 13, 11, 11, 10, 9, 
                   9, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 4, 4)
solved_male <- c(22, 21, 21, 21, 19, 19, 17, 17, 17, 17, 17, 18, 16, 16, 16,
                 15, 15, 15, 15, 13, 13, 13, 12, 11, 11, 11, 10, 10, 8, 5)

Rechne den \(U\)-Test für den Vergleich zwischen Männern und Frauen. Stell hier die Korrektur aus, sodass der gleiche Wert herauskommt wie bei der Berechnung per Hand in der Übung.

wilcox.test(solved_female, solved_male, correct = FALSE)

Wilcoxon-Test (signed-rank)

Wir haben an Tag 4 auch einen \(t\)-Test gerechnet um zu prüfen, ob sich die Deutsch- und Mathepunkte unterscheiden. Die Voraussetzungen für diesen Test sind jedoch sehr wahrscheinlich nicht erfüllt. Eine Alternative ist hier der Wilcoxon-Test (signed-rank):

wilcox.test(data$mathe, data$deutsch, paired = TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  data$mathe and data$deutsch
## V = 1773, p-value = 0.202
## alternative hypothesis: true location shift is not equal to 0

Die Test-Statistik ist hier \(V\), wir kennen sie unter \(T\).

Die entscheidende Frage ist, wie weit der \(p\)-Wert dieses Tests vom \(p\)-Wert des \(t\)-Tests entfernt ist:

t.test(data$mathe, data$deutsch, paired = TRUE)
## 
##  Paired t-test
## 
## data:  data$mathe and data$deutsch
## t = -1.623, df = 106, p-value = 0.1076
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.9965795  0.0993832
## sample estimates:
## mean difference 
##      -0.4485981

Tatsächlich macht das hier einen großen Unterschied (0.202 versus 0.1076)! Das ist jedoch nicht immer so, man sollte sich aber auf jeden Fall Gedanken machen auf welchem Skalenniveau man misst und ob die Voraussetzungen für einen Test näherungsweise gegeben sind. Ist dies nicht der Fall, so ist man mit den nonparametrischen Tests auf der sicheren Seite. Wie schon erwähnt, die eigentlich bessere Alternative ist Bootstrapping und Randomisierungstests, aber die behandeln wir erst im Master.

Interessant wäre hier noch zu schauen, wie groß der Unterschied zu einem unabhängigen Wilcoxon-Test ist:

wilcox.test(data$mathe, data$deutsch, paired = FALSE)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data$mathe and data$deutsch
## W = 5396, p-value = 0.4648
## alternative hypothesis: true location shift is not equal to 0

Macht einen noch größeren Unterschied (0.202 versus 0.4648)! Das sollte klar sein, denn bei abhängigen Stichproben sind die Effekte immer größer und man sollte das hier nicht verschenken. Im umgekehrten Fall sollte man auch nie unabhängige Stichproben als abhängige behandeln; das wäre grob fahrlässig!

In der Übung hatten wir auch einen Wilcoxon-Test gerechnet. Es ging um eine neue Therapieform bei Anorexie-Patienten. Gemessen wurde das Körpergewicht in kg vor und nach der Therapie. Die Daten waren:

gewicht_t1 <- c(58, 55, 60, 52, 53, 53, 55, 49, 50, 55)
gewicht_t2 <- c(61, 60, 64, 56, 59, 60, 59, 47, 52, 56)

Prüfe nach ob Du auf das selbe Ergebnis kommst wie in der Übung (wir haben einseitig getestet, auf eine Erhöhung des Gewichts):

wilcox.test(gewicht_t2, gewicht_t1, paired = TRUE, alternative="greater", correct = FALSE)

Sehr gut! Jetzt haben wir aber wirklich alle Signifikanztests der Vorlesung behandelt. Ich hoffe ich habe nichts vergessen.

Das war harte Arbeit und jetzt kommt die Belohnung. Poweranalyse :). Klingt wohl eher nach noch mehr Arbeit, aber bevor Du zu nicht-freier, unmoralischer Software wie G*Power greifst, investiere lieber noch 30 Minuten in R!

Power-Analyse

Wir müssen zunächst ein Paket installieren. Ich benutze gerne WebPower, da es viele Fälle abdeckt und auch ein Online-Interface hat (https://webpower.psychstat.org/login.php). Installation:

install.packages("WebPower")

Nun laden wir das Paket:

library(WebPower)

Mit WebPower sind wir zu einer Supermacht aufgestiegen. Power für Korrelationen:

wp.correlation(r = .3, power = .8, alternative = "greater")
## Power for correlation
## 
##            n   r alpha power
##     66.55538 0.3  0.05   0.8
## 
## URL: http://psychstat.org/correlation
wp.correlation(r = .3, power = .8, alternative = "two.sided")
## Power for correlation
## 
##            n   r alpha power
##     83.94932 0.3  0.05   0.8
## 
## URL: http://psychstat.org/correlation

Wir geben also die Effektgröße, die gewünschte Power und die Alternativhypothese (einseitig, zweiseitig) an und bekommen die Stichprobengröße heraus. Nicht übel.

Diese Syntax funktioniert für alle möglichen Tests. Zum Beispiel \(t\)-Test:

wp.t(d = 0.5, power = .9, type = "two.sample")
## Two-sample t-test
## 
##            n   d alpha power
##     85.03126 0.5  0.05   0.9
## 
## NOTE: n is number in *each* group
## URL: http://psychstat.org/ttest

Bei der ANOVA müssen wir allerdings die Effektgröße, die uns bekannt ist (\(\eta^2\)) in \(f\) umwandeln. Ist ja aber in R ein Klacks:

eta_squared <- .1
f <- sqrt(eta_squared / (1 - eta_squared))
wp.anova(k = 3, power = .7, f = f)
## Power for One-way ANOVA
## 
##     k       n         f alpha power
##     3 72.3917 0.3333333  0.05   0.7
## 
## NOTE: n is the total sample size (overall)
## URL: http://psychstat.org/anova

Wobei hier \(k\) die Anzahl der Gruppen ist.

Okay, dann kannst Du jetzt mal ein bisschen üben:

In der Vorlesung zur ANOVA wurde eine Tabelle für die Power präsentiert:

n \(\eta^2_p = 0.01\) \(\eta^2_p = 0.06\) \(\eta^2_p = 0.14\)
10 0.09 0.34 0.69
20 0.14 0.60 0.94
30 0.19 0.78 0.99
40 0.25 0.88 *
50 0.29 0.94 *

Prüfe die Werte nach für ein \(2\times2\)-Design, ein \(\eta^2=0.06\), und die Stichprobengrößen (pro Gruppe) 30, 40 und 50 (In der Vorlesung sind diese als 0.78, 0.88 und 0.94 angegeben). Hinweis: Du musst in diesem Fall die Funktion wp.kanova benutzen. Lies in der Hilfe oder im Web nach, was die Parameter bedeuten.

eta_squared <- .06
f <- sqrt(eta_squared / (1 - eta_squared))
wp.kanova(n = 120, ndf = 1, f = f, ng = 4)
wp.kanova(n = 160, ndf = 1, f = f, ng = 4)
wp.kanova(n = 200, ndf = 1, f = f, ng = 4)

Du planst eine eigene Studie und gehst von einem Effekt von \(d=0.6\) aus, das \(\alpha\)-Niveau legst Du auf 0.10 fest. Du würdest gerne wissen wie groß die Power Deiner Studie in Abhängigkeit der Stichprobengröße ist. Dafür willst Du eine Abbildung erstellen; auf der \(x\)-Achse die Stichprobengröße und auf der \(y\)-Achse die Power. Die minimale Stichprobengröße soll 10 sein, die maximale 500 (in 10er-Schritten, schau Dir hierfür die Hilfs-Funktion seq an). Erstelle mit Hilfe der Funktion wp.t zunächst die Daten, greife dann auf die Stichprobengröße und Power zu und erstelle, die Abbildung.

power <- wp.t(d = 0.6, n1 = seq(10, 500, 10), type = "two.sample")
plot(power$n, power$power, type = "b")

Nicht schlecht! Spürst Du jetzt die Power der Power-Analyse?

Für heute bleiben uns noch zwei kleine Themen, die irgendwie nirgends richtig dazugehören: die Lowess-Kurve und die Partialkorrelation.

Lowess-Kurve

Wir schauen uns nochmal den Zusammenhang zwischen den Mathe- und Deutschpunkten an, diesmal jedoch mit einer Lowess-Kurve:

sunflowerplot(data$mathe, data$deutsch)
lines(lowess(x = data$mathe, y = data$deutsch))

Standardmäßig ist \(f=2/3\), je kleiner der Wert ist, desto zackiger (lokaler) wird das Smoothing:

sunflowerplot(data$mathe, data$deutsch)
lines(lowess(x = data$mathe, y = data$deutsch))
lines(lowess(x = data$mathe, y = data$deutsch, f = 0.2), col = "green")
lines(lowess(x = data$mathe, y = data$deutsch, f = 0.1), col = "blue")

Erstelle nun Du einen Plot für den Zusammenhang zwischen dem Mögen von Hunden (\(x\)-Achse) und dem Mögen von Katzen (\(y\)-Achse) mit einer Lowess-Kurve (\(f=2/3\)). Die Variablen heißen hunde_m u. katzen_m.

sunflowerplot(data$hunde_m, data$katzen_m)
lines(lowess(x = data$hunde_m, y = data$katzen_m))

Partialkorrelation

Hierzu gibt es nicht viel zu sagen. Wir verwenden einfach die Funktion pcor.test aus dem Paket ppcor:

install.packages("ppcor")

Wir laden das Paket:

library(ppcor)

Fürs IQ-Beispiel aus der Vorlesung:

verb_test  <- c(4, 5, 6, 7, 8, 10, 11, 12, 14, 15)
num_test <- c(5, 6.5, 7, 5.5, 9, 6.5, 8, 10.5, 8, 10)
iq <- c(98, 98, 103, 101, 100, 106, 111, 125, 120, 115)
pcor.test(verb_test, num_test, iq) 

Herauspartialisiert wird der IQ aus dem Zusammenhang zwischen einen verbalen und numerischem Test.

Probier Du es nun. Berechne die Partialkorrelation zwischen Deutsch- und Mathepunkten ohne den Einfluss der Abiturnote:

pcor.test(data$deutsch, data$mathe, data$abi)
# im Vergleich zu
cor.test(data$deutsch, data$mathe)

Super! Nun kannst Du entweder chillen oder Du festigst Dein Wissen mit ein paar weiteren Übungsaufgaben.

Übungsaufgaben Tag 5

Aufgabe \(\chi^2\)-Test

Die folgenden Daten sind aus einer Übungsaufgabe der Methodenlehre zum Simpsons Paradox.

a <- c(820, 80, 680, 20)
b <- c(20, 80, 100, 200)
simpsonsparadox <- matrix(a+b, nrow=2, byrow=T)
colnames(simpsonsparadox) <- c("Accept", "Reject")
simpsonsparadox <- as.data.frame(simpsonsparadox, row.names = c("Male", "Female"))
simpsonsparadox

Dargestellt sind hier die Häufigkeiten der Annahme und Ablehnung von Studienplatzbewerbungen aufgeteilt nach Geschlecht. Wir haben in der Übung den \(p\)-Wert aus dem Artikel von Kievit hierzu nachgerechnet. Das haben wir per Hand gemacht, also ohne Yates Korrektur und kamen auf einen Wert von 11.696. Welcher Wert kommt in R (mit Korrektur) heraus?

chisq.test(simpsonsparadox, correct = T)

Das ist auch der Wert, der im Paper angegeben ist.

Aufgabe Partialkorrelation

Schreibe eine Funktion mit der Du aus Korrelationen Partialkorrelationen berechnen kannst. Benutze die Funktion anschließend für das Beispiel aus der Vorlesung:

vl <- matrix(c(1, .43, -.51, -.54, .43, 1, .41, -.12, -.51, .41, 1, .23, -.54, -.12, .23, 1), ncol = 4)
row.names(vl) <- c("Gewichtsverlust", "Sport", "Kalorienaufnahme", "Alter")
colnames(vl) <- c("Gewichtsverlust", "Sport", "Kalorienaufnahme", "Alter")
vl 
##                  Gewichtsverlust Sport Kalorienaufnahme Alter
## Gewichtsverlust             1.00  0.43            -0.51 -0.54
## Sport                       0.43  1.00             0.41 -0.12
## Kalorienaufnahme           -0.51  0.41             1.00  0.23
## Alter                      -0.54 -0.12             0.23  1.00

Der Zusammenhang zwischen Gewichtsverlust und Sport, ohne den Einfluss von Kalorienaufnahme soll berechnet werden. Achtung: Das ist für Anfänger eine recht anspruchsvolle Aufgabe, für die Du etwas Zeit einplanen solltest.

parcor <- function(xy, xz, yz){
  (xy - xz * yz) / (sqrt(1 - xz^2) * sqrt(1 - yz^2))
}
parcor(vl[1, 2], vl[1, 3], vl[2, 3])

Aufgabe Power

Wie groß ist die Power bei einer Regression mit zwei Prädiktoren, einer Stichprobengröße von 88, einem \(\alpha\) von 5% und einem \(R^2\) von 0.13? Benutze die Funktion wp.regression und finde heraus, wie Du \(R^2\) in \(f^2\) umrechnen kannst.

f2 = 0.13 / (1 - 0.13)
wp.regression(n = 88, p1 = 2, f2 = f2, alpha = .05)

Gut gemacht. Du hast diesen Tag erfolgreich überstanden. Nur noch zwei Tage, dann bist Du ein R-Profi!

5. Nonparametrische