Modul 4 - Bivariate deskriptive Statistik

Sebastian Schnettler und Svenja Heinecke 2024

 

Bivariate statistische Maße in R

Im letzten Modul zur univariaten deskriptiven Statistik haben wir bereits einige statistische Funktionen behandelt. Das Prinzip ist immer dasselbe: Das statistische Maß wird durch eine Funktion in R erzeugt, die auf ein R-Objekt angewendet werden kann. Dieses Objekt ist in der Regel eine Variable in einem Datensatz. In ähnlicher Weise gibt es eine Reihe von Funktionen für die bivariate deskriptive Statistik.

Kreuztabellen mit absoluten und relativen Häufigkeitstabellen in R

In diesem Abschnitt werden wir einige dieser Funktionen zunächst anhand eines fiktiven Beispieldatensatzes veranschaulichen, den wir hier spontan mit der Funktion data.frame() erzeugen. Der Datensatz wird aus zwei Variablen mit insgesamt 25 Fällen zur sozialen Mobilität bestehen. Die beiden Variablen enthalten zum einen die Information über die soziale Schicht der Herkunftsfamilie, zum anderen die Angabe über die eigene soziale Schichtzugehörigkeit im Erwachsenenalter. Die Merkmalsausprägungen sind jeweils “Unterschicht”,“Mittelschicht” und “Oberschicht”.

mydata <- data.frame(Herkunftsschicht =   
                       c("2 Mittel","3 Ober", "3 Ober","1 Unter","2 Mittel", 
                         "3 Ober","1 Unter","1 Unter", "1 Unter","1 Unter",
                         "1 Unter","2 Mittel","3 Ober","1 Unter","1 Unter", 
                         "1 Unter","1 Unter","2 Mittel", "3 Ober","3 Ober",
                         "3 Ober", "2 Mittel","2 Mittel","3 Ober","2 Mittel"),
                     Eigene.Schichtzugehörigkeit = 
                       c("3 Ober","3 Ober","2 Mittel","1 Unter","2 Mittel", 
                         "3 Ober","1 Unter","2 Mittel","1 Unter","1 Unter", 
                         "2 Mittel", "1 Unter","1 Unter","1 Unter","1 Unter", 
                         "3 Ober","2 Mittel", "2 Mittel","2 Mittel","3 Ober",
                         "3 Ober", "2 Mittel","2 Mittel","3 Ober","3 Ober"))

Wir können einen schnellen Blick auf die Verteilung der absoluten Häufigkeiten der Fälle werfen, indem wir die beiden Variablen mit dem xtabs()-Befehl in R kreuztabellieren. Dieser Befehl ist vor allem bei Tabellen nützlich, die, wie hier, dieselbe Anzahl an Zeilen und Spalten haben, da die Variablennamen als Zeilen- bzw. Spaltenbeschriftung mit anzeigt werden:

xtabs(~ Herkunftsschicht + Eigene.Schichtzugehörigkeit, data = mydata)

##                 Eigene.Schichtzugehörigkeit
## Herkunftsschicht 1 Unter 2 Mittel 3 Ober
##         1 Unter        6        3      1
##         2 Mittel       1        4      2
##         3 Ober         1        2      5

Auch kann der bereits in Modul 2 vorgestellte table()- Befehl verwendet werden, um die absoluten Häufigkeiten anzeigen zu lassen. Hierzu schreiben wir, mit Komma getrennt, statt nur einer Variable eine zweite Variable mit in den Befehl. Mit dem zusätzlichen Argument dnn = lässt sich außerdem eine eigene Zeilen- und Spaltenbeschriftung wählen:

table(mydata$Herkunftsschicht, mydata$Eigene.Schichtzugehörigkeit, 
      dnn = c("Herkunftsschicht", "Eigene Schichtzugehörigkeit"))

##                 Eigene Schichtzugehörigkeit
## Herkunftsschicht 1 Unter 2 Mittel 3 Ober
##         1 Unter        6        3      1
##         2 Mittel       1        4      2
##         3 Ober         1        2      5

Ein einfacher Blick auf die Tabelle zeigt eine starke Konzentration von Fällen in der Diagonalen, was deutlich macht, dass viele Personen in ihrer sozialen Herkunftsschicht verbleiben. Es gibt aber auch Fälle außerhalb der Diagonalen, die zeigen, dass ein gewisses Maß an sozialer Mobilität besteht.

Um die relativen Häufigkeiten in Prozentwerten zu ermitteln, kann die Funktion prop.table() auf die Tabelle angewendet werden. Auch diesen Befehl hatten wir bereits in Modul 2 im Kapitel “Haufigkeiten” gezeigt. Zu diesem Zweck speichern wir die Tabelle zunächst als R-Objekt my.table und wenden dann die Funktion prop.table() auf sie an.

my.table <- table(mydata$Herkunftsschicht, mydata$Eigene.Schichtzugehörigkeit, 
      dnn = c("Herkunftsschicht", "Eigene Schichtzugehörigkeit"))
prop.table(my.table)

##                 Eigene Schichtzugehörigkeit
## Herkunftsschicht 1 Unter 2 Mittel 3 Ober
##         1 Unter     0.24     0.12   0.04
##         2 Mittel    0.04     0.16   0.08
##         3 Ober      0.04     0.08   0.20

Wenn Sie die Funktion mit den Standardeinstellungen ausführen, erhalten Sie Prozentsätze auf der Grundlage der Gesamtzahl von 25 Fällen, die zu 1 bzw. zu 100% in der Tabelle aufsummiert sind. So lässt sich hier bspw. aussagen, dass ein Anteil von 0,04 bzw. 4% aller Fälle als Information zur Herkunftsschicht die Angabe “Mittelschicht” hat und als eigene Schichtzugehörigkeit die Information “Unterschicht”. 0,2 bzw. 20% haben als Information zur Herkunftsfamilie “Oberschicht” und als Information zur eigenen Schichtzugehörigkeit ebenfalls “Oberschicht”.

Für das Lernen über Assoziationen ist dies jedoch noch nicht sehr hilfreich. Stattdessen ist es informativer, sich die bedingten relativen Häufigkeiten in Zeilen- und Spaltenprozenten anzusehen. Diese können durch Setzen des Arguments margin ermittelt werden: margin = 1 liefert dabei die Zeilen- und margin = 2 die Spaltenprozente.

prop.table(my.table, margin = 1)

##                 Eigene Schichtzugehörigkeit
## Herkunftsschicht   1 Unter  2 Mittel    3 Ober
##         1 Unter  0.6000000 0.3000000 0.1000000
##         2 Mittel 0.1428571 0.5714286 0.2857143
##         3 Ober   0.1250000 0.2500000 0.6250000

Umklammern wir dies zusätzlich mit dem round()-Befehl, lässt sich die Darstellung der Nachkommastellen begrenzen und die Tabelle wird noch übersichtlicher:

round(prop.table(my.table, margin = 1),2)

##                 Eigene Schichtzugehörigkeit
## Herkunftsschicht 1 Unter 2 Mittel 3 Ober
##         1 Unter     0.60     0.30   0.10
##         2 Mittel    0.14     0.57   0.29
##         3 Ober      0.12     0.25   0.62

Zu sehen sind hier mit margin = 1 zunächst die relativen Häufigkeiten in Zeilenprozenten, die die Anteile für die eigene Schichtzugehörigkeit bedingt für die Herkunftsschicht angeben. So lässt sich erkennen, dass ein Anteil von 0,6 bzw. 60% derjenigen Personen, deren Herkunftsschicht die “Unterschicht” abbildet, auch in der “Unterschicht” verbleiben. Ein Anteil von 0,3 bzw. 30% erleben einen sozialen Aufstieg in die “Mittelschicht”, während lediglich ein Anteil von 0,1 bzw. 10% einen Aufstieg in die “Oberschicht” schaffen. Die Tabelle ist hier nun also zeilenweise zu je 1 bzw. 100% aufsummiert.

round(prop.table(my.table, margin = 2),2)

##                 Eigene Schichtzugehörigkeit
## Herkunftsschicht 1 Unter 2 Mittel 3 Ober
##         1 Unter     0.75     0.33   0.12
##         2 Mittel    0.12     0.44   0.25
##         3 Ober      0.12     0.22   0.62

Mit margin = 2 sind hier nun die relativen Häufigkeiten in Spaltenprozenten abgetragen, die die Anteile für die Herkunftsschicht bedingt für die eigene Schichtzugehörigkeit angeben. Zu sehen ist bspw., dass von allen Personen, die als Information für die eigene Schichtzugehörigkeit die Ausprägung “Oberschicht” haben, ein Anteil von 0,62 bzw. 62% ebenfalls “Oberschicht” als Angabe ihrer Herkunftsfamilie hat. Ein Anteil von 0,25 bzw. 25% stammen ursprünglich aus der “Mittelschicht” und ein Anteil von 0,12 bzw. 12% aus der “Unterschicht”. Die Tabelle ist hier nun also spaltenweise zu je 1 bzw. 100% aufsummiert. (Dass sich dies nicht genau zu 1 aufsummiert, sondern zu 0,99, liegt hier an Rundungsdifferenzen.)

Wie gezeigt, unterscheiden sich die Werte innerhalb derselben Zeile und zwischen den Spalten recht deutlich voneinander (Zeilenprozente). Vice versa unterscheiden sich auch die Werte innerhalb derselben Spalte und zwischen den Zeilen (Spaltenprozente). Dies lässt den Schluss zu, dass es einen Zusammenhang zwischen den beiden Variablen gibt.

Erweiterungen für schönere Tabellen

Die in R eingebaute Tabellenfunktion erzeugt sehr einfache Tabellen. Wenn Sie schönere, vielleicht sogar publikationsreife Tabellen erstellen möchten, sollten Sie sich die Funktionen ansehen, die einige externe Zusatzpakete anbieten:

• das Paket {gt}, https://gt.rstudio.com/ (und eine Erweiterung namens {gtextras}, mit der Sie sogar Datenvisualisierungen in Ihre Tabellen einbetten können)
• die Funktion kable im Paket {knitr}: https://bookdown.org/yihui/rmarkdown-cookbook/kable.html
• das formattable-Paket: https://renkun-ken.github.io/formattable/
• das flextable-Paket: https://davidgohel.github.io/flextable/index.html
• letzte Abrufe: März 2024

Zusammenhangsmaße für nominal-skalierte Variablen

Um den vermuteten Zusammenhang zwischen den beiden Variablen für die Schichtzugehörigkeit aus unserem fiktiven Beispiel-Datensatz zu testen, lassen sich Zusammenhangsmaße berechnen. Für die Schichtzugehörigkeiten nehmen wir an dieser Stelle einmal an, dass es sich um nominal-skalierte Variablen handelt, wobei eine Rangabfolge der Werte (für ein ordinales Skalenniveau) damit ausgeschlossen wird.

“Cramer’s V” und “Phi” (als Spezialfall von Cramer´s V für 2x2-Tabllen) sind Beispiele für Zusammenhangsmaße für nominal-skalierte Variablen. Die Berechnungen dazu können über das Zusatzpaket DescTools bezogen werden. (Insofern das Paket noch nicht installiert wurde, muss dieses mit dem install.packages("DescTools")-Befehl zunächst heruntergeladen und installiert werden.) Angewandt auf die Tabelle von oben, können wir “Cramer’s V” direkt auf unser R-Objekt my.table anwenden und berechnen, oder wir können die Variablennamen unseres Datensatzes mydata für die Argumente x und y verwenden:

library(DescTools) # Paket "DescTools" aufrufen
CramerV(my.table) 

## [1] 0.4286954

# ODER:
CramerV(x = mydata$Herkunftsschicht, y = mydata$Eigene.Schichtzugehörigkeit)

## [1] 0.4286954

Das zwischen 0 und 1 normierte Zusammenhangsmaß “Cramer´s V” zeigt hier mit einem Wert von gerundet 0,43 eine mittlere Zusammenhangsstärke zwischen der Schichtzugehörigkeit der Herkunftsfamilie und der eigenen späteren Schichtzugehörigkeit nach Verlassen des Elternhauses an. Getestet wird der Zusammenhang allerdings ungerichtet, was bedeutet, dass wir nicht gestestet haben, ob eine der Variablen einen Einfluss auf die andere hat, auch wenn wir aufgrund der Zeitlichkeit hier annehmen können, dass die Herkunftsschicht einen Einfluss auf die Schichtzugehörigkeit im Erwachsenenalter hat und nicht umgekehrt. Wir können also nur aussagen, dass ein grundsätzlicher Zusammenhang zwischen den beiden Variablen gefunden wurde.

Zusammenhangsmaße für ordinal-skalierte Variablen

Geläufige und häufig genutzte Korrelationsmaße für ordinale Variablen sind:

• der Rangkorrelationskoeffizient nach Spearman (auch: Spearman´s Rho),
• das Konkordanzmaß Kendall’s Tau,
• das Konkordanzmaß Goodman und Kruskal´s Gamma.

Um die aufgeführten Maße beispielhaft zu berechnen, nutzen wir den Übungsdatensatz des Sozioökonomischen Panels (SOEP) der Jahre 2015-2019, den wir bereits kennengelernt haben. Dieser kann, insofern noch nicht geschehen, hier als practice_dataset.dta heruntergeladen werden:

https://www.diw.de/de/diw_01.c.603154.de/soep_in_der_lehre.html#c_603164

Sobald Sie die Daten heruntergeladen und in ihrem Modul-Ordner abgespeichert haben, gibt es wieder zwei Möglichkeiten, den Pfad zum Arbeitsverzeichnis für R zu setzen, in dem ihre Materialen liegen. Diese Möglichkeiten hatten wir Ihnen im Modul 2 Teil 1 ausführlich vorgestellt. Sie setzen entweder den Pfad mit dem setwd()-Befehl oder Sie haben ein Projekt in R angelegt und (wieder) geöffnet, womit der Pfad zum Arbeitsverzeichnis bereits intern gesetzt ist. Danach können Sie die SOEP-Daten mit der Funktion read_stata() aus dem Paket haven importieren. Wir nutzen hier wieder den Pfad mit dem vorgeschlagenen Unterordner “Daten”.

library(haven)
d1 <- read_stata("./Daten/practice_dataset.dta")

Der Befehl head() zeigt uns dann die ersten Zeilen des Datensatzes an. Mit dem Befehl View() ließe sich die komplette Datenbank anzeigen (hier nicht gezeigt).

head(d1)

## # A tibble: 6 x 15
##       id syear sex             alter anz_pers anz_kind bildung  erwerb branche
##    <dbl> <dbl> <dbl+lbl>       <dbl>    <dbl>    <dbl> <dbl+lb> <dbl+l> <dbl+lb>
## 1    194  2015 1 [[1] weiblich]   59        2        0 10.5     2 [[2]~ 84 [[84~
## 2    194  2016 1 [[1] weiblich]   60        2        0 10.5     2 [[2]~ 84 [[84~
## 3    194  2017 1 [[1] weiblich]   61        2        0 10.5     2 [[2]~ 84 [[84~
## 4    194  2018 1 [[1] weiblich]   62        2        0 10.5     2 [[2]~ 84 [[84~
## 5    194  2019 1 [[1] weiblich]   63        2        0 10.5     2 [[2]~ 84 [[84~
## 6  19052  2015 1 [[1] weiblich]   74        2        0 10 5       [[5]~ NA
## #  i 6 more variables: gesund_org <dbl+lbl>, lebensz_org <dbl+lbl>,
## #    einkommenj1 <dbl>, einkommenj2 <dbl>, einkommenm1 <dbl>, einkommenm2 <dbl>

Im Folgenden wird uns eine mögliche Korrelation zwischen den Variablen gesund_org für die eigene, subjektiv empfundene Gesundheit und lebensz_org für die eigene, subjektiv empfundene Lebenszufriedenheit interessieren. Mit dem table()-Befehl können wir uns zunächst einen Überblick über die Wertebereiche der Variablen verschaffen, wobei wir mit dem Zusatz exclude = NULL erzwingen, dass uns auch ggf. enthaltene Werte mit NA angezeigt werden, die einen fehlenden Wert markieren.

table(d1$lebensz_org, exclude = NULL)

## 
##    0    1    2    3    4    5    6    7    8    9   10 <NA> 
##   62   83  202  418  574 1837 1868 4645 7963 3909 1350  611

Der Wertebereich der Variablen für die Lebenszufriedenheit einer Person liegt zwischen “0 - ganz und gar unzufrieden” und “10 - ganz und gar zufrieden”. Außerdem sind 102 fehlende NA-Werte enthalten.

table(d1$gesund_org, exclude = NULL)

## 
##    1    2    3    4    5 <NA> 
## 2535 9488 7576 3089  732  102

Die Merkmalsausprägungen für die Werte der Variablen für das subjektive Gesundheitsempfinden einer Person geben Folgendes an: 1 = “Sehr gut”; 2 = “Gut”; 3 = “Zufriedenstellend”; 4 = “Weniger gut”; 5 = “Schlecht”. Es sind außerdem 102 fehlende NA-Werte enthalten.

Da für beide Variablen die subjektiven Einschätzungen der Teilnehmenden ausschlaggebend sind, kann nicht davon ausgegangen werden, dass zwischen den einzelnen Ausprägungen feste Abstände herrschen; daher sollten die Variablen nicht als metrisch, sondern als lediglich ordinal-skaliert behandelt werden.

Eine Kreuztabelle erzeugen wir dann mit dem xtabs()-Befehl, wobei die Anwendung des table()-Befehls (siehe oben, unter “Zusammenhangsmaße für nominal-skalierte Variablen”) ebenso möglich ist:

xtabs(~ lebensz_org + gesund_org , data = d1)

##            gesund_org
## lebensz_org    1    2    3    4    5
##          0     3    7    8   16   28
##          1     4   10   10   20   38
##          2     9   25   28   69   71
##          3     6   53   98  173   86
##          4    11   89  196  205   73
##          5    44  345  738  539  168
##          6    51  447  863  436   70
##          7   247 1632 2031  649   83
##          8   736 3955 2495  683   85
##          9   786 2122  792  194   13
##          10  446  566  248   75   13

Schaut man sich die Tabelle an, lässt sich bereits erahnen, dass eine höhere Lebenszufriedenheit mit geringeren Werten auf der Skala für die Gesundheit einer Person einhergeht. Geringere Werte bedeuten hier - etwas kontraintuitiv -, dass die Person eine bessere Einschätzung ihrer Gesundheit angegeben hat.

Da die Variable lebensz_org recht viele Ausprägungen hat, lohnt es sich der Anschaulichkeit halber an dieser Stelle einige der Ausprägungen jeweils in Klassen bzw. Kategorien zusammenzufassen, um einen deutlicheren Überblick zu erhalten. Dafür kann der bereits vorgestellte mutate()-Befehl verwendet werden, wobei in Verbindung mit der Funktion case_when() Bedingungen darüber formuliert werden können, welche Werte zu welcher Klasse zusammengefasst werden sollen. Beide Befehle sind im Zusatzpaket dplyr enthalten. (Sollte dieses noch nicht installiert sein, ist dies an dieser Stelle mit dem Befehl install.packages("dplyr") nachzuholen.) Die Werte werden im Folgenden so eingeteilt, dass etwa gleich viele Ausprägungen in einer Kategorie zusammengefasst werden, ohne dass dies unter einem vertieften Wissen zu der Forschung zu diesem Thema geschieht. Hier dient es lediglich dem Zweck der Veranschaulichung.

library(dplyr)
d1 <- mutate(d1,lebensz_kat = 
               case_when(d1$lebensz_org >= 0 & d1$lebensz_org <= 3 ~ "1 niedrig",
                         d1$lebensz_org >= 4 & d1$lebensz_org <= 6 ~ "2 mittel",
                         d1$lebensz_org >= 7 & d1$lebensz_org <= 10 ~ "3 hoch"))

Nun können wir die Kreuztabelle mit der neu erstellten Variable lebensz_kat betrachten:

xtabs(~ lebensz_kat + gesund_org, data = d1)

##            gesund_org
## lebensz_kat    1    2    3    4    5
##   1 niedrig   22   95  144  278  223
##   2 mittel   106  881 1797 1180  311
##   3 hoch    2215 8275 5566 1601  194

Neben dieser Kontingenztabelle mit absoluten Häufigkeiten lassen sich außerdem die Tabellen mit den relativen und den bedingten relativen Häufigkeiten in Zeilen- und Spaltenprozenten anzeigen. Der Einfachheit halber speichern wir die Tabelle zunächst wieder als R-Objekt ab.

t1 <- xtabs(~ lebensz_kat + gesund_org, data = d1)

Relative Häufigkeiten mit prop.table(), wobei round( ,2) die Darstellung auf zwei Nachkommastellen begrenzt:

round(prop.table(t1),2)

##            gesund_org
## lebensz_kat    1    2    3    4    5
##   1 niedrig 0.00 0.00 0.01 0.01 0.01
##   2 mittel  0.00 0.04 0.08 0.05 0.01
##   3 hoch    0.10 0.36 0.24 0.07 0.01

Bedingte relative Häufigkeiten in Zeilenprozenten mit margin = 1:

round(prop.table(t1, margin = 1),2)

##            gesund_org
## lebensz_kat    1    2    3    4    5
##   1 niedrig 0.03 0.12 0.19 0.36 0.29
##   2 mittel  0.02 0.21 0.42 0.28 0.07
##   3 hoch    0.12 0.46 0.31 0.09 0.01

Bedingte relative Häufigkeiten in Spaltenprozenten mit margin = 2:

round(prop.table(t1, margin = 2),2)

##            gesund_org
## lebensz_kat    1    2    3    4    5
##   1 niedrig 0.01 0.01 0.02 0.09 0.31
##   2 mittel  0.05 0.10 0.24 0.39 0.43
##   3 hoch    0.95 0.89 0.74 0.52 0.27

Wie Sie sehen können, unterscheiden sich die Werte innerhalb derselben Zeile und zwischen den Spalten (Zeilenprozente) z.T. recht deutlich voneinander. Vice versa unterscheiden sich auch die Werte innerhalb derselben Spalte und zwischen den Zeilen (Spaltenprozente). Dies lässt den Schluss zu, dass es einen Zusammenhang zwischen den beiden Variablen gibt. Wie bereits oben, deutet sich hier an, dass eine höhere Lebenszufriedenheit mit einem besseren gesundheitlichen Befinden einhergeht.

Schauen wir uns dazu die oben angegebenen Korrelationsmaße näher an, wobei wir hier wieder auf die Original-Variablen zurückgreifen. Die Funktion cor() aus den Grundbefehlen von R ermöglicht bereits die Berechnung einiger dieser Assoziationsmaße. Welcher Koeffizient berechnet wird, hängt dann von dem Namen ab, den wir in das Argument method = ... schreiben, das wir innerhalb der Funktion angeben können. Da die Variable gesund_org auch fehlende Werte enthält, fügen wir hier außerdem das Argument use = "pairwise.complete.obs" mit ein. Anderenfalls würde ein NA als Wert zurückgegegeben werden. Dieses Argument schließt auf die gleiche Weise wie das bereits vorgestellte na.rm = T die fehlenden Werte innerhalb der Variablen aus der Berechnung aus.

• Rangkorrelationskoeffizient nach Spearman:

cor(x = d1$lebensz_org, y = d1$gesund_org,
    use = "pairwise.complete.obs", method = "spearman")

## [1] -0.4306209

• Kendall’s Tau:

cor(x = d1$lebensz_org, y = d1$gesund_org,
    use = "pairwise.complete.obs", method = "kendall")

## [1] -0.3669886

Laut der Hilfeseite ?cor wird bei der Anwendung der Methode method = "kendall" im Falle des Vorhandenseins von Bindungen automatisch “Kendall’s Tau b” berechnet, ansonsten würde “Kendall´s Tau a” berechnet werden.

Weitere Korrelationsmaße wie “Goodman und Kruskal´s Gamma” und spezifischere Berechnungen von “Kendall´s Tau” finden sich wiederum im Zusatzpaket DescTools, das wir bereits oben für die Berechnung von “Cramer´s V” verwendet haben. Die Befehle innerhalb dieses Paketes kontrollieren selbstständig auf fehlende Werte und lassen diese bei Vorhandensein automatisch in der Berechnung außen vor. Auch hier ließe sich ein Tabellenobjekt, wie oben für “Cramer´s V” gezeigt, einfügen, auf dessen Darstellung an dieser Stelle allerdings verzichtet wird.

• Goodman und Kruskal´s Gamma:

library(DescTools) # falls oben noch nicht aufgerufen
GoodmanKruskalGamma(x = d1$lebensz_org, y = d1$gesund_org)

## [1] -0.4831976

• Kendall´s Tau a und b

KendallTauA(x = d1$lebensz_org, y = d1$gesund_org)

## [1] -0.273096

KendallTauB(x = d1$lebensz_org, y = d1$gesund_org)

## [1] -0.3669886

Wie zu sehen ist, unterscheiden sich hier die Werte für Kendall´s Tau a und b (und Goodman und Kruskal´s Gamma) deutlich. Dies lässt schlussfolgern, dass einige Bindungen innerhalb der Variablen-Werte vorhanden sind.

Insgesamt betrachtet liegen die Werte der zwischen -1 und 1 normierten Zusammenhangsmaße zwischen -0,27 und -0,48, was auf einen negativen, schwachen bis allenfalls mittleren Zusammenhang schließen lässt. Höhere Werte auf der Skala der Lebenszufriedenheit gehen also, wie bereits erkannt, mit geringen Werten auf der Gesundheitsskala (und damit einer besseren Gesundheit!) einher. Auch diese verwendeten Maße testen dabei ungerichtet.

Zusammenhangsmaße für metrisch-skalierte Variablen

Um beispielhaft einen möglichen Zusammenhang zwischen zwei metrisch-skalierten Variablen zu berechnen, nutzen wir nun den Datensatz Prestige, den wir bereits kennengelernt haben. Um diesen Datensatz - der Teil des Zusatzpaketes car ist - verwenden zu können, müssen wir dieses Paket zunächst wieder in den Arbeitsspeicher laden. Danach können wir diesen mit der Funktion data() aufrufen. Der Befehl head() zeigt uns die ersten Zeilen des Datensatzes an. Mit dem Befehl View() ließe sich wiederum die komplette Datenbank anzeigen.

library(car)
data(Prestige)
head(Prestige)

##                     education income women prestige census type
## gov.administrators      13.11  12351 11.16     68.8   1113 prof
## general.managers        12.26  25879  4.02     69.1   1130 prof
## accountants             12.77   9271 15.70     63.4   1171 prof
## purchasing.officers     11.42   8865  9.11     56.8   1175 prof
## chemists                14.62   8403 11.68     73.5   2111 prof
## physicists              15.64  11030  5.13     77.6   2113 prof

Im Folgenden wird uns eine mögliche Korrelation zwischen den beiden metrischen Variablen prestige und education interessieren. Die Variable prestige gibt dabei den “Pineo-Porter-Prestigewert für Berufe” an, der aus einer Mitte der 1960er Jahre durchgeführten Sozialerhebung stammt. Die Variable education gibt die “durchschnittliche Ausbildung der Berufsinhaber in Jahren” an und wurde im Jahr 1971 erhoben. Diese Informationen kann man der Hilfeseite entnehmen, die mit dem Befehl ?Prestige aufrufbar ist.

Mit der Funktion summary() können wir uns zunächst einen Überblick über einige gängige univariate Lagemaße verschaffen. Auch ließe sich der table() oder xtabs()-Befehl verwenden (nicht gezeigt), wobei hier aufgrund der vielen Merkmalsausprägungen die Übersichtlichkeit leidet.

summary(Prestige$prestige)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   14.80   35.23   43.60   46.83   59.27   87.20

Der Befehl zeigt folgende Informationen: Das Minimum der Variable prestige liegt bei 14,8%; der 1. Quartilswert liegt bei 35,23%; der Median (2. Quartilswert) liegt bei 43,6%; das arithmetische Mittel liegt bei 46,83%, der 3. Quartilswert liegt bei 59,27% und das Maximum bei 87,2%.

summary(Prestige$education)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   6.380   8.445  10.540  10.738  12.648  15.970

Auch die Darstellung innerhalb einer Kreuztabelle wäre aufgrund der vielen Merkmalsausprägungen sehr unübersichtlich. So kann es auch hier, wie bereits im Kapitel zu den ordinalen Zusammenhangsmaßen gezeigt, sinnvoll sein, die Variable in Kategorien bzw. Klassen einzuteilen. Hierzu eignet sich aufgrund der metrischen Skalierung z.B. die Einteilung in Quartile, die wir als zusammenfassende Darstellung von Lagemaßen oben mit dem summary()-Befehl zur Verfügung gestellt bekommen haben.

Neben der bereits gezeigten Verwendung der Befehle mutate und case_when aus dem Paket dplyr eignet sich zur schnellen Einteilung für die tabellarische Darstellung von Variablen der cut()-Befehl. (Siehe ?cut für weitere Informationen). Die darüber erstellte Variable lässt sich dann jedoch aufgrund der besonderen Speicherung in Klassen-Abschnitte nicht mehr für weitere Berechnungen verwenden. Die Cut-Off-Werte können manuell unter breaks angegeben werden, wobei zu beachten ist, dass auch das Minimum und das Maximum als untere bzw. obere Grenze anzugeben sind. Außerdem muss aufgrund der speziellen Ziehung der Grenzen innerhalb der Standardeinstellungen der Funktion darauf geachtet werden, einen etwas geringeren als den Minimum-Wert anzugeben:

Prestige$educ_kat <- cut(Prestige$education, 
                        breaks = c(6.3,8.445,10.54,12.648,15.97))

Die Verteilung der neuen Variable educ_kat sieht dann folgendermaßen aus:

table(Prestige$educ_kat)

## 
##  (6.3,8.45] (8.45,10.5] (10.5,12.6]   (12.6,16] 
##          26          25          25          26

Aufgrund der Einteilung in Quartile hat jede Werteklasse nun (annähernd) dieselbe Fallanzahl. “(” bedeutet dabei ausgeschlossen, “]” inklusive. Aufgrund eben dieser Einteilung mussten wir oben ein etwas kleineres Minimum wählen, da der kleinste Wert ansonsten bei der Ziehung der Grenze der ersten Werteklasse mit einer “(” ausgeschlossen worden wäre. So gibt es nun also 26 Fälle im Datensatz, bei denen die Befragten mehr als 6,3 Bildungsjahre und bis zu 8,45 Bildungsjahre aufweisen. Außerdem 25 Befragte, die mehr als 8,45 Bildungsjahre und bis zu 10,5 Bildungsjahren angegeben haben, usw.

Dieser Vorgang lässt sich mit den Werten für die Lagemaße der Variablen prestige wiederholen:

Prestige$prestige_kat <- cut(Prestige$prestige, 
                            breaks = c(14.7,35.23,43.6,59.27,87.2))

Die Verteilung der neuen Variable prestige_kat sieht dann folgendermaßen aus:

table(Prestige$prestige_kat)

## 
## (14.7,35.2] (35.2,43.6] (43.6,59.3] (59.3,87.2] 
##          26          25          25          26

Nun lassen sich die Variablen übersichtlich innerhalb einer Kreuztabelle abbilden:

table(Prestige$prestige_kat,Prestige$educ_kat, 
      dnn = c("Prestige-Kategorien", "Bildungsjahre-Kategorien"))

##                    Bildungsjahre-Kategorien
## Prestige-Kategorien (6.3,8.45] (8.45,10.5] (10.5,12.6] (12.6,16]
##         (14.7,35.2]         15           9           2         0
##         (35.2,43.6]          9          10           6         0
##         (43.6,59.3]          2           6          12         5
##         (59.3,87.2]          0           0           5        21

Zu erkennen ist nun eine deutliche Häufung von Fällen in der Diagonalen, was darauf hindeutet, dass das Berufsprestige steigt, je mehr Bildungsjahre eine befragte Person aufweist.

Für die Berechnung ungerichteter metrischer Zusammenhänge eignet sich der gängige “Korrelationskoeffizient nach Bravais-Pearson” (auch: “Pearson´s r”). Auch hierfür können Sie den bereits oben verwendeten cor-Befehl nutzen und müssen nur die Berechnungsmethode unter method ändern:

• der Korrelationskoeffizient nach Bravais-Pearson:

cor(x = Prestige$prestige, y = Prestige$education, 
    use = "pairwise.complete.obs", method = "pearson")

## [1] 0.8501769

Der zwischen -1 und 1 normierte Korrelationskoeffizient nach Bravais-Pearson zeigt hier einen Wert von 0,85, was einen starken, positiven Zusammenhang schlussfolgern lässt. Somit bestätigt sich unsere Vermutung von oben.

Möglich ist es übrigens ebenso, mehr als zwei Variablen gleichzeitig zu berechnen, wobei dann eine Korrelationsmatrix für jeden paarweisen Vergleich bereitgestellt wird:

cor(Prestige[, c("prestige","education","income","women")] , 
     method = c("pearson"))

##             prestige  education     income       women
## prestige   1.0000000 0.85017689  0.7149057 -0.11833419
## education  0.8501769 1.00000000  0.5775802  0.06185286
## income     0.7149057 0.57758023  1.0000000 -0.44105927
## women     -0.1183342 0.06185286 -0.4410593  1.00000000

• Kovarianz

Sollten Sie nur die Kovarianz benötigen, können Sie diese analog mit der Funktion cov() ermitteln. Weitere Einzelheiten finden Sie auf der Hilfeseite ?cov.

cov(x = Prestige$prestige, y = Prestige$education, 
    use = "pairwise.complete.obs", method = "pearson")

## [1] 39.90856

• lineare Regression

Auch die lineare Regressionsanalyse kann für die Berechnung von Zusammenhängen metrischer Variablen eingesetzt werden. Anders als alle anderen hier gezeigten Korrelationsmaße lassen sich Zusammenhänge damit gerichtet anstatt nur ungerichtet testen. Da sie ein eigenes, wichtiges und oft genutztes Testverfahren innerhalb der induktiven (schließenden) Statistik ist, sind ihr die letzten drei Module dieses Kurses gewidmet.

Zusammenhangsmaße bei unterschiedlichen Skalenniveaus

Nicht in allen Fällen werden Variablen, deren Korrelation zueinander berechnet werden soll, auch dasselbe Skalenniveau aufweisen. Hier ist es zum einen möglich, die Variable mit dem höheren Niveau zu dem niedrigeren Niveau abzustufen. Bspw. lässt sich für eine metrische Variable durch die Einteilung in Klassen eine Rangbildung vornehmen, sodass die Variable nun ein ordinales Skalenniveau aufweist und ein potenzieller Zusammenhang mit einer weiteren ordinal-skalierten Variable berechnet werden kann. In diesem Fall werden Maße verwendet, die (mindestens) für das Ordinal-Skalenniveau geeignet sind. Zum anderen gibt es neuere Verfahren, wie bspw. die “Biseriale Rangkorrelation” (Variablen: dichotom/ordinal), die “Punktbiseriale Korrelation” (Variablen: dichotom/metrisch) und “Eta” (Variablen: nominal/metrisch), die an dieser Stelle nur kurz erwähnt werden sollen, da sie die Ausmaße dieses Kurses sprengen würden.

---

Dieses Modul ist Teil des Kurses: Dein SoWi-Statistik-Kurs: Von der Theorie zur Praxis in R - Materialien zum Selbstlernen und Unterrichten. Carl von Ossietzky Universität Oldenburg. Für eine eventuelle Vervielfältigung, Veränderung oder Weitergabe der Kursinhalte beachten Sie bitte die Lizenzbestimmungen nach CC BY 4.0 (siehe dazu auch die Modulbeschreibung dieses Kurses).