Kapitel 12 Deskriptive Statistik

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.

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

library(readr)
qmsurvey  <- readRDS("data/qm_survey_ss2017.rds")
gotdeaths <- read_csv("data/got_deaths.csv", col_types = cols())

12.1.1 Mittelwert, Median, Modus

Mittelwert und Median sind einfach:

mean(qmsurvey$zufrieden)
[1] 3.539216
median(qmsurvey$zufrieden)
[1] 4

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:

table(qmsurvey$zufrieden)

 1  2  3  4  5 
 2 13 25 52 10 

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:

library(tadaatoolbox)

modus(qmsurvey$zufrieden)
[1] "4"

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:

library(sjmisc)

frq(qmsurvey$zufrieden, out = "viewer")
Wie zufrieden bist du insgesamt mit dem Psychologiestudium (x) <integer>
val frq raw.prc valid.prc cum.prc
1 2 1.96 1.96 1.96
2 13 12.75 12.75 14.71
3 25 24.51 24.51 39.22
4 52 50.98 50.98 90.2
5 10 9.8 9.8 100
NA 0 0 NA NA
total N=102 · valid N=102 · x̄=3.54 · σ=0.91

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.

12.1.2 Variabilität

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

var(qmsurvey$zufrieden)
[1] 0.8251796
sd(qmsurvey$zufrieden)
[1] 0.908394

Wenig überraschend. Solange ihr euch merken könnt, dass Standardabweichung auf englisch standard deviation 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:

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!

variation(qmsurvey$zufrieden)
[1] 83.34314

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.

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

median(qmsurvey$zufrieden)
[1] 4

Wenig beeindruckend.
Für Quantile:

quantile(qmsurvey$zufrieden)
  0%  25%  50%  75% 100% 
   1    3    4    4    5 

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

quantile(qmsurvey$alter, probs = c(.333, .666))
33.3% 66.6% 
   19    21 

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

quantile(qmsurvey$alter, probs = seq(from = 0, to = 1, by = 0.1))
  0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
18.0 19.0 19.0 19.0 20.0 20.0 21.0 22.0 24.8 27.0 37.0 

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:

min(qmsurvey$alter)
[1] 18
max(qmsurvey$alter)
[1] 37
range(qmsurvey$alter)
[1] 18 37
max(qmsurvey$alter) - min(qmsurvey$alter)
[1] 19

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

12.2 Bivariate Statistiken

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

library(sjPlot)

sjt.xtab(qmsurvey$cannabis, qmsurvey$behandlung,
         show.exp = TRUE, show.legend = TRUE)
Wie oft konsumierst
Du Cannabis
Warst du schon
einmal in
psychotherapeutischer
Behandlung
Total
Ja Nein
Nie 26
22
44
48
70
70
Selten 4
8
22
18
26
26
Regelmäßig 2
2
4
4
6
6
Total 32
32
70
70
102
102
χ2=4.180 · df=2 · Cramer’s V=0.202 · Fisher’s p=0.094

observed values
expected values

Die Berechnung in R ist einfach:

chisq.test(qmsurvey$cannabis, qmsurvey$behandlung)

    Pearson's Chi-squared test

data:  qmsurvey$cannabis and qmsurvey$behandlung
X-squared = 4.1801, df = 2, p-value = 0.1237
# Reihenfolge egal
chisq.test(qmsurvey$behandlung, qmsurvey$cannabis)

    Pearson's Chi-squared test

data:  qmsurvey$behandlung and qmsurvey$cannabis
X-squared = 4.1801, df = 2, p-value = 0.1237

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

library(tadaatoolbox)

nom_chisqu(qmsurvey$behandlung, qmsurvey$cannabis)
[1] 4.180078

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:

nom_phi(qmsurvey$behandlung, qmsurvey$cannabis)
[1] 0.202438
nom_v(qmsurvey$behandlung, qmsurvey$cannabis)
[1] 0.202438

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:

nom_c(qmsurvey$behandlung, qmsurvey$cannabis)
[1] 0.1984133

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
nom_lambda(qmsurvey$behandlung, qmsurvey$cannabis)
[1] 0
nom_lambda(qmsurvey$behandlung, qmsurvey$cannabis, reverse = TRUE)
[1] 0
nom_lambda(qmsurvey$behandlung, qmsurvey$cannabis, symmetric = TRUE)
[1] 0

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:

library(vcd)

freq_table <- table(qmsurvey$behandlung, qmsurvey$cannabis)
freq_table
      
       Nie Selten Regelmäßig
  Ja    26      4          2
  Nein  44     22          4
assocstats(freq_table)
                    X^2 df P(> X^2)
Likelihood Ratio 4.5754  2  0.10150
Pearson          4.1801  2  0.12368

Phi-Coefficient   : NA 
Contingency Coeff.: 0.198 
Cramer's V        : 0.202 

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

tadaa_nom(qmsurvey$behandlung, qmsurvey$cannabis, print = "markdown")
\(\chi^2\) Cramer’s V \(\lambda_x\) \(\lambda_y\) \(\lambda_{xy}\) c
4.18 0.2 0 0 0 0.2



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.

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

ord_gamma(qmsurvey$lernen, qmsurvey$zufrieden)
[1] 0.1170716
ord_somers_d(qmsurvey$lernen, qmsurvey$zufrieden)
[1] 0.1090856
ord_somers_d(qmsurvey$lernen, qmsurvey$zufrieden, reverse = TRUE)
[1] 0.07730777
ord_somers_d(qmsurvey$lernen, qmsurvey$zufrieden, symmetric = TRUE)
[1] 0.09319667

Wobei Somers’ D analog Lambda drei Geschmacksrichtungen hat.

Oder auch hier, alle zusammen:

tadaa_ord(qmsurvey$lernen, qmsurvey$zufrieden, print = "markdown")
\(\gamma\) \(D_x\) \(D_y\) \(D_{xy}\) \(\tau_A\) \(\tau_B\) \(\tau_C\)
0.12 0.11 0.08 0.09 0.07 0.09 0.09



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

ord_tau(qmsurvey$lernen, qmsurvey$zufrieden, tau = "a")
[1] 0.07202485
ord_tau(qmsurvey$lernen, qmsurvey$zufrieden, tau = "b")
[1] 0.09183225
ord_tau(qmsurvey$lernen, qmsurvey$zufrieden, tau = "c")
[1] 0.0891484

Da \(\tau_b\) auch als Korrelationskoeffizient durchgeht, gibts’s den auch via cor:

cor(qmsurvey$lernen, qmsurvey$zufrieden, method = "kendall")
[1] 0.09183225

12.2.3 Intervallskala

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

cov(qmsurvey$alter, qmsurvey$beziehungen, use = "complete.obs")
[1] 2.810619

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"

cor(qmsurvey$alter, qmsurvey$beziehungen, use = "complete.obs")
[1] 0.5225806
cor(qmsurvey$alter, qmsurvey$beziehungen, use = "complete.obs", method = "spearman")
[1] 0.5303142

Über method = "kendall" bekommen wir auch Kendall’s \(\tau_b\). Just in case.

12.2.4 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

12.2.4.1 \(\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 Anwendung, weshalb auch R eher diesen Kontext erwartet.

Wenn ihr \(\eta^2\) aber “einfach so” haben wollt, braucht ihr etwas Umweg:

library(sjstats)

anova <- aov(alter ~ alkohol, data = qmsurvey)
eta_sq(anova)
     term etasq
1 alkohol  0.09

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:

library(tadaatoolbox)

etasq(alter ~ alkohol, data = qmsurvey)
[1] 0.08971183

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

12.2.4.2 Cohen’s Kappa (\(\kappa\))

library(irr)
data(anxiety) # Beispieldatensatz aus irr package laden
head(anxiety)
  rater1 rater2 rater3
1      3      3      2
2      3      6      1
3      3      4      4
4      4      6      4
5      5      2      3
6      5      4      2
# Nur zwei Spalten / Rater für Kappa
anxiety_2 <- anxiety[c(1, 2)]

kappa2(anxiety_2)
 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 20 
   Raters = 2 
    Kappa = 0.119 

        z = 1.16 
  p-value = 0.245 

12.2.4.3 Kendall’s \(W\)

library(irr)

kendall(anxiety)
 Kendall's coefficient of concordance W

 Subjects = 20 
   Raters = 3 
        W = 0.502 

Chisq(19) = 28.6 
  p-value = 0.0724 

 Coefficient may be incorrect due to ties