Diese Webseite ist archiviert. Unter Umständen sind nicht alle Inhalte barrierefrei. Weitere Informationen finden Sie in unserer Erklärung zur Barrierefreiheit.

Vorbemerkungen

Dieses kurze Dokument ersetzt keine Einführung in das R Grafiksystem, im Folgenden geht es lediglich um ein paar kurze Anregungen und Beispiele, die für die Erstellung eigener Grafiken nützlich sein könnten.

Es handelt sich auch um keine systematische Übersicht, sonderen es werden lediglich einige ad hoc Beispiele gezeigt.
Für eine erste Einführung und einfache Beispiele siehe z.B. Quick-R https://www.statmethods.net/graphs/index.html. Für allgemeine Hinweise zur Visalisierung von Daten siehe z.B. Schwabish (2014).

R Grafiken (allgemeine Grundlagen)

Einer der vielen Vorzüge von R besteht darin, dass man damit (relativ) einfach hochwertige Grafiken erstellen kann, die auch professionellen Ansprüchen genügen.

Die Art, wie die Grafiken erstellt werden, unterscheidet sich jedoch deutlich von Programmen wie z.B. Excel, denn auch Grafiken werden in R ‘programmiert’.

Mittlerweile existieren mehrere sehr mächtige R-packages zur Erstellung von Grafiken (z.B. ggplot2), denen zum Teil eine etwas andere Philosophie zugrunde liegt, doch hier werden wir uns auf die Fähigkeiten von base R beschränken.1

In base R erfolgt der Aufbau von Grafiken ähnlich wie mit Papier und Bleistift, man zeichnet z.B. Linien und Punkte auf ein Papier, und sobald sie gezeichnet sind sind diese fixiert, aber man kann andere Elemente dazu- und darüberzeichnen.

Prinzipiell wird in base R zwischen high level und low level Plot-Funktionen unterschieden.

High level Grafik Funktionen

High level Funktionen legen die Grundstruktur der Grafiken fest, z.B. die Länge der Achsen, und erzeugen in der Regel vollständige Grafiken. Die wichtigste Funktion hier ist plot(), die je nach Typ des Datenobjekts eine Grafik an das aktive Grafik-Fenster (device) schickt.

Zur Illustration erzeugen wir anfangs zwei gleichverteilte Zufallsvariablen x und y und ploten diese aud vier verrschiedene Arten

set.seed(12345)
x <- floor(runif(n = 20, min = 1, max = 10))
y <- 0.5*x + runif(n = 20, min = 1, max = 10)

par(mfrow = c(2,2))
plot(x, main = "plot(x)")
plot(x, y, main = "plot(x, y)")
plot(factor(x), main = "plot(factor(x))")
plot(factor(x), y, main = "plot(factor(x), y)")

Was ist hier passiert? Mit set.seed(12345) wir der Zufallsgenerator initialisiert, so dass bei einem neuerlichen Aufruf wieder die gleichen Zufallszahlen generiert werden und somit die Ergebnisse reproduziert werden können. Mit runif() werden gleichverteilte Zufallszahlen erzeugt, und floor() gibt die nächst kleinere ganze Zahl zurück, wir erhalten also einen Integer.

Mit par() können Grafikparameter wie z.B. die Ränderbreite eingestellt und übergeben werden, hier wird mit der Option par(mfrow = c(2,2)) eingestellt, dass die nächsten vier Grafiken matrixartig in einer Grafik mit zwei Zeilen (hoizontal) und zwei Spalten (vertikal) zusammengefasst werden.

x ist eine numerische Variable, und mit plot(x) wird diese einfach gegen den Index (d.h. die Beobachtungen 1:20) aufgetragen. Mit plot(x, y) wird ein Streudiagramm (Scatterplot) erzeugt, d.h. die Punkte von x werden gegen y aufgetragen.

Mit factor(x) wird aus x eine kategoriale Variable erzeugt, und wird plot() mit einer (nominalskalierten) kategorialen Variable aufgerufen erzeugt R automatisch ein Balkendiagramm mit den absoluten Häufigkeiten.

Mit plot(factor(x), y) wird für jede Ausprägung der kategorialen Variable factor(x) ein ensprechender Boxplot für y erzeugt.

Neben der einfachen plot() Funktion gibt es zahlreiche weitere high level Plot Funktionen, z.B. barplot() (Säulendiagramm), dotchart() (Punkteplot), hist() (Histogramm), boxplot(), pie(), …

Für alle diese Funktionen gibt es Optionen und Einstellungen, mit denen die Grafiken modifiziert werden können. Keine Sorge, Sie müssen all die folgenden Optionen und Grafikparameter nicht auswendig lernen, sie werden Ihnen z.B. von RStudio nach drücken der <Tab> Taste vorgeschlagen, Sie sollten aber eine Ahnung haben, was die Optionen bewirken. Im Folgenden bezeichnet string jeweils eine beliebige Zeichenkette.

Option Beschreibung
main = "string" Titel (Überschrift) der Grafik
sub = "string" Untertitel
xlab = "string" x-Achsenbeschriftung (label)
ylab = "string" y-Achsenbeschriftung
xlim = c(x1, x2) Wertebereich der x-Achse (von - bis)
ylim = c(y1, y2) Wertebereich der y-Achse; c(y2, y1) um die Achse “umzudrehen”
axes = TRUE/FALSE Achsen anzeigen oder unterdrücken
ann = TRUE/FALSE unterdrückt Beschriftungen (Titel, labels), für “annotation”
frame.plot = TRUE/FALSE Rahmen anzeigen oder nicht
type = "." l = Linien, p = Punkte, b = Linien + Punkte, n = keine (nur Achsen), o = Überlagerung, h = vertikale Verbindungen zur x-Achse (histogrammartig), s = Treppenfunktion

Die wichtigsten Grafiktypen sind sind der folgenden Abbildung dargestellt

par(mfrow = c(3,3))
plot(x, type = "l", main = "type = \"l\"")
plot(x, type = "p", main = "type = \"p\"")
plot(x, type = "b", main = "type = \"b\"")
plot(x, type = "o", main = "type = \"o\"")
plot(x, type = "h", main = "type = \"h\"")
plot(x, type = "c", main = "type = \"c\"")
plot(x, type = "s", main = "type = \"s\"")
plot(x, type = "S", main = "type = \"S\"")
plot(y, type = "n", main = "type = \"n\"")

Low level Plot Funktionen

Mit Low level Plot Funktionen können Elemente zu einer high level Funktion hinzugefügt werden, oder es können damit auch neue Grafiken von Grund auf neu erzeugt werden.

Einige der wichtigsten low level Plot Funktionen sind

Funktion Beschreibung
points(x, y) zeichnet Punkte, äquivalent zu plot(x,y, type = "p")
lines(x, y) zeichnet Linien, äquivalent zu plot(x,y, type = "l")
abline() zeichnet Linien wie z.B. Regressionsgeraden
curve() zeichnet eine Funktion
arrows() Pfeile
text(x, y, labels) Text in der Grafik
mtext() Text in den Rändern zwischen den Bereichen
grid() Gitternetz
legend() Legende
title(main, sub) ezeugt Titel und Untertitel
locator() Ermittelt durch Klick in eine Grafik die Koordinaten
identify() Bei Klicken in die Nähe des definierten Punktes wird eine Ausgabe erzeugt

Die jeweiligen Optionen etc. erhalten Sie mit der Hilfefunktion von R, z.B. ?abline.

Wenn wir z.B. die Punkte beschriften möchten können wir dies mit der text() Funktion machen. Mit LETTERS[1:20] erhalten wir die ersten 20 Großbuchstaben, also

plot(x, y)
text(x, y, labels = LETTERS[1:20], pos = 4)
abline(h = 4, lty = 2, lwd = 2, col = "blue")
abline(v = 6, lty = 3, lwd = 3, col = "darkblue")
abline(lm(y ~ x), col = "red")

Die vier Seiten von Grafiken werden in R meist von unten im Uhrzeigersinn durchnummeriert, die Option pos = 4 der text() Funktion gibt an, dass die Beschiftung rechts erfolgen soll (1 = unten, 2 = links, 3 = oben, 4 = rechts).

Die beiden Funktionen abline() zeichnen ein horizontale (bzw. vertikale) Line vom Linientypus 2 (= strichliert; bzw. 3 = Punkte) in der Linenstärke (linewidth) 2 und in der Farbe (dunkel-) blau. Die dritte abline() Funktion zeichnet schließlich eine rote Regressionsgerade (lm()-Befehl) ein.

Einige wichtige Grafikparameter

Darüber sind viele weitere Einstellungen verfügbar, z.B. (siehe auch https://www.statmethods.net/advgraphs/parameters.html)

Argument Beschreibung
pch = n Plotsymbol für Punkte 1 bis 25 (plotcharacter); eine Übersicht finden Sie unten
lty = n (für linetype) eine Zahl 0-6 oder die engl. Bezeichnung in Anführungszeichen; 1 = solid, 2 = dashed (gestrichelt), 3 = dotted (gepunktet), 4 = dotdash, 5 = longdash, 6 = twodash, 0 = blank (unsichtbar)
lwd = n (für linewidth) Linienstärke, eine Zahl größer Null für das Vielfache der Standardstrichstärke
adj = n Textausrichtung: 0 = links, 0.5 = zentriert, 1 = rechts
las = n Ausrichtung der Tickmarks (Label)-Beschriftung: 0 = parallel zur Achse, 1 = horizontal, 2 = senkrecht zur Achse, 3 = vertikal
tck = n Länge der Tickmarks, je nach Vorzeichen weisen diese nach innen oder nach außen
font = n Code: 1 = normal, 2=fett, 3=kursiv, 4 = fett + kursiv, 5 = Symbol
font.axis = n Font-Code Achse
font.main = n Font-Code Titel
font.sub = n Font-Code Subtitel
lfont.lab = n Font-Code für die x- und y-Achsen-Labels
family = ... Name der Fontfamilie (z.B. serif, symbol) für Text
col = ... Wahl der Farbe; die Auswahl erfolgt entweder per englischer Bezeichnung (red, lightred, blue, darkblue, yellow, …) oder wie in html: RRGGBB mit R/G/B-Anteil oder mit der Angabe einer Zahl (1:n)
col.axis = ... Farbe für die Achsenerklärung
col.lab = ... Farbe für die x- und y-Achsen-Labels
col.main = ... Farbe für den Titel
col.sub = ... Farbe des Subtitels
cex = n n = Vergrößerungs- (>1) bzw. Verkleinerungsfaktor (<1) für Text, d.h. mit n < 1 wird der Text entsprechend verkleinert, mit n > 1 vergrößert
cex.main = n Größe Titel
cex.sub = n Größe Font-Code Achse
cex.lab = n Größe Achsenbeschriftungen
cex.axis = n Größe Tickmark-Labels

Eine Übersicht über plotcharachters (pch):

Die obigen Parameter können bei fast allen Grafikfunktionen direkt als Option eingestellt werden, oder sie können vorab mit einem eigenen Befehl par() eingestellt werden.

Werden die Optionen mit par() eingestellt bleiben diese wirksam bis der Grafik-Device wieder geschlossen wird.

Einige Grafik-Optionen können nur mit par() eingestellt werden, z.B. die Ränderbreite von Grafiken oder die Anordnung bei Mehrfachgrafiken (par(mfrow = c(...))

Die Ränder von einfachen Grafiken können z.B. mit par(mar = c(2,4,1,2) verändert werden, wobei die Zahlen bedeuten, dass der untere Rand 2 Zeilenstärken beträgt, der linke Rand 4 Zeilenstärken, und der obere Abstand eine und der rechte Rand 2 Zeilenstärken betragen soll.

Der Befehl par(mai = c(0.1, 0.2, 0.1, 0.2)) funktioniert im wesentlichen gleich, nur werden die Abstände hier in inch angegeben, also unten und oben jeweils 0.1 Inch und links und rechts jeweils 0.2 Inch.

Mit par(oma = c(unten, links, oben, rechts)) können die äußeren Ränder von Mehrfachgrafiken eingestellt werden.

Der folgende Code, welcher die Abbildung unten erzeugt, wurde Braun and Murdoch (2008) entnommen

par(mar = c(5, 5, 5, 5) + 0.1)
plot(c(1, 9), c(0, 50), type = "n", xlab = "", ylab = "")
text(6, 40, "Plot region")
points(6, 20)
text(6, 20, "(6, 20)", adj = c(0.5, 2))
mtext(paste("Margin", 1:4), side = 1:4, line = 3)
mtext(paste("Line", 0:4), side = 1, line = 0:4, at = 3, cex = 0.6)
mtext(paste("Line", 0:4), side = 2, line = 0:4, at = 15, cex = 0.6)
mtext(paste("Line", 0:4), side = 3, line = 0:4, at = 3, cex = 0.6)
mtext(paste("Line", 0:4), side = 4, line = 0:4, at = 15, cex = 0.6)

Für weitere Einstellungen siehe z.B. in der R Hilfe ?par oder https://www.rdocumentation.org/packages/graphics/versions/3.6.2/topics/par.

Wenn man nachsehen möchte, welche Einstellung gerade wirksam ist, kann man par() mit der jeweiligen Option in Anführungszeichen aufrufen, z.B.

par("mar")
## [1] 5.1 4.1 4.1 2.1

Mit par(no.readonly = TRUE) erhält man alle verstellbaren Grafikeinstellungen.

Die mit der par() Funktion getroffenen Einstellungen bleiben wirksam, bis der Grafik device wieder geschlossen wird.

Möchte die Einstellungen vorher wieder zurücksetzen kann man die verstellbaren Einstellungen am Beginn der Sitzung speichern und später wieder zurücksetzen

oldpar <- par(no.readonly = TRUE) # speichern
[...]                             # diverse Befehle
par(oldpar)                       # Zurücksetzen

Beispiel: Fragebogen

Zur Einstimmung zeigen wir, wie man einige typische Grafiken für die Auswertung eines Fragebogens erstellt.

Dazu laden wie die Datendatei im (deutschen) csv-Format einmalig aus dem Internet und speichern sie im default-Verzeichis.

## Auswertung eines Fragebogens
rm(list = ls())
# setwd(/mydirectory)

if (!file.exists("fragebogen.csv")) {
  mydf <- read.csv2("https://www.hsto.info/tmp/fragebogen.csv")
  write.csv2(mydf, file = "fragebogen.csv", row.names = FALSE)
} else {
  mydf <- read.csv2("fragebogen.csv")
}
str(mydf)  
## 'data.frame':    75 obs. of  17 variables:
##  $ Geschlecht          : chr  "männlich" "weiblich" NA "weiblich" ...
##  $ Alter               : chr  "21 Jahre oder älter" "18 Jahre oder jünger" NA "19 oder 20 Jahre" ...
##  $ Kgr                 : num  181 168 NA 1.65 165 173 162 176 170 179 ...
##  $ Kgr_vater           : num  177 173 NA 1.72 195 172 186 182 180 180 ...
##  $ Eink_erw            : int  1600 2500 NA 2800 1600 1700 1700 1900 1800 1800 ...
##  $ Herkunft            : chr  "Südtirol" "Süddeutschland (Bayern, Baden-Württemberg) oder Schweiz" NA "Aus einer anderen Region" ...
##  $ Info_uibk           : chr  "Ja" "Ja" NA "Nicht Gewählt" ...
##  $ Info_umit           : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_interet        : chr  "Ja" "Nicht Gewählt" NA "Ja" ...
##  $ Info_medien         : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_schule         : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_freunde        : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_messen         : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_anderes        : chr  NA NA NA NA ...
##  $ Interesse_wirtschaft: int  2 2 NA 3 3 2 2 2 4 5 ...
##  $ Interesse_gesundheit: int  3 5 NA 2 4 1 1 2 3 3 ...
##  $ Interesse_sport     : int  1 5 NA 1 1 1 1 1 5 2 ...
str(mydf)
## 'data.frame':    75 obs. of  17 variables:
##  $ Geschlecht          : chr  "männlich" "weiblich" NA "weiblich" ...
##  $ Alter               : chr  "21 Jahre oder älter" "18 Jahre oder jünger" NA "19 oder 20 Jahre" ...
##  $ Kgr                 : num  181 168 NA 1.65 165 173 162 176 170 179 ...
##  $ Kgr_vater           : num  177 173 NA 1.72 195 172 186 182 180 180 ...
##  $ Eink_erw            : int  1600 2500 NA 2800 1600 1700 1700 1900 1800 1800 ...
##  $ Herkunft            : chr  "Südtirol" "Süddeutschland (Bayern, Baden-Württemberg) oder Schweiz" NA "Aus einer anderen Region" ...
##  $ Info_uibk           : chr  "Ja" "Ja" NA "Nicht Gewählt" ...
##  $ Info_umit           : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_interet        : chr  "Ja" "Nicht Gewählt" NA "Ja" ...
##  $ Info_medien         : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_schule         : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_freunde        : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_messen         : chr  "Nicht Gewählt" "Nicht Gewählt" NA "Nicht Gewählt" ...
##  $ Info_anderes        : chr  NA NA NA NA ...
##  $ Interesse_wirtschaft: int  2 2 NA 3 3 2 2 2 4 5 ...
##  $ Interesse_gesundheit: int  3 5 NA 2 4 1 1 2 3 3 ...
##  $ Interesse_sport     : int  1 5 NA 1 1 1 1 1 5 2 ...

Offensichtlich haben alle Variablen den gewünschten Typ und den erwarteten Inhalt.

Bevor wir mit der Auswertung beginnen speichern wir die Vorgabewerte der Grafikparameter

oldpar <- par(no.readonly = TRUE) # ursprüngliche Grafikparameter speichern

Nun können wir mit der Auswertung beginnen. Normalerweise würde man mit der zentralen Botschaft der Untersuchung beginnen, um die Aufmerksamkeit der Leserin zu gewinnen, aber hier werden wir aus eher didaktischen Gründen mit kategorialen (nominal skalierten) Fragen beginnen.

Kategoriale (nominal skalierte) Variablen

Für kategoriale Variablen wird in R der Datentyp factorverwendet. Die Unterschiedlichen Ausprägungen können mit levels() abgefragt werden.

Im folgenden einfachen Beispiel wird eine Variable vom Datentyp character angelegt, dann wird daraus ein factor-Objekt erzeugt, und mit labels werden die ursprünglichen levels überschrieben. Gespeichert werden nur die (neuen) levels, diese können wir mit dem levels() Befehl abfragen.

x <- c("w", "w", "m", "w", "w", "m")
str(x)
##  chr [1:6] "w" "w" "m" "w" "w" "m"
x <- factor(x, levels = c("w", "m"), labels = c("female", "male"))
str(x)
##  Factor w/ 2 levels "female","male": 1 1 2 1 1 2
levels(x) <- c("weiblich", "männlich")
str(x)
##  Factor w/ 2 levels "weiblich","männlich": 1 1 2 1 1 2

Man beachte, dass die Reihenfolge der levels und labels übereinstimmen muss!

Geschlecht

Mit str(mydf$Geschlecht) oder class(mydf$Geschlecht) können wir überprüfen, ob es sich tatsächlich um ein factor Objekt (also eine kategoriale, bzw. nominal skalierte) Variable mit zwei Ausprägungen handelt. Auf die Ausprägungen von factor Objekten kann mit levels() zugegriffen werden,

levels(mydf$Geschlecht)
## NULL

und mit dem table() Befehl können wir die Auspräungen auszählen

table(mydf$Geschlecht)
## 
## männlich weiblich 
##       36       36

Ein einfaches Tortendiagramm könnten wir mit pie(table(mydf$Geschlecht)) erzeugen. Wenn wir es auch noch beschriften möchten können wir uns die prozentuellen Anteile berechnen (pct) und diese mit dem paste() Befehl an die levels anfügen und als lbls speichern. Außerdem müssen wir mit par(mar = c(unten, links, oben, rechts)) die Ränder anpassen

pct <- round(prop.table(table(mydf$Geschlecht))*100)     # proz. Anteil
lbls <- paste(levels(mydf$Geschlecht), pct,"%", sep=" ") # Label
# x11()
# pdf(file = "G_Geschlecht.pdf")
pie(table(mydf$Geschlecht), labels = lbls, cex=0.75, cex.main = 0.8,
    main = "Geschlecht")

# dev.off()

Mit x11() könnte auf Windows Rechnern ein eigenes Fenster als Grafik-device geöffnet werden, für macOS muss Quartz installiert werden, dann kann mit quartz() die Ausgabe in ein Grafikfenster umgeleitet werden (siehe ?quartz). Mit pdf(file = "...") kann die Grafik automatisch im pdf-Format im default-Verzeichnis gespeichert werden, mit dev.off() muss diese Ausgabeeinheit wieder geschlossen werden.

Ähnlich können die Grafiken mit jpeg(), bmp(), png() und tiff() im entsprechenden Format gespeichert werden.

Alter

Ebenso können wir die Altersverteilung (Variable mydf$Alter) darstellen. Zur Gegenüberstellung lassen wir uns hier ein Torten- und ein Balkendiagramm ausgeben, und einer Publikation würden Sie sich natürlich für eine der beiden Darstellungen entscheiden.

Vorher überschreiben wir die bestehenden levels, da diese zu lang sind. Mit dem Steuerzeichen "\n" können wir in eine Zeichenkette einen Zeilenumbruch einfügen.

# Alter
# ==================================================================
levels(mydf$Alter)
## NULL
levels(mydf$Alter)[1] <- "18 Jahre \noder \njünger"
levels(mydf$Alter)[2] <- "19 Jahre \noder 20\nJahre"
levels(mydf$Alter)[3] <- "21 Jahre \noder \nälter"

pct <- round(prop.table(table(mydf$Alter))*100)   # proz. Anteile
lbls <- paste(levels(mydf$Alter), " (", pct,"%)",sep="")  # Labels

par(mfrow = c(1, 2), cex = 0.8, cex.axis = 0.7)
# pdf(file = "G_Alter.pdf")
pie(table(mydf$Alter), labels = lbls, main = "Alter")
barplot(table(mydf$Alter), main = "Alter")

Generell wird von der Verwendung von Tortendiagrammen (pie) abgeraten, da sich das menschliche Auge/Gehirn diese visuelle Information nicht so gut verarbeiten kann, Balkendiagramme eignen sich in der Regel besser.

Wo befand sich ihr Lebensmittelpunkt vor Aufnahme des Studiums?

# L01: nominal skaliert
levels(mydf$Herkunft)   # bestehende levels
## NULL
# ==================================================================
# levels durch kürzere mit Zeilenumbrüchen ersetzen
levels(mydf$Herkunft)[1] <- "Landeck und \nUmgebung (30km)"
levels(mydf$Herkunft)[2] <- "Innsbruck \nund Umgebung"
levels(mydf$Herkunft)[3] <- "Rest Nord- & \nOsttirol"
levels(mydf$Herkunft)[4] <- "Rest Österreich"
levels(mydf$Herkunft)[5] <- "Südtirol"
levels(mydf$Herkunft)[6] <- "Süddeutschland \nSchweiz"
levels(mydf$Herkunft)[7] <- "Rest der Welt"
levels(mydf$Herkunft)
## [1] "Landeck und \nUmgebung (30km)" "Innsbruck \nund Umgebung"     
## [3] "Rest Nord- & \nOsttirol"       "Rest Österreich"              
## [5] "Südtirol"                      "Süddeutschland \nSchweiz"     
## [7] "Rest der Welt"
# x11()
# pdf(file = "G_woher.pdf")
par(mai = c(0.5, 2.2, 0.8, 0.5)) # Ränder setzen
xlim_neu <- c(0, 1.2*max(table(mydf$Herkunft))) # x-Achse verlängern
pct <- paste(round(prop.table(table(mydf$Herkunft)),2)*100, "%") # proz. Anteile

BP_Herkunft <- barplot(rev(table(mydf$Herkunft)), horiz = TRUE, las = 1, 
                  xlim = xlim_neu, cex.lab = 0.8, cex.main = 1, 
                  main = "Wo befand sich ihr Lebensmittelpunkt \nvor Aufnahme des Studiums?")

text(x = rev(table(mydf$Herkunft)), y = BP_Herkunft, label = rev(pct), 
     pos = 4, cex = 0.8, col = "black")

# dev.off()
mit rev() wird die Reihenfolge umgekehrt (also mit Landeck begonnen), horiz = TRUE werden horizontale Balken gezeichnet, las = 1 sorgt für horizontale Labels.

 

Der Zusammenhang zwischen nominal skalierten Variaben kann grafisch mit einem spineplot() oder mosaicplot() dargestellt werden

par(las = 1, cex = 0.8, mar = c(2, 8, 4, 2)+0.1)
spineplot(table(mydf$Geschlecht, mydf$Herkunft), xlab = "", ylab = "", 
          yaxlabels = c("Landeck", "Innsbruck", "Rest Tirol", "Rest Österr.", "Südtirol", 
                        "Süddtld, Schweiz", "Andere"))

Um zu erkennen, ob es sich bei den Unterschieden nur um Stichprobenrauschen handelt oder ob sich dahinter möglicherweise eine Systematik verbirgt, könnte z.B. ein \(\chi^2\)-Test (chi-squared Test) durchgeführt werden. Der \(\chi^2\)-Test ist ein asymptotischer Test, also nur in großen Stichproben näherungsweise gültig. Er erfordert, dass die Einträge der Zellen zumindest annähernd normalverteilt sind.

chisq.test(table(mydf$Geschlecht, mydf$Herkunft))
## Warning in chisq.test(table(mydf$Geschlecht, mydf$Herkunft)): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(mydf$Geschlecht, mydf$Herkunft)
## X-squared = 3.0108, df = 6, p-value = 0.8075

R warnt uns, dass die Approximation möglicherweise ungenau ist. Der Grund dafür ist, dass einzelne Zellen zu wenige Beobachtungen enthalten. Aber wenn wir dem Ergebnis vertrauen würden, würden wir es folgendermaßen interpretieren:

Interpretation: Wenn die Nullhypothese wahr ist (also es keine systematischen Unterschiede in der Herkunft zwischen Männern und Frauen gibt), und wenn alle dem Test zugrunde liegenden Annahmen erfüllt sind, würden wir bei hypothetisch wiederholten Stichprobenziehungen in ca. 80% der Fälle einen ähnlich extremen oder noch extremeren Wert der Teststatistik erwarten. Nachdem dies kein sehr unwahrscheinliches Ergebnis ist würden wir die Nullhypothese nicht verwerfen und schließen, dass es keine systematischen Unterschiede in der Herkunft zwischen Männern und Frauen gibt.

Hinweis: bei wenigen Beobachtungen kann Fisher’s exakter Test der Nullhypothese, dass Zeilen und Spalten einer Kontingenztabelle unabbhängig sind, herangezogen werden

fisher.test(table(mydf$Geschlecht, mydf$Herkunft))
## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(mydf$Geschlecht, mydf$Herkunft)
## p-value = 0.8239
## alternative hypothesis: two.sided

Die Interpretation ist gleich wie bei dem \(\chi^2\)-Test, da der p-Wert (0.8239) viel größer als 0.05 ist kann die Nullhypothese, dass es in der Grundgesamtheit keine Unterschiede der Herkunft nach Geschlechtern gibt, nicht verworfen werden.

Als deskriptive Kennzahl könnten wir auch ein Kontingenzmaß berechnen, z.B. Cramers V. Dies geht am einfachsten mit dem Package vcd (für visualizing categorical data) und den darin definierten Befehl assocstats(), der auch gleich zwei Teststatistiken berechnet

library(vcd)
assocstats(table(mydf$Geschlecht, mydf$Herkunft))
##                     X^2 df P(> X^2)
## Likelihood Ratio 3.0692  6  0.80012
## Pearson          3.0108  6  0.80750
## 
## Phi-Coefficient   : NA 
## Contingency Coeff.: 0.202 
## Cramer's V        : 0.206

Fragen mit Mehrfachantworten

Wenn bei einer Frage Mehrfachantworten zugelassen wurden, so wird jede Antwortmöglichkeit in einer eigenen Variable (Spalte) gespeichert mit den zwei levels "Ja" und "Nicht gewählt". Deshalb müssen wir zuerst mit Hilfe des table() Befehls jede Antwortmöglichkeit einzeln auszählen und die Ergebnisse in einem Vektor ablegen, wobei wir gleich den einzelnen Elementen Namen zuordnen.

# Frage mit Mehrfachantworten: zuerst auszählen
wie_erf <- c("Homepage der Universität \nInnsbruck (UIBK)" =
               table(mydf$Info_uibk)[[1]],
          "Homepage der UMIT" = table(mydf$Info_umit)[[1]],
          "Andere Internetquellen" = table(mydf$Info_interet)[[1]],
          "Massenmedien (Zeitung, \nRadio, Fernsehen, etc)" =
               table(mydf$Info_medien)[[1]],
          "In der Schule" = table(mydf$Info_schule)[[1]],
          "Von Freunden, \nVerwandten, etc" =
               table(mydf$Info_freunde)[[1]],
          "Informationsveranstaltungen, \nMessen, etc" =
               table(mydf$Info_messen)[[1]],
          "Sonstiges" = length(table(mydf$Info_anderes))
)

Der wesentliche Unterschied zwischen doppelten eckigen Klammern [[ ]] und einzelnen eckigen Klammern [ ] (bzw. dem $ Zeichen) besteht darin, dass erstere immer nur ein Element zurückgeben (in diesem Fall, wie oft "Ja" angekreuzt wurde), während letztere mehrere Elemente zurückgeben können.

Unter "Sonstiges" waren freie Eingaben möglich, deshalb zählen wir hier mit length() nur die Anzahl der Eintragungen.

Das Ergebnis können wir nun z.B. in einem Balkendiagramm ausgeben

par(mai = c(1.1, 2.5, 0.8, 0.5)) # Ränder einstellen, von unten im Uhrzeigersinn
xlim_neu <- c(0, 1.1*max(wie_erf))  # Länge der x-Achse

BP_wie_erf <- barplot(rev(wie_erf), horiz = TRUE, xlim = xlim_neu, las = 1, cex.names = 0.9, 
                   main = "Wie haben Sie von diesem Bachelor-Studium erfahren?",
                   sub = "Mehrfachnennungen möglich")
text(x = rev(wie_erf), y = BP_wie_erf, label = rev(wie_erf), pos = 4, cex = 0.9, col = "black")

Ordinal skalierte (kategoriale) Fragen

Welcher Bereich hat Sie bei der Studienwahl besonderns angesprochen?

Hier wurden drei Fragen gestellt, die Ausgabe soll nebeneinander in einer Grafik erfolgen.
Um einen sinnvollen Vergleich zu ermöglichen ist darauf zu achten, dass die y-Achse in allen Grafiken gleich lang ist (ylim).

# pdf(file = "G_Interesse.pdf")
par(mfrow=c(1,3), mar = c(3.5, 1.8, 5, 1.8), oma = c(4,1,4,0))

BP_Interesse_a <- barplot(table(mydf$Interesse_wirtschaft), 
                    ylim = c(0,55), main = "Wirtschaft")
text(x = BP_Interesse_a, y = table(mydf$Interesse_wirtschaft), 
     label = table(mydf$Interesse_wirtschaft), 
     pos = 3, cex = 1, col = "black")

BP_Interesse_b <- barplot(table(mydf$Interesse_gesundheit), 
                    ylim = c(0,55), main = "Gesundheitstourismus")
text(x = BP_Interesse_b, y = table(mydf$Interesse_gesundheit), 
     label = table(mydf$Interesse_gesundheit), 
     pos = 3, cex = 1, col = "black")

BP_Interesse_c <- barplot(table(mydf$Interesse_sport), ylim = c(0,55), 
                    main = "Sporttourismus")
text(x = BP_Interesse_c, y = table(mydf$Interesse_sport), 
     label = table(mydf$Interesse_sport), 
     pos = 3, cex = 1, col = "black")

title(main = "Welcher Bereich hat Sie bei der Studienwahl besonderns angesprochen?", 
      line=-0.5, cex.main = 1.5, outer = TRUE, 
      sub = "1 = sehr stark; 5 = überhaupt nicht")

# dev.off()

Einr Rangkorrelation zwischen dem Interesse an Gesundheit und Sport können wir mit

cor(mydf$Interesse_gesundheit, mydf$Interesse_sport, method = "spearman", use = "complete.obs")
## [1] 0.6834573

Offensichtlich gibt es einen relativ starken positiven Zusammenhang zwischen diesen Interessen. Da der Test des Spearmanschen Korrelationskoeffizient bei Bindungen (ties) nicht gültig ist verwenden wir für den Test der Korrelation zwischen dem Interesse an Gesundheit und Sport Kendalls Tau

cor.test(mydf$Interesse_gesundheit, mydf$Interesse_sport, method = "kendall")
## 
##  Kendall's rank correlation tau
## 
## data:  mydf$Interesse_gesundheit and mydf$Interesse_sport
## z = 5.7428, p-value = 9.315e-09
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.5804467

Interpretation: Wenn die Nullhypothese wahr ist (d.h. es in der Grundgesamtheit keinen Zusammenhang zwischen diesen beiden Interessen gibt) und wenn alle erforderlichen Annahmen erfüllt sind, würden wir bei wiederholter Durchführung des zugrunde liegenden Zufallsexperiments (d.h. Stichprobenziehungen) in 0.0000..:% der Fälle einen ähnlich extremen oder noch extremeren Wert der Teststatistik erwarten. Da dies ein extrem unwahscheinliches Ereignis ist werden wir die Nullhypothese verwerfen und von einem systematischen Zusammenhang ausgehen.

In einer Publikation würden wir allerdings darauf vertrauen, dass die Leserin weiß wie ein Hypothesentest zu interpretieren ist, und einfach schreiben, dass die Nullhypothese eines gleich großen Interesses an Gesundheit und Sport auf dem 1% Niveau verworfen werden kann, oder noch einfacher, aber deutlich schlampiger, dass sich die Interessen statistisch hoch signifikant unterscheiden.

Metrisch skalierte Fragen

Deskriptive Übersichtsstatistiken für alle metrischen Varablen erhält man z.B. mit dem stargazer package

library(stargazer)  # erstmalig installieren!
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
stargazer(mydf, type = "html", median = TRUE, omit.summary.stat = c("p25", "p75"))
Statistic N Mean St. Dev. Min Median Max
Kgr 73 167.644 35.560 1.650 175.000 197.000
Kgr_vater 73 175.419 30.205 1.720 180.000 200.000
Eink_erw 72 2,322.222 3,930.576 600.000 1,800.000 35,000.000
Interesse_wirtschaft 73 2.753 0.983 1.000 3.000 5.000
Interesse_gesundheit 72 2.569 1.372 1.000 2.000 5.000
Interesse_sport 73 2.411 1.786 1.000 1.000 5.000

Wir interessieren uns im Moment für die Körpergröße in cm (Kgr).

Nachdem offensichtich niemand 1.65cm (!) groß hat offensichtlich jemand (oder mehrere) ihre Körpergröße in Metern angegeben. Um dies zu ändern verwenden wir den ifelse() Befehl um alle Eingaben unter 2.5 mit 100 zu multiplizieren, die restlichen Eingaben werden mit sich selbst überschrieben (bleiben also erhalten)

mydf$Kgr = ifelse(mydf$Kgr < 2.5, mydf$Kgr*100, mydf$Kgr)

Nun können wir die Körpergrößen z.B. in einem Histogramm oder einem Boxplot darstellen.

par(mfrow = c(2,1))
par(mar = c(0.5, 3, 2, 3))
hist(mydf$Kgr, freq = FALSE, main = "Körpergröße in cm", 
     col = "lightgrey", xlab = "")
par(mar = c(2, 3, 0.5, 3))
boxplot(mydf$Kgr, horizontal = TRUE, frame = FALSE, axes = FALSE, col = grey(0.9))

Gehaltsvorstellungen von Frauen und Männern

Wie wir der obigen Tabelle mit den deskriptiven Statistiken entnehmen können hat irgend ein Scherzkeks eine Einkommenserwartun von 35000 Euro. bevor wir mit der Auswertung beginnen ersetzen wir alle Werte größer 20000 mit NA

mydf$Eink_erw[mydf$Eink_erw > 20000] <- NA

Um die Einkommenserwartungen nach Geschlecht zu vergleichen zeigen wir zwei Boxplots

boxplot(mydf$Eink_erw ~ mydf$Geschlecht, main = "Gehaltsvorstellungen nach Geschlecht", cex.main = 0.8, xlab = "", ylab = "")

Offensichtlich haben Männer höhere Einkommenserwartungen als Frauen. Aber ist dies nur dem Stichprobenfehler geschuldet, oder können wir dahinter einen systematischen Zusammenhang in der Grundgesamtheit vermuten?

Ein einfacher t.test() gibt folgendes Ergebnis

t.test(mydf$Eink_erw ~ mydf$Geschlecht)
## 
##  Welch Two Sample t-test
## 
## data:  mydf$Eink_erw by mydf$Geschlecht
## t = 1.7952, df = 66.189, p-value = 0.07719
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -21.23194 399.99011
## sample estimates:
## mean in group männlich mean in group weiblich 
##               1958.824               1769.444

Wenn die Nullhypothese wahr ist (d.h. Männer und Frauen in der Grundgesamtheit die gleichen Gehaltsvorstellungen haben) und wenn alle erforderlichen Annahmen erfüllt sind (u.a., dass es sich um einen echte i.i.d. Zufallsstichprobe handelt), würden wir bei wiederholten Stichprobenziehungen in ca. 7% der Fälle einen ähnlich extremen oder noch extremeren Wert der Teststatistik erwarten. Da dies größer ist als die üblicherweise geforderten 5% werden wir die Nullhypothese nicht verwerfen, und können deshalb von keinem systematischen Zusammenhang ausgehen.

Körpergröße von Söhnen und Vätern

Als nächstes wollen wir die Körpergröße von Söhnen (Kgr) und Vätern (Kgr_vater) vergleichen. Da einige Eingaben in Metern erfolgten müssen wir zuerst diese Eingaben korrigieren.

Mit Hilfe der Indexfunktion in eckigen Klammern können wir die Auswahl auf Männer einschränken.

mydf$Kgr_vater = ifelse(mydf$Kgr_vater < 2.5, mydf$Kgr_vater*100, mydf$Kgr_vater)

boxplot(mydf$Kgr[mydf$Geschlecht == "männlich"], mydf$Kgr_vater[mydf$Geschlecht == "männlich"], names = c("Söhne", "Väter"))

t.test(mydf$Kgr[mydf$Geschlecht == "männlich"], mydf$Kgr_vater[mydf$Geschlecht == "männlich"])
## 
##  Welch Two Sample t-test
## 
## data:  mydf$Kgr[mydf$Geschlecht == "männlich"] and mydf$Kgr_vater[mydf$Geschlecht == "männlich"]
## t = 0.8569, df = 67.495, p-value = 0.3945
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.698210  4.253765
## sample estimates:
## mean of x mean of y 
##  180.5278  179.2500
t.test(mydf$Kgr[mydf$Geschlecht == "männlich"], mydf$Kgr_vater[mydf$Geschlecht == "männlich"], paired = TRUE)
## 
##  Paired t-test
## 
## data:  mydf$Kgr[mydf$Geschlecht == "männlich"] and mydf$Kgr_vater[mydf$Geschlecht == "männlich"]
## t = 1.3175, df = 35, p-value = 0.1962
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.691168  3.246724
## sample estimates:
## mean of the differences 
##                1.277778

Gibt es eine Assoziation zwischen Gehaltsvorstellungen und der Körpergröße?

Zahlreiche empirische Untersuchungen finden einen positiven Zusammenhang zwischen Körpergröße und Einkommen (z.B. Schick and Steckel 2015). Wir könnten uns fragen, ob es auch einen Zusammenhang zwischen den Einkommenserwartungen Studierender und deren Körpergröße gibt.

Da sich die durchschnittliche Körpergröße von Männern und Frauen unterscheidet ist es wichtig auch für das Geschlecht zu kontrollieren, um einen omitted variables bias zu vermeiden.

eq1 <- lm(Eink_erw ~ Kgr + Geschlecht, data = mydf)
summary(eq1)
## 
## Call:
## lm(formula = Eink_erw ~ Kgr + Geschlecht, data = mydf)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1173.52  -270.24   -58.06    76.91  1260.99 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)         399.978   1639.480   0.244    0.808
## Kgr                   8.628      9.064   0.952    0.345
## Geschlechtweiblich  -84.562    153.019  -0.553    0.582
## 
## Residual standard error: 444.3 on 67 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.05744,    Adjusted R-squared:  0.02931 
## F-statistic: 2.042 on 2 and 67 DF,  p-value: 0.1378

Allerdings würde man diesen eher unübersichtlichen Output nie direkt publizieren, man könnte z.B. wieder das stargazer package verwenden

# library(stargazer) # schon vorher aufgerufen
stargazer(eq1, dep.var.labels = "Einkommenserwartungen", dep.var.caption = "",
          covariate.labels = c("Körpergröße", "weiblich", "Interzept"), 
          type = "html", omit.stat = c("f", "ser"), header = FALSE)
Einkommenserwartungen
Körpergröße 8.628
(9.064)
weiblich -84.562
(153.019)
Interzept 399.978
(1,639.480)
Observations 70
R2 0.057
Adjusted R2 0.029
Note: p<0.1; p<0.05; p<0.01

Mit dieser Regression versuchen wir zu überprüfen, ob Einkommenserwartungen einer Gruppe Studierender mit deren Körpergröße korreliert ist. Aufgrund früherer Untersuchungen (u.a. Schick and Steckel 2015) erwarten wir einen positiven Zusammenhang.

Allerdings ist keiner der Koeffizienten statistisch signifikant von Null verschieden (alle p-Werte sind weit größer als 0.05). Auch das Bestimmtheitsmaß ist sehr klein.

Deshalb werden wir die Nullhypothese nicht verwerfen, wir finden also keine empirische Evidenz für einen systematischen Zusammenhang zwischen Einkommenserwartungen und Körpergröße.

Einen Scatterplot (Streudiagramm) mit Regressionsgrade erhalten Sie mit

plot(mydf$Eink_erw ~ mydf$Kgr, main = "Einkommenserwartungen und Körpergröße", 
    xlab = "Körpergröße (cm)", ylab = "Eink.-Erwartungen (Euro)")
abline(lm(mydf$Eink_erw ~ mydf$Kgr), col = "blue", lty = 2, lwd = 2)

 

Beispiel: Stundenlöhne in Österreich


Im diesem Beispiel werden die Brutto-Stundenlöhne unselbständig Beschäftigter in Österreich (EU-Silc 2018, Statistik Austria) in Abhängigkeit von der potentiellen Bildungsdauer (in Jahren) dargestellt.

Die potentiellen Bildungsdauer (potBildg) ist definiert als Alter bei höchstem Bildungsabschluss minus sechs.

Getrennte Regressionen der logarithmierten Stundenlöhne auf die potentielle Bildung für Männer und Frauen ergeben

## Stundenlöhne, Eu-Silc 2018
## Datum: Sommer 2020 

rm(list = ls())
silc <- read.csv2("https://www.uibk.ac.at/econometrics/data/stdl_silc2018.csv")

# Ausreisser
silc$StdL[silc$StdL < 1 | silc$StdL > 200] <- NA 

###############################################
# Regression: Stundenlöhne und pot. Bildung
###############################################

eq_w <- lm(log(StdL) ~ potBildg, data = silc, subset = weibl == 1)
eq_m <- lm(log(StdL) ~ potBildg, data = silc, subset = weibl == 0)

# library(stargazer)
stargazer(eq_w, eq_m, type = "html",
          dep.var.labels = "log(Stundenlöhne)",
          column.labels = c("weibl.", "männl."),
          intercept.bottom = FALSE, 
          omit.stat = c("adj.rsq", "f", "ser"),
          notes = "Daten: EU-Silc 2018 (Statistik Austria)"
          )
Dependent variable:
log(Stundenlöhne)
weibl. männl.
(1) (2)
Constant 2.395*** 2.475***
(0.025) (0.029)
potBildg 0.018*** 0.023***
(0.001) (0.002)
Observations 1,709 1,562
R2 0.081 0.106
Note: p<0.1; p<0.05; p<0.01
Daten: EU-Silc 2018 (Statistik Austria)

 

Die Grafik wird mit Hilfe von low-level Befehlen aufgebaut

###############################################
# Abbildung: Stundenlöhne und pot. Bildung
###############################################

# x11()
plot.new()
plot.window(xlim = c(0, 51),
            ylim = c(0, 5.5), xaxs = "i")
# box()
axis(1, col.axis = "grey30")
axis(2, col.axis = "grey30", las = 1)
title(main = "log Stundenlöhne & Bildung für Österreich", cex.main = 0.9,
      # col.main = "green4",
      sub = "(Datenquelle: EU-Silc 2018)", cex.sub = 0.8,
      # col.sub = "green4",
      xlab = "pot. Bildung (Jahre)", 
      ylab = "log(Stundenlöhne)", 
      col.lab = "blue", font.lab = 3)
points(jitter(x = silc$potBildg[silc$weibl == 0],factor=2), 
       y = log(silc$StdL[silc$weibl == 0]), 
       pch=24, col="blue", cex=0.7)
points(jitter(x = silc$potBildg[silc$weibl == 1],factor=2), 
       y = log(silc$StdL[silc$weibl == 1]), 
       pch=21, col="red", cex=0.8)
# getrennte Regressionen fuer Maenner und Frauen
abline(lm(log(StdL) ~ potBildg, data = silc, subset = weibl == 0), col="blue", lwd=2, lty=3)
abline(lm(log(StdL) ~ potBildg, data = silc, subset = weibl == 1), col="red",  lwd=2, lty=4)
abline(v = mean(silc$potBildg[silc$weibl == 0 & !is.na(silc$StdL)], na.rm = TRUE), lty=1, col="blue")
abline(v = mean(silc$potBildg[silc$weibl == 1 & !is.na(silc$StdL)], na.rm = TRUE), lty=2, col="red")
legend("topright", legend=c("Männer", "Frauen"),
       col=c("blue", "red"), pch=c(24,21), cex=1)
text(x = mean(silc$potBildg[silc$weibl == 1], na.rm = TRUE), 
     y = 5.4, 
     labels = "Mittelw. von Bildg. \nFrauen & Männer", 
     cex = 0.7, pos = 4)

 

 

Literatur

Braun, W. John, and Duncan J. Murdoch. 2008. A First Course in Statistical Programming with R. 1st ed. USA: Cambridge University Press.

Schick, Andreas, and Richard H. Steckel. 2015. “Height, Human Capital, and Earnings: The Contributions of Cognitive and Noncognitive Ability.” Journal of Human Capital 9 (1): 94–115. https://doi.org/10.1086/679675.

Schwabish, Jonathan A. 2014. “An Economist’s Guide to Visualizing Data.” Journal of Economic Perspectives 28 (1): 209–34. https://doi.org/10.1257/jep.28.1.209.

  1. Das R Base Package wird automatisch mit R geladen und umfasst die wesentlichen Funktionen der Programmiersprache R.↩︎