Kapitel 8 Statistik mit R und RStudio

Ebenso wie die Statistik, ist dieses Kapitel in folgende Bereiche unterteilt:

  1. Deskriptive Statistik (Abschnitt 8.1)

  2. Schließende Statistik (Abschnitt 8.2)

Die Funktionen sind alle in klassischem R-Code dargestellt, da sie so auch in einer Tidyverse-Befehlskette stehen würden.

Zu Demonstrationszwecken werden wir in diesem Kapitel auf Beispieldatensätze epa, mma , nw und pf8 zurückgreifen, die Sie sich aus dem Internet in Ihre R-Session laden können.

# lade Datensatz "epa"
load(url("https://www.produnis.de/R/epa.RData"))
# lade Datensatz "mma"
load(url("https://www.produnis.de/R/mma.RData"))
# lade Datensatz "pf8"
load(url("https://www.produnis.de/R/pf8.RData"))
# lade Datensatz "nw"
load(url("https://www.produnis.de/R/nw.RData"))

8.1 Deskriptive Statistik

Zur Beschreibung von Verteilungen eigenen sich Tabellen, Lage- und Streuungskenngrößen.

8.1.1 Häufigkeiten

Eine Häufigkeitstabelle wird mit der Funktion table() erstellt

# Daten erzeugen
x <- c(1, 3, 5, 7, 7, 7, 6, 34)

# Häufigkeitstabelle
table(x)
## x
##  1  3  5  6  7 34 
##  1  1  1  1  3  1

Für relative Häufigkeiten teilt man einfach durch die Länge der Variable x. So ändert sich der Befehl in:

# relative Häufigkeitstabelle
table(x)/length(x)
## x
##     1     3     5     6     7    34 
## 0.125 0.125 0.125 0.125 0.375 0.125

Kumulierte Häufigkeiten erhält man mit der Funktion cumsum()

# kumulierte Häufigkeitstabelle
cumsum(table(x))
##  1  3  5  6  7 34 
##  1  2  3  4  7  8

Und entsprechend auch kumulierte relative Häufigkeiten

# kumulierte relative Häufigkeitstabelle
cumsum(table(x)/length(x))
##     1     3     5     6     7    34 
## 0.125 0.250 0.375 0.500 0.875 1.000

Eine weitere Option bietet die Funktion freq() aus dem rettyR-Paket. Sie erstellt eine Tabelle mit absoluten und relativen Häufigkeiten sowie relativen Häufikeiten ohne fehlende Wert (NA).

# Paket "prettyR" installieren
install.packages("prettyR", dependencies=TRUE)
# Paket laden
library(prettyR)

# Häufigkeitstabelle mit absoluten und relative Werten erstellen
freq(x)
## 
## Frequencies for x 
##         7    1    3    5    6   34   NA
##         3    1    1    1    1    1    0
## %    37.5 12.5 12.5 12.5 12.5 12.5    0 
## %!NA 37.5 12.5 12.5 12.5 12.5 12.5

8.1.2 Lagekenngrößen

Mit dem Befehl summary() erhält man eine erste Zusammenfassung der Werteverteilung mit den wichtigsten Lagekenngrößen.

# Daten erzeugen
x <- c(1, 3, 5, 7, 7, 7, 6, 34)

# Zusammenfassung
summary(x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    4.50    6.50    8.75    7.00   34.00

Die Lagewerte lassen sich auch einzeln bestimmen.

Mit der Funktion min() wird das Minimum der Verteilung ausgegeben.

# Minimum bestimmen
min(x)
## [1] 1

Dementsprechend liefert max() das Maximum.

# Maximum bestimmen
max(x)
## [1] 34

Die Spannweite wird mit range() abgefragt.

# Spannweite bestimmen
range(x)
## [1]  1 34

Das arithmetische Mittel wird mit der Funktion mean() bestimmt

# x-quer bestimmen
mean(x)
## [1] 8.75

Sind fehlende Werte in der Reihe enthalten, gibt R nur NA zurück.

#  mean für "Gewicht" im Datensatz "pf8"
mean(pf8$Gewicht)
## [1] NA

Beheben kann man dies, indem der Parameter na.rm auf TRUE gesetzt wird.

#  mean für "Gewicht" im Datensatz "pf8"
# entferne NAs
mean(pf8$Gewicht, na.rm=TRUE)
## [1] 75.5325

Der Median wird mit der Funktion median() ermittelt.

# Median bestimmen
median(x)
## [1] 6.5

Auch hier müssen enthaltene NAs ausgeblendet werden:

# Median bestimmen, OHNE NAs
median(pf8$Gewicht, na.rm=TRUE)
## [1] 73

In den R-Hausmitteln existiert keine eigene Funktion zur Bestimmung des Modalwertes. Wir können ihn daher “auf Umwegen” bestimmen. Wir lassen mittels table() eine Häufigkeitstabelle der Variable ausgeben. Diese sortieren wir mittels sort() absteigend. Der Modus ist dann der erste Wert der sortierten Tabelle.

# Daten als Häufigkeitstabelle
table(x)
## x
##  1  3  5  6  7 34 
##  1  1  1  1  3  1
# absteigend sortieren
sort(table(x), decreasing=TRUE)
## x
##  7  1  3  5  6 34 
##  3  1  1  1  1  1
# Modalwert
sort(table(x), decreasing=TRUE)[1]
## 7 
## 3

Es ist der Wert 7 und er kommt 3 mal vor.

Eine weitere Möglichkeit bietet das Zusatzpaket statip. Es enthält die Funktion mfv(), mit welcher der Modalwert bestimmt werden kann.

# Daten als Häufigkeitstabelle
library(statip)

# Modalwert bestimmen
mfv(x)
## [1] 7

Die vier Quartile werden mit der Funktion quantile() ermittelt (ja, der Befehl für Quartile heisst in R quantile()).

# Quartile bestimmen, ohne NAs
quantile(x, na.rm=T)
##   0%  25%  50%  75% 100% 
##  1.0  4.5  6.5  7.0 34.0

Es lassen sich beliebige Quantile berechnen.

# Quartile bestimmen, ohne NAs
quantile(x, probs = c(.333, .666))
## 33.3% 66.6% 
## 5.331 7.000

8.1.3 Streuungskenngrößen

Die wichtigste Streuungskenngröße ist die Standardabweichung. Sie wird bestimmt mittels der Funktion sd().

# Standardabweichung bestimmen
sd(x, na.rm=T)
## [1] 10.43004

Die Varianz der Stichprobe wird mit var() berechnet.

# Varianz von x bestimmen
var(x, na.rm=T)
## [1] 108.7857

Der Interquartilsabstand (oder einfach nur Quartilsabstand) wird mit der Funktion IQR() bestimmt.

# Quartilsabstand bestimmen
IQR(x, na.rm=T)
## [1] 2.5

8.1.4 Kreuztabellen

Mit der Funktion table() lassen sich auch Kreuztabellen erstellen, indem einfach beide Variablen übergeben werden.

# erzeuge Daten
Punktwert <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)
Geschlecht <- rep(c("maennlich" , "weiblich") ,5)

# Kreuztabelle
table(Punktwert, Geschlecht)
##          Geschlecht
## Punktwert maennlich weiblich
##        2          1        0
##        4          0        1
##        7          1        0
##        9          0        1
##        10         1        0
##        12         0        1
##        13         1        0
##        15         0        1
##        16         1        0
##        19         0        1

Ebenso steht die Funktion xtabs() zur Verfügung.

# Kreuztabelle mit "xtabs()"
xtabs(~ Punktwert + Geschlecht)
##          Geschlecht
## Punktwert maennlich weiblich
##        2          1        0
##        4          0        1
##        7          1        0
##        9          0        1
##        10         1        0
##        12         0        1
##        13         1        0
##        15         0        1
##        16         1        0
##        19         0        1

8.1.5 z-Transformation

Die z-Transformation erfolgt mit der Funktion scale().

# Wertereihe
Punktwerte <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)

# z-Transformation
scale(Punktwerte)
##             [,1]
##  [1,] -1.6183421
##  [2,]  0.2418212
##  [3,] -0.6882604
##  [4,]  0.7998702
##  [5,] -0.1302114
##  [6,] -1.2463094
##  [7,]  0.4278376
##  [8,] -0.3162278
##  [9,]  0.9858866
## [10,]  1.5439356
## attr(,"scaled:center")
## [1] 10.7
## attr(,"scaled:scale")
## [1] 5.375872

8.1.6 Korrelation

Die allgemeine Funktion zur Berechnung von Korrelationen heisst cor(). Die Korrelation nach Pearson berechnet sich mit der Funktion cor() und deren Parameter method="pearson".

# erzeuge Daten
Lesetest <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)
Rechtschreibung <- c(3 , 14 , 9 , 17 , 12 , 4 , 16 , 12 , 18 , 20)

# Pearson Maßkorrelationskoeffizient
cor(Lesetest, Rechtschreibung, method="pearson")
## [1] 0.9884616

Mit der Funktion cor.test() kann direkt ein Siginifkanztest mitgefahren werden

# Pearson Maßkorrelationskoeffizient
# mit Signifikanztest
cor.test(Lesetest, Rechtschreibung, method="pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  Lesetest and Rechtschreibung
## t = 18.458, df = 8, p-value = 7.648e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9502085 0.9973659
## sample estimates:
##       cor 
## 0.9884616

Sollten fehlende Wert (NA) enthalten sein, gibt der Befehl nur ein NA zurück.

# Pearson Maßkorrelationskoeffizient
# mit Datensatz "pf8", da sind NAs enthalten!!!
with(pf8, cor(Alter,Größe, method="pearson"))
## [1] NA

Mit dem Parameter use="complete.obs" (nutze nur vollständige Beobachtungen) können wir dies beheben.

# Pearson Maßkorrelationskoeffizient
# mit Datensatz "pf8", verwende nur "komplette" Daten (ohne NA)
with(pf8, cor(Alter,Größe, method="pearson", use="complete.obs"))
## [1] -0.1067484

Für die Korrelationsberechnung nach Spearmen (Spearman’s \(\rho\)) wird der Parameter method auf spearman gesetzt.

# Spearman's Korrelationsberechnung
cor(Lesetest, Rechtschreibung, method="spearman")
## [1] 0.9969651

Auch hier müssen NAs mit dem Parameter use="complete.obs" unterdrückt werden.

# Spearman's Rho
# mit Datensatz "pf8", verwende nur "komplette" Daten (ohne NA)
with(pf8, cor(Alter,Größe, method="spearman", use="complete.obs"))
## [1] -0.05920716

Auch für Spearman’s Rho kann ein Signifikanztest mitgefahren werden, indem cor.test() aufgerufen wird.

# Spearman's Korrelationsberechnung
# mit Signifikanztest
cor.test(Lesetest, Rechtschreibung, method="spearman")
## Warning in cor.test.default(Lesetest, Rechtschreibung, method = "spearman"):
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  Lesetest and Rechtschreibung
## S = 0.50076, p-value = 3.698e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9969651

Die Warnmeldung kann ignoriert werden, der p-Wert wird in der Ausgabe angezeigt.

Da cor.test() ein Objekt der Klasse Liste zurückliefert, kann auf die einzelnen Komponenten des Tests zugegriffen werden.

# schreibe Testergebnis in Variable "test"
test <- cor.test(Lesetest, Rechtschreibung, method="spearman")

# zeige nur p-Wert an
test$p.value
## [1] 3.698093e-10
# zeige nur Teststatistik an
test$statistic
##         S 
## 0.5007599
# zeige nur Spearmans's Rho an
test$estimate
##       rho 
## 0.9969651

8.1.7 Regressionen

Die Funktion zur Berechnung linearer Regressionen ist lm() (für “linear model”). Sie folgt der Syntax x erklärt durch y, was durch eine Tilde ~ ausgedrückt wird.

# erzeuge Daten
Lesetest <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)
Rechtschreibung <- c(3 , 14 , 9 , 17 , 12 , 4 , 16 , 12 , 18 , 20)

# lineare Regression "Rechtschreibung erklärt durch Lesetest"
lm(Rechtschreibung~Lesetest)
## 
## Call:
## lm(formula = Rechtschreibung ~ Lesetest)
## 
## Coefficients:
## (Intercept)     Lesetest  
##       1.208        1.055

Hierbei ist Intercept der Schnitt durch die Y-Achse.

Wenn Sie über den Parameter data das Datenframe angeben, können Sie auf die Variablen referenzieren, ohne $ wenden zu müssen.

# lade Testdaten
load(url("https://www.produnis.de/R/nw.RData"))

# entweder Variablen per "$" referenzieren
lm(nw$workyear ~ nw$worknight)
## 
## Call:
## lm(formula = nw$workyear ~ nw$worknight)
## 
## Coefficients:
##  (Intercept)  nw$worknight  
##       9.9143        0.8291
# oder Datenframe mit angeben
lm(workyear ~ worknight, data=nw)
## 
## Call:
## lm(formula = workyear ~ worknight, data = nw)
## 
## Coefficients:
## (Intercept)    worknight  
##      9.9143       0.8291

Es ist üblich, Modelle in ein Objekt zu speichern und über die summary()- Funktion aufzurufen. Die Informationen der Modellzusammenfassung sind so detailierter.

# speichere Modell in "fit"
fit <- lm(workyear~worknight, data=nw)
# schaue mit "summary()" an
summary(fit)
## 
## Call:
## lm(formula = workyear ~ worknight, data = nw)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.889  -5.743  -2.351   3.660  26.427 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.91425    0.74088   13.38   <2e-16 ***
## worknight    0.82912    0.05958   13.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.463 on 274 degrees of freedom
## Multiple R-squared:  0.4141, Adjusted R-squared:  0.4119 
## F-statistic: 193.6 on 1 and 274 DF,  p-value: < 2.2e-16

Außerdem lässt sich das Modell so leicht plotten.

plot(fit)

Dies erzeugt insgesamt 4 Plots (Residualplot, Q-Q, Scale-Location und Einfluss). Durch drücken der [Enter]-Taste schalten Sie zum nächsten Plot weiter.

Sie können aber auch alle 4 Plots in einem ausgeben lassen.

# plotte alle 4 auf einer Seite
par(mfrow = c(2, 2), oma = c(0, 0, 2, 0))
plot(fit)

8.1.7.1 multiple lineare Regressionen

Multiple lineare Regressionen werden ebenfalls mit der Funktion lm() berechnet. Hierbei werden die erklärenden Variablen per + Zeichen aneinandergereiht.

# multiple lineare Regression 
# "workyear" erklärt durch "worknight" UND "age"
fit <- lm(workyear ~ worknight+age, data=nw)
summary(fit)#
## 
## Call:
## lm(formula = workyear ~ worknight + age, data = nw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7077  -3.9748  -0.3066   3.5936  19.9569 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.99856    1.72557  -3.476 0.000591 ***
## worknight    0.57707    0.05714  10.100  < 2e-16 ***
## age          0.41614    0.04195   9.921  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.41 on 273 degrees of freedom
## Multiple R-squared:  0.5693, Adjusted R-squared:  0.5662 
## F-statistic: 180.5 on 2 and 273 DF,  p-value: < 2.2e-16

Ebenfalls können hier die 4 Plots erzeugt werden.

# plotte alle 4 auf einer Seite
par(mfrow = c(2, 2), oma = c(0, 0, 2, 0))
plot(fit)

8.1.7.2 Logistische Regression

Als Beispiel dienen Daten der Challenger-Katastrophe (am 28. Januar 1986 explodierte die Raumfähre Challenger 73 Sekunden nach dem Start).

Der Ausfall eines oder mehrerer Dichtungsringe in einer der seitlichen Feststoffrakten wurde als Grund der Katastrophe ermittelt. Probleme mit den Dichtungsringen sind auch vor- her am Boden aufgetreten. So liegen Erhebungsdaten vor, die den Defekt ei- nes Rings in Zusammenhang mit der Außentemperatur darstellen können: Je kälter es draußen ist, desto höher die Ausfallwahrscheinlichkeit eines Rings. Bedauerlicher Weise wurden diese Daten erst rückwirkend ausgewertet. Diese Daten übertragen wir wie folgt in R:

Temp <- c(66, 67, 68, 70, 72, 75, 76, 79, 53, 58, 70, 75, 67, 67, 69, 70, 73, 76, 78, 81, 57, 63, 70)
Ausfall <- factor(c(0, 0,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1 ,1 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1, 1))
shuttle <- data.frame(Temp, Ausfall)
colnames(shuttle) <- c("Temp", "Ausfall")

Die Kodierung des factor Ausfall lautet: 0=kein Ausfall | 1=Ausfall Die Einheit für Temp lautet Fahrenheit.

Der Aufruf zur logistischen Regression (Ausfall erklärt durch Temp) lautet:

fit <- glm(Ausfall~Temp,data=shuttle, family=binomial)
summary(fit)
## 
## Call:
## glm(formula = Ausfall ~ Temp, family = binomial, data = shuttle)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0611  -0.7613  -0.3783   0.4524   2.2175  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  15.0429     7.3786   2.039   0.0415 *
## Temp         -0.2322     0.1082  -2.145   0.0320 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.267  on 22  degrees of freedom
## Residual deviance: 20.315  on 21  degrees of freedom
## AIC: 24.315
## 
## Number of Fisher Scoring iterations: 5

Die Odds Ratio (OR) für \(\beta\)Temp= -0,2322 errechnet sich mit der exp()-Funktion:

exp(-0.2322)
## [1] 0.7927875

Dies liefert das Ergebnis OR = 0.7927875. Die Interpretation lautet: Mit jedem höhren Grad an Temparatur ändert sich die Chance auf einem Ausfall um den Faktor 0.7927. Das bedeutet, dass die Chancen eines Ausfalls niedriger werden, je höher die Außentemperatur steigt.

In der Nacht vor dem Challengerstart konnte eine Temperatur von 37° Fahrenheit gemessen werden. Auf Grundlage der logistischen Regression lässt sich nun die Ausfallwahrscheinlichkeit bei 37° Außentemperatur berechnen:

predict(fit, data.frame(Temp=37),type="response")
##         1 
## 0.9984265

Das Modell sagt bei einer Außentemperatur von 37° eine Ausfallwahrscheinlichkeit von 99,84% voraus.

Jetzt kann man dies für jede mögliche Temperatur durchmachen, und schon erhält man dieses schöne Plot:

y <- predict(fit, data.frame(Temp=30:90), type="response")
plot(30:90, y, type="l", ylab="Ausfallswahrscheinlichkeit", xlab="Temperatur in F")
abline(v=seq(20,100,1), h=seq(-0.2,1.2,0.05),col="gray") #Mathegitter
abline(v=37, col="red") # 37-Grad
abline(h=0.5, col="black") # Mittellinie
points(30:90, y, type="l",lwd=2, col="blue") # Schätzwerte in "blau"

8.1.7.3 Ordinale Regression

Zur Bestimmung des Einflusses von Prädiktoren auf ein ordinales Ziel stehen verschiedene Modelle der ordinalen Regression zur Verfügung. Namentlich werden hier vor allem Proportional Odds Modelle und Continuation Ratio Modelle verwendet.

Eine verständliche Einführung in ordinale Modelle findet sich im Open Access Artikel von große Schlarmann und Galatsch6.

Für die ordinalen Regressionen verwenden wir das “VGAM,” welches zunächst in R installiert werden muss

install.packages("VGAM", dependencies=T)

Als Beispiel dient der folgende Datensatz:

load(url("http://www.produnis.de/R/OrdinalSample.RData"))
mydata <- ordinalSample
head(mydata)
##    Konflikt Zufriedenh Geschlecht Stimmung
## 1       4.2       2.50          1  maessig
## 6       1.8       3.00          1 sehr gut
## 15      3.4       1.75          1  maessig
## 16      4.0       2.50          2  maessig
## 22      2.2       2.50          1      gut
## 23      3.4       2.25          1  maessig
  • Die Variable “Stimmung” dient hier als ordinales Ziel mit den Ausprägungen “schlecht” < “maessig” < “gut” < “sehr gut.”

  • Die Variable “Konflikt” dient als Prädiktor. Sie gibt an, wieviele Konflikte derzeit im Arbeitsleben vorliegen.

  • Die Variable “Zufriedenh” beschreibt, wie zufrieden Probanden mit ihrem Job sind.

  • Die Variable “Stimmung” soll nun durch “Konflikt” und “Zufriedenh” beschrieben werden.

Wir wollen wissen, wie stark Stimmung durch Konflikt und Zufriedenh beinflusst wird. Das heisst, in unserem Modell soll die Variable Stimmung durch Konflikt und Zufriedenh beschrieben werden.

8.1.7.3.1 Proportional Odds Modell

Ein Proportional Odds Modell errechnet sich leicht mit dem VGAM-Paket und der Funktion vglm():

library(VGAM)
pom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T))

# oder abgekürzt
pom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=propodds)
summary(pom)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = propodds, 
##     data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1   0.6210     0.6482   0.958   0.3381    
## (Intercept):2  -1.7703     0.6540  -2.707   0.0068 ** 
## (Intercept):3  -4.0578     0.6759  -6.004 1.93e-09 ***
## Konflikt       -0.5762     0.0970  -5.940 2.86e-09 ***
## Zufriedenh      1.3643     0.2021   6.751 1.47e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3]), 
## logitlink(P[Y>=4])
## 
## Residual deviance: 920.5493 on 1240 degrees of freedom
## 
## Log-likelihood: -460.2746 on 1240 degrees of freedom
## 
## Number of Fisher scoring iterations: 5 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
##   Konflikt Zufriedenh 
##  0.5620549  3.9129275

Mit dem Modell können nun weitere Parameter errechnet werden:

# Koeffizienten
pom.coef <- (coef(summary(pom)));pom.coef
##                 Estimate Std. Error    z value     Pr(>|z|)
## (Intercept):1  0.6209962 0.64821816  0.9580049 3.380603e-01
## (Intercept):2 -1.7702692 0.65404383 -2.7066523 6.796539e-03
## (Intercept):3 -4.0577657 0.67586741 -6.0037895 1.927645e-09
## Konflikt      -0.5761558 0.09700358 -5.9395317 2.858375e-09
## Zufriedenh     1.3642858 0.20208974  6.7508909 1.469400e-11

Die Odds Ratio errechnet sich mit der exp()-Funktion:

# Odds Ratio
pom.odds <- exp(coef(pom));pom.odds
## (Intercept):1 (Intercept):2 (Intercept):3      Konflikt    Zufriedenh 
##     1.8607808     0.1702871     0.0172876     0.5620549     3.9129275
# Devianz und AIC
pom.devi <- deviance(pom);pom.devi
## [1] 920.5493
pom.aic  <- AIC(pom);pom.aic
## [1] 930.5493
# logLikelihood
pom.ll <- logLik(pom);pom.ll
## [1] -460.2746
# 0-Modell (fuer pseudo R^2)
p0 <- vglm(Stimmung ~ 1, data=mydata, family=propodds)
p0.ll <- logLik(p0)
# R^2 McFadden
pom.mcfad <- as.vector(1- (pom.ll/p0.ll));pom.mcfad
## [1] 0.1240613
# R^2 Cox&Snell
N <- length(mydata[,1]) # Anzahl der Fälle
pom.cox <- as.vector(1 - exp((2/N) * (p0.ll - pom.ll)));pom.cox
## [1] 0.2696034
# R^2 Nagelkerke
pom.nagel <- as.vector((1 - exp((2/N) * (p0.ll - pom.ll))) / (1 - exp(p0.ll)^(2/N))) ;pom.nagel
## [1] 0.2928789

Das Proportional Odds Modell geht von der “equal slopes assumption” (auch: “proportional odds assumption”) aus. Diese Annahme muss geprüft werden, bevor das Modell als gültig angesehen werden kann. Dies geht mit dem VGAM-Paket recht einfach. Der Befehl vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=propodds) ist eine Abkürzung für vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T)). Der Parameter parallel=TRUE/FALSE stellt ein, ob das Modell mit equal slopes (=TRUE), oder ohne equal slopes assumption (=FALSE) erstellt werden soll. Zur Überprüfung der “Equal Slopes Assumption” erstellt man jeweils ein TRUE und ein FALSE-Modell, und führt dann einen Likelihood-Test durch. Die Nullhypothese lautet, dass equal slopes vorliegen.

# Modell OHNE equal slopes
npom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=F))

# log-likelihood-Test auf Equal Slopes Assumption
lrtest(pom,npom)# Test
## Likelihood ratio test
## 
## Model 1: Stimmung ~ Konflikt + Zufriedenh
## Model 2: Stimmung ~ Konflikt + Zufriedenh
##    #Df  LogLik Df  Chisq Pr(>Chisq)
## 1 1240 -460.27                     
## 2 1236 -456.46 -4 7.6202     0.1065

Der Test ist mit p=0,1065 nicht signifikant. Das heisst, es liegen equal slopes vor. Das Modell darf also verwendet werden.

Alternativ kann dieser Test auch auf die Devianz der Modelle bestimmt werden:

# Devianz-Test auf Equal Slopes Assumption
pom.pdevi = (1 - pchisq(deviance(pom) - deviance(npom), df=df.residual(pom)-df.residual(npom)));pom.pdevi
## [1] 0.106526

Ist der Test signifikant, liegen keine “equal slopes” vor. Hier kann ein Partial Proportional Odds Modell erstellt werden, welches die “equal slopes” nur für bestimmte Prädiktoren annimmt. Mit dem VGAM-Paket kann recht einfach bestimmt werden, wie ein solches Modell erstellt werden soll. Hierfür erweitertet man den “parallel”-Switch wie folgt:

# Equal Slopes NUR für "Konflikt"
ppom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T~Konflikt-1));ppom
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cumulative(parallel = T ~ 
##     Konflikt - 1), data = mydata)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3      Konflikt  Zufriedenh:1 
##     0.6614816     1.5373699     2.7337389     0.5705897    -1.9684047 
##  Zufriedenh:2  Zufriedenh:3 
##    -1.2805318    -0.8865225 
## 
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.4929 
## Log-likelihood: -457.2464
# Nochmals EqualSlopes NUR für "Konflikt"
ppom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=F~Zufriedenh));ppom
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cumulative(parallel = F ~ 
##     Zufriedenh), data = mydata)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3      Konflikt  Zufriedenh:1 
##     0.6614816     1.5373699     2.7337389     0.5705897    -1.9684047 
##  Zufriedenh:2  Zufriedenh:3 
##    -1.2805318    -0.8865225 
## 
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.4929 
## Log-likelihood: -457.2464
  • parallel=T~Konflikt-1 bedeutet, dass equal slopes nur für Konflikt angenommen wird.

  • parallel=F~Zufriedenh bedeutet, dass equal slopes nur für Zufriedenh nicht angenommen wird.

Beide Befehle bedeuten also das selbe. Daher ist die R-Ausgabe bei beiden Befehlen gleich.

Eine Koeffizientenübersicht erhält man per

ppom.ce <- coef(summary(ppom));ppom.ce
##                 Estimate Std. Error    z value     Pr(>|z|)
## (Intercept):1  0.6614816 0.85866412  0.7703613 4.410856e-01
## (Intercept):2  1.5373699 0.71997824  2.1353005 3.273647e-02
## (Intercept):3  2.7337389 0.98084390  2.7871294 5.317723e-03
## Konflikt       0.5705897 0.09715456  5.8730094 4.279541e-09
## Zufriedenh:1  -1.9684047 0.34535219 -5.6997024 1.200167e-08
## Zufriedenh:2  -1.2805318 0.23585229 -5.4293805 5.654999e-08
## Zufriedenh:3  -0.8865225 0.32542569 -2.7241933 6.445877e-03

Mit einem kleinen Script können Konfidenzintervalle und andere Werte ausgegeben werden:

get.ci <- function(x){
  back <-cbind( x[1], # estimate
                (x[1] - (1.96 * x[2])), #lCI
                (x[1] + (1.96 * x[2])), #uCI
                x[2], x[3], # SD und z
                (2*(1 -pnorm(abs(x[3]))) ),# p-wert
                exp(x[1])) 
  colnames(back) <- c("Estimate", "lCI", "uCI","SD","z","p-value", "OR")
  return(back)
}

Der Aufruf erfolgt z.B. so:

get.ci(as.numeric(ppom.ce[4,]))
##       Estimate       lCI       uCI         SD        z      p-value      OR
## [1,] 0.5705897 0.3801667 0.7610126 0.09715456 5.873009 4.279541e-09 1.76931
8.1.7.3.2 Continuation Ratio Model

Um ein Continuation Ratio Modell zu rechnen, muss der Parameter family auf cratio gesetzt werden.

# CR-Modell MIT PO-Assumption  (Vorwärts)
crm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=T,reverse=F))
summary(crm)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = T, 
##     reverse = F), data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1  0.92610    0.57405   1.613   0.1067    
## (Intercept):2 -1.14132    0.57924  -1.970   0.0488 *  
## (Intercept):3 -3.01660    0.60706  -4.969 6.72e-07 ***
## Konflikt      -0.53645    0.08639  -6.209 5.32e-10 ***
## Zufriedenh     1.16010    0.17944   6.465 1.01e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y>1|Y>=1]), logitlink(P[Y>2|Y>=2]), 
## logitlink(P[Y>3|Y>=3])
## 
## Residual deviance: 924.3824 on 1240 degrees of freedom
## 
## Log-likelihood: -462.1912 on 1240 degrees of freedom
## 
## Number of Fisher scoring iterations: 6 
## 
## No Hauck-Donner effect found in any of the estimates

Dies berechnet standardmäßig die Vorwärts-Methode. Möchte man die Rückwärts-Methode rechnen, setzt man den Parameter reverse auf TRUE.

# CR-Modell MIT PO-Assumption (Rückwärts)
crm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=T,reverse=TRUE))
summary(crm)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = T, 
##     reverse = TRUE), data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1  0.12822    0.57259   0.224 0.822809    
## (Intercept):2  2.08751    0.58094   3.593 0.000326 ***
## (Intercept):3  4.04245    0.60636   6.667 2.62e-11 ***
## Konflikt       0.45234    0.08488   5.329 9.88e-08 ***
## Zufriedenh    -1.25631    0.18108  -6.938 3.98e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3]), 
## logitlink(P[Y<4|Y<=4])
## 
## Residual deviance: 920.4627 on 1240 degrees of freedom
## 
## Log-likelihood: -460.2314 on 1240 degrees of freedom
## 
## Number of Fisher scoring iterations: 5 
## 
## No Hauck-Donner effect found in any of the estimates

Hat man die passende Methode gewählt, folgen die selben Befehle wie beim Proportional Odds Modell. Zunächst muss überprüft werden, ob “equal slopes” vorliegen:

# CR-Modell OHNE PO-Assumption (Rückwärts)
ncrm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=F,reverse=T));summary(ncrm)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = F, 
##     reverse = T), data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1   2.3649     1.1740   2.014 0.043969 *  
## (Intercept):2   1.9998     0.8117   2.464 0.013753 *  
## (Intercept):3   2.4670     1.1190   2.205 0.027474 *  
## Konflikt:1      0.1465     0.1811   0.809 0.418633    
## Konflikt:2      0.4864     0.1194   4.073 4.65e-05 ***
## Konflikt:3      0.6308     0.1689   3.736 0.000187 ***
## Zufriedenh:1   -1.8008     0.3919  -4.596 4.32e-06 ***
## Zufriedenh:2   -1.2602     0.2578  -4.889 1.01e-06 ***
## Zufriedenh:3   -0.8351     0.3382  -2.469 0.013543 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3]), 
## logitlink(P[Y<4|Y<=4])
## 
## Residual deviance: 914.7924 on 1236 degrees of freedom
## 
## Log-likelihood: -457.3962 on 1236 degrees of freedom
## 
## Number of Fisher scoring iterations: 6 
## 
## No Hauck-Donner effect found in any of the estimates
# loglikelihood test auf equal slopes assumption
lrtest(ncrm,crm)# Test
## Likelihood ratio test
## 
## Model 1: Stimmung ~ Konflikt + Zufriedenh
## Model 2: Stimmung ~ Konflikt + Zufriedenh
##    #Df  LogLik Df  Chisq Pr(>Chisq)
## 1 1236 -457.40                     
## 2 1240 -460.23  4 5.6703     0.2252

Der Test ist nicht signifikant (p=0,23). Das bedeutet, dass die Annahme der equal slopes für die hier getesteten Rückwärts-Modelle beibehalten werden kann.

Nun können weitere Modellparameter bestimmt werden:

# 0-Modell (fuer pseudo R^2)
c0 <- vglm(Stimmung ~ 1, data=mydata, family=cratio(parallel=T))
c0.ll <- logLik(c0)
crm.ll <- logLik(crm)

# R^2 Nagelkerke
crm.nagel <- as.vector((1 - exp((2/N) * (c0.ll - crm.ll))) / (1 - exp(c0.ll)^(2/N))) ;crm.nagel
## [1] 0.2930443
# Devianz
crm.devi <- deviance(crm); crm.devi
## [1] 920.4627
# AIC
crm.aic  <- AIC(crm); crm.aic
## [1] 930.4627

Mit unserer Funktion von oben können wir uns die Modellparameter der Koeffizienten ausgeben lassen.

# Konfidenzinztervalle
crm.ce <- coef(summary(crm));crm.ce
##                 Estimate Std. Error    z value     Pr(>|z|)
## (Intercept):1  0.1282227 0.57259420  0.2239329 8.228095e-01
## (Intercept):2  2.0875126 0.58094099  3.5933299 3.264788e-04
## (Intercept):3  4.0424508 0.60636000  6.6667504 2.615293e-11
## Konflikt       0.4523353 0.08488357  5.3288915 9.881398e-08
## Zufriedenh    -1.2563060 0.18108066 -6.9378255 3.981811e-12
get.ci(as.numeric(crm.ce[4,]))
##       Estimate       lCI       uCI         SD        z      p-value       OR
## [1,] 0.4523353 0.2859635 0.6187071 0.08488357 5.328892 9.881398e-08 1.571979
get.ci(as.numeric(crm.ce[5,]))
##       Estimate       lCI        uCI        SD         z      p-value        OR
## [1,] -1.256306 -1.611224 -0.9013879 0.1810807 -6.937825 3.981704e-12 0.2847038

8.2 Schließende Statistik

8.2.1 Konfidenzintervalle

Konfidenzintervalle folgen der Logik, dass zum Punktschätzwert ein gewisser Fehlerbereich addiert und subtrahiert wird. Hierdurch wird ein Intervall mit Unter- und Obergrenze aufgezogen, das mit einer bestimmten Irrtumswahrscheinlichkeit den “wahren Wert” der Grundgesamtheit enthält.

8.2.1.1 Normalverteilung

Um ein Konfidenzintervall für einen Mittelwert aus einer normalverteilten Wertemenge zu berechnen wird folgende Formel angewendet:

\[\begin{equation} \label{formel:konfidenznormal} \begin{array}{l l} \mu_{u}\\\mu_{o}\\ \end{array} \Big\}= \bar{x} \pm t \cdot \frac{s}{\sqrt{n}} \end{equation}\]

  • \(\bar{x}\) = arithmetisches Mittel

  • \(n\) = Stichprobenumfang

  • \(s\) = Standardabweichung

  • \(t\) = Wert der \(t\)-Verteilung, (mit Freiheitsgraden \(df = n-1\))

Nehmen wir an, eine Zufalls-Erhebung über die jährlichen Ausgaben für Pfle- gehilfsmittel ergab bei 100 Befragten (n = 100) einen Mittelwert \(\bar{x}\)=400,- bei einer Standardabweichung von s=20,-. Das Konfidenzintervall für \(\alpha\)=0.05 wird nun wie folgt bestimmt.

# schreibe gegebene Werte in Variablen
n     <- 100
xquer <- 400
s     <- 20

# 95%KI-Fehler berechnen mit (alpha/2) 
fehler95 <- qt(0.975, df=n-1)*s/sqrt(n)

# Fehlerwert ausgeben
fehler95
## [1] 3.968434

Beachten Sie, dass in der Funktion qt() zur Bestimmung des kritischen t-Wertes die Irrtumswahrscheinlichkeit \(\alpha\) halbiert werden muss, da es sich um eine zweiseitige Fragestellung handelt. So muss qt() nicht mit dem Wert 0.95 (\(\alpha=0.05\)) aufgerufen werden, sondern mit 0.975 (\(\frac{\alpha}{2}=0.025\)).

# untere Grenze des 95%KI
xquer - fehler95
## [1] 396.0316

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 396.0316 Euro.

# obere Grenze des 95%KI
xquer + fehler95
## [1] 403.9684

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 403.9684 Euro.

# Länge des 95%KI
(xquer + fehler95) - (xquer - fehler95)
## [1] 7.936868
# oder schlauer
2*fehler95
## [1] 7.936868

Das 95%-Konfidenzintervall hat eine Länge von 7.936868 Euro.

Für das 99% Konfidenzintervall ändern wir entsprechend um. Auch hier muss, da es sich um eine zweiseitige Fragestellung handelt, mit \(\frac{\alpha}{2}\), also 0.005, gerechnet werden. In der Funktion qt() ändert sich also der Wert auf 0,995.

# 99%KI-Fehler berechnen mit (alpha/2) 
fehler99 <- qt(0.995, df=n-1)*s/sqrt(n)

# zeige fehler99 an
fehler99
## [1] 5.252811
# untere Grenze des 99%KI
xquer - fehler99
## [1] 394.7472

Die untere Grenze des 99%-Konfidenzintervalls liegt bei 394.7472 Euro.

# obere Grenze des 99%KI
xquer + fehler99
## [1] 405.2528

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 405.2528 Euro.

# Länge des 99%KI
(xquer + fehler99) - (xquer - fehler99)
## [1] 10.50562
# oder schlauer
2*fehler99
## [1] 10.50562

Das 99%-Konfidenzintervall hat eine Länge von 10.50562 Euro.

8.2.1.1.1 Unterschied von Mittelwerten

Das Konfidenzintervall für den Unterschied von Mittelwerten normalverteilter Daten berechnet sich mit der Formel

\[\begin{equation} \begin{array}{l l} \mu_{u}\\\mu_{o}\\ \end{array} \Big\}= \hat{d} \pm t \cdot \sqrt{\frac{s_{pool}^2}{n_x}+\frac{s_{pool}^2}{n_y}} \end{equation}\]

  • \(\hat{d}\) = Differenz der Mittelwerte \(\bar{x}\),\(\bar{y}\)

  • \(n_{x}, n_{y}\) = die Umfänge der beiden Stichproben

  • \(s_{pool}\) = gepoolte Standardabweichungen = \(\frac{(n_x-1)\cdot s_x^2\ +\ (n_y-1)\cdot s_y^2}{n_x + n_y -2}\)

  • \(t\) = Wert der \(t\)-Verteilung (mit Freiheitsgraden \(df = n_x+n_y -1\))

Nehmen wir an, in einer Studie wurde die Wirkung von zwei Medikamenten auf den Magnesiumgehalt des Blutserums (in mg/dl) untersucht. 13 Personen erhielten Medikament A und 10 Personen Medikament B. Nach 6 Wochen wurden folgende Daten bestimmt:

-\(\bar{x}_{A}\) = 2,22 mg/dl

  • \(\bar{x}_{B}\) = 2,70 mg/dl

  • \(s_{A}\) = 0,255 mg/dl

  • \(s_{B}\) = 0,306 mg/dl

  • \(n_{A}\) = 13

  • \(n_{B}\) = 10

Das Konfidenzintervall für den Unterschied der Mittelwerte bei \(\alpha=0,05\) wird nun wie folgt bestimmt (Achtung, da es sich um eine zweiseitige Fragestellung handelt, müssen wir mit \(\frac{\alpha}{2} =0,025\) rechnen).

# erzeuge die Daten
xa <- 2.22
xb <- 2.7
sa <- 0.255
sb <- 0.306
na <- 13
nb <- 10

# Differenz der Mittelwerte
d   <- xb - xa

# ausgeben
d
## [1] 0.48
# berechne die gepoolten Standardabweichungen
s_pool = ((na-1)*sa^2 + (nb-1)*sb^2) / (na+nb-2)

# ausgeben
s_pool
## [1] 0.07728686
#Fehlerquote berechnen
fehler95 <- qt(0.975, df=na+nb-1) * sqrt(s_pool^2/na + s_pool^2/nb)

# ausgeben
fehler95
## [1] 0.06741865
# untere Grenze des 95%KI
d - fehler95
## [1] 0.4125813

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.4125813 mg/dl.

# obere Grenze des 95%KI
d + fehler95
## [1] 0.5474187

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.5474187 mg/dl.

# Länge des Intervalls
2*fehler95
## [1] 0.1348373

Das 95%-Konfidenzintervall hat eine Länge von 0.1348373 mg/dl.

Für das 99%-Konfidenzintervall wird der Wert für die Funktion qt() auf 0.995 gesetzt (da es sich um eine zweiseitige Fragestellung handelt, müssen wir mit \(\frac{\alpha}{2} = 0,005\) rechnen).

#Fehlerquote berechnen
fehler99 <- qt(0.995, df=na+nb-1) * sqrt(s_pool^2/na + s_pool^2/nb)

# ausgeben
fehler99
## [1] 0.09163373
# untere Grenze des 99%KI
d - fehler99
## [1] 0.3883663

Die untere Grenze des 99%-Konfidenzintervalls liegt bei 0.3883663 mg/dl.

# obere Grenze des 99%KI
d + fehler99
## [1] 0.5716337

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.5716337 mg/dl.

# Länge des Intervalls
2*fehler99
## [1] 0.1832675

Das 99%-Konfidenzintervall hat eine Länge von 0.1832675 mg/dl.

8.2.1.2 Binomialverteilung

Zur Berechnung von Konfidenzintervallen für Anteilswerte wird folgende Formel verwendet.

\[\begin{equation} \label{formel:konfidenzbinomial} \begin{array}{l l} p_{o}\\p_{u}\\ \end{array} \Big\}= k \pm c \cdot \sqrt{\frac{k \cdot (1-k)}{n}} \end{equation}\]

  • \(n\) = Umfang der Stichprobe

  • \(k\) = geschätzter Anteil

  • \(c\) = Quartilswert für Irrtumswahrscheinlichkeit \(\alpha\)

Nehmen wir an, in einer klinischen Studie wurde an 150 Patienten die Dekubitusprävalenz ermittelt. Dabei wurde eine Prävalenz von 35% festgestellt. Die Konfidenzintervalle für den “wahren Wert” der Grundgesamtheit wird wie folgt bestimmt.

# Daten erzeugen
n <- 150
k <- 0.35

# Fehler berechnen mit (alpha/2)!
fehler95 <- qnorm(0.975)*sqrt(k*(1-k)/n)

# anzeigen
fehler95
## [1] 0.07632963

Beachten Sie, dass in der Funktion qnorm() die Irrtumswahrscheinlichkeit \(\alpha\) halbiert werden muss, da es sich um eine zweiseitige Fragestellung handelt. So muss qnorm() nicht mit dem Wert 0.95 (\(\alpha=0.05\)) aufgerufen werden, sondern mit 0.975 (\(\frac{\alpha}{2}=0.025\)).

# untere Grenze des KI
k - fehler95
## [1] 0.2736704

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.2736704 = 27,36704%.

# obere Grenze des KI
k + fehler95
## [1] 0.4263296

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.4263296 = 42,63296%.

# Länge des KI
2*fehler95
## [1] 0.1526593

Das 95%-Konfidenzintervall hat eine Länge von 0.1526593 = 15,26593%.

Für das 99%-Konfidenzintervall muss der Wert der Funktion qnorm() auf 0.995 (\(\frac{\alpha}{2}=0,005\)) angepasst werden.

# Fehler berechnen mit (alpha/2)!
fehler99 <- qnorm(0.995)*sqrt(k*(1-k)/n)

# anzeigen
fehler99
## [1] 0.1003141
# untere Grenze des KI
k - fehler99
## [1] 0.2496859

Die untere Grenze des 99%-Konfidenzintervalls liegt bei 0.2496859 = 24,96859%.

# obere Grenze des KI
k + fehler99
## [1] 0.4503141

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.4503141 = 45,03141%.

# Länge des KI
2*fehler99
## [1] 0.2006283

Das 99%-Konfidenzintervall hat eine Länge von 0.2006283 = 20,06283%.

8.2.1.2.1 Unterschied von Anteilswerten

Zur Berechnung der Konfidenzintervalle bei Unterschieden von Anteils- werten wird die folgende Formel verwendet.

\[\begin{equation} \label{formel:konfidenzunterschiedbinomial} \begin{array}{l l} d_{o}\\d_{u}\\ \end{array} \Big\}= \hat{d} \pm c \cdot \sqrt{\frac{\hat{p}_{1} \cdot (1 - \hat{p}_{1})}{n_{1}} + \frac{\hat{p}_{2} \cdot (1 - \hat{p}_{2})}{n_{2}}} \end{equation}\]

  • \(\hat{d}\) = die Differenz der Anteile

  • \(n_{1}, n_{2}\) = die Umfänge der beiden Stichproben

  • \(\hat{p}_{1}\), \(\hat{p}_{2}\) = die geschätzten Anteile in beiden Stichproben

  • \(c\) = Quartilswert für Irrtumswahrscheinlichkeit \(\alpha\)

Nehmen wir an, in Bundesland A gaben 25 von 100 Befragten an, Zuzahlungen bei Pflegehilfsmitteln leisten zu müssen. In Bundesland B wurden 150 Personen befragt, von denen 60 von Zuzahlungen betroffen waren. Bundesland A hat

dementsprechend einen Anteil von 0,25 (25%) und Bundesland B einen Anteil von 0,40 (40%) Zuzahlern. Der Punktschätzwert für die Differenz beträgt folglich 0,40-0,25 = 0,15.

Die Konfidenzintervalle für die wahre Differenz d kann nun wie folgt bestimmt werden.

# Daten erzeugen
n1 <- 100
p1 <- 0.25
  
n2 <- 150
p2 <- 0.4

# Differenz der Anteile
d <- p2-p1

Beachten Sie, dass in der Funktion qnorm() die Irrtumswahrscheinlichkeit \(\alpha\) halbiert werden muss, da es sich um eine zweiseitige Fragestellung handelt. So muss qnorm() nicht mit dem Wert 0.95 (\(\alpha=0.05\)) aufgerufen werden, sondern mit 0.975 (\(\frac{\alpha}{2}=0.025\)).

# 95%KI-Fehler berechnen mit (alpha/2) 
fehler95 <- qnorm(0.975)*sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)

#anzeigen
fehler95
## [1] 0.1155382
# untere Grenze des 95%KI
d - fehler95
## [1] 0.03446183

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.03446183 = 3,446183%.

# obere Grenze des 95%KI
d + fehler95
## [1] 0.2655382

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.2655382 = 26,55382%.

# Länge des KI
2*fehler95
## [1] 0.2310763

Die Länge des 95%-Konfidenzintervalls liegt bei 0.2310763 = 23,10763%.

Für das 99%-Konfidenzintervall muss der Wert der Funktion qnorm() auf 0.995 (\(\frac{\alpha}{2}=0,005\)) angepasst werden.

# 99%KI-Fehler berechnen mit (alpha/2) 
fehler99 <- qnorm(0.995)*sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)

#anzeigen
fehler99
## [1] 0.1518429
# untere Grenze des 99%KI
d - fehler99
## [1] -0.001842898

Die untere Grenze des 99%-Konfidenzintervalls liegt bei -0.0018429 = -0,18429%.

# obere Grenze des 99%KI
d + fehler99
## [1] 0.3018429

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.3018429 = 30,18429%.

# Länge des KI
2*fehler99
## [1] 0.3036858

Die Länge des 99%-Konfidenzintervalls liegt bei 0.3036858 = 30,36858%.

Wie Sie sehen, liegt die untere Grenze des 99%-Konfidenzintervalls im negativen Bereich. Errechnet sich ein Intervall, welches die 0 einschließt, bedeutet dies, dass der wahre Unterschied auch 0 sein könnte.

8.2.2 Varianzanalysen

8.2.2.1 einfaktoriell (ANOVA)

Eine einfaktorielle Varianzanalyse (ANOVA) wird mit der Funktion aov() durchgeführt. Sie wird nach der Syntax x erklärt durch y aufgerufen. Wie bei anderen Modellen ist es ratsam, die Ergebnisse der Berechnung in ein Objekt zu speichern, und dieses per summary() anzuschauen.

# Lade Testdatensatz
load(url("https://www.produnis.de/R/pf8.RData"))

# erstelle ANOVA  "Alter" erklärt durch "Standort"
fit <- aov(Alter ~ Standort, data=pf8)

# schaue Modellzusammenfassung an
summary(fit)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## Standort      4  10557  2639.3    8.48 1.09e-06 ***
## Residuals   725 225633   311.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness

Das Ergebnis ist signifikant, wir haben unterschiedliche Varianzen im Alter erklärt durch Standort.

Wir wissen aber noch nicht, innerhalb welche Standorte die Unterschiede bestehen. Dazu führen wir post hoc Tests durch, in denen alle Gruppen paarweise verglichen werden. Wegen dieses multiplen Testens muss eine Fehlerkorrektur für Alpha durchgeführt werden. Dazu steht die Funktion pairwise.t.test() zur Verfügung, die mit den selben Angaben aufgerufen wird.

pairwise.t.test(pf8$Alter, pf8$Standort)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pf8$Alter and pf8$Standort 
## 
##           Rheine Münster Bahn   Ladbergen
## Münster   0.5042 -       -      -        
## Bahn      0.4354 1.0000  -      -        
## Ladbergen 0.5432 1.0000  1.0000 -        
## Internet  0.0090 1.8e-05 0.0002 0.0059   
## 
## P value adjustment method: holm

Die gewählte Fehlerkorrektur ist standardmäßig holm. Möchten Sie auf die Bonferroni-Korrekturmethode wechseln, übergeben Sie den Parameter p.adjust="bonferroni".

pairwise.t.test(pf8$Alter, pf8$Standort, p.adjust="bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pf8$Alter and pf8$Standort 
## 
##           Rheine  Münster Bahn    Ladbergen
## Münster   1.00000 -       -       -        
## Bahn      0.72572 1.00000 -       -        
## Ladbergen 1.00000 1.00000 1.00000 -        
## Internet  0.01292 1.8e-05 0.00022 0.00733  
## 
## P value adjustment method: bonferroni

Die signifikanten Unterschiede befinden sich zwischen dem Standort Internet und allen anderen Standorten.

Ein anschauliches Video zur ANOVA in R findet sich bei Youtube: https://www.youtube.com/watch?v=GctGncQrJew

8.2.3 Signifikanztests

8.2.3.1 \(\chi^2\)-Test

Der \(\chi^2\)-Test wird mit der Funktion chisq.test() durchgeführt. Hierfür müssen die Daten als Kreuztabelle vorliegen.

# lade Testdaten
load(url("https://www.produnis.de/R/nw.RData"))

# erstelle Kreuztabelle aus Alter und Geschlecht
xtabelle <- xtabs(~ nw$age + nw$sex)

# Führe Chiquadrat-Test durch
chisq.test(xtabelle)
## Warning in chisq.test(xtabelle): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  xtabelle
## X-squared = 43.886, df = 44, p-value = 0.4765

Der p-Wert ist nicht signifikant, das heisst, das kein Zusammenhang gefunden werden konnte.

Die Warnung bedeutet, dass der p-Wert aus sehr kleinen Zellenhäufigkeiten (kleiner 5) geschätzt wurde, und daher fehlerhaft sein kann (wir sollten daher einen exakten Test durchführen, zum Beispiel den Fisher-Test).

Ein gutes Anleitungsvideo zum \(\chi^2\)-Test findet sich bei Youtube: https://www.youtube.com/watch?v=4lFyRcXoJB8

8.2.3.2 Fisher’s Exakttest

Der exakte Fisher-Test wird mit der Funktion fisher.test() durchgeführt.

# Führe Fisher-Test durch
fisher.test(xtabelle)
Fisher's Exact Test for Count Data with simulated p-value (based on 2000 replicates)
data: nw$sex and nw$age
p-value = 0.3848
alternative hypothesis: two.sided


Der exakte Fisher-Test bestätigt, dass es keinen Unterschied gibt.

8.2.3.3 t-Test

Der t-Test ist wohl einer der bekanntesten Signifikanztests. In R wird er mit der Funktion t.test() ausgeführt.

8.2.3.3.1 Einstichproben t-Test

Beim Einstichproben t-Test wird der Mittelwert der Grundgesamtheit über den Parameter mu mitgegeben. Schauen wir, ob im Nachtwachendatensatz das Alter von einem gedachten Mittelwert 40 in der Grundgesamtheit abweicht.

# lade Testdaten
load(url("https://www.produnis.de/R/nw.RData"))

# Führe t-Test für "Alter" auf Mittelwert 40 durch
t.test(nw$age, mu=40)
## 
##  One Sample t-test
## 
## data:  nw$age
## t = 6.8274, df = 275, p-value = 5.524e-11
## alternative hypothesis: true mean is not equal to 40
## 95 percent confidence interval:
##  43.00908 45.44744
## sample estimates:
## mean of x 
##  44.22826

Der Test ist signifikant, das bedeutet, das Alter im Datensatz nw weicht von einem Mittelwert in der Grundgesamtheit von 40 ab. Es liegt ein Unterschied vor.

Zur Beschreibung der Effektstärke des Unterschieds können wir Cohen’s d berechnen. Dies kann über das Zusatzpaket lsr und der Funktion cohensD() erfolgen.

# aktiviere Zusatzpaket "lsr"
library(lsr)

# berechne Effektsärke Cohen's D
cohensD(nw$age, mu=40)
## [1] 0.4109627

Dieser Wert kann nun nach Cohen’s Angaben7 interpretiert werden (0.2=kleiner Effekt, 0.5=mittel, 0.8=stark).

Unsere Effektgröße ist mit 0.4109627 eher “mittel.”

Ein gutes Anleitungsvideo zum Einstichproben t-Test findet sich bei Youtube: https://www.youtube.com/watch?v=7vcyjuSwfzI

8.2.3.3.2 unabhängiger t-Test

Für den Vergleich der Mittelwerte bei unabhängige Stichproben werden die Werte beider Gruppen als Variable übergeben. So können wir z.B. schauen, ob sich die Datensätze nw und epa hinsichtlich des Alters unterscheiden

# t-test auf Alter in Datensatz "epa" und "nw"
t.test(nw$age, epa$age)
## 
##  Welch Two Sample t-test
## 
## data:  nw$age and epa$age
## t = -19.316, df = 854.37, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -20.65610 -16.84545
## sample estimates:
## mean of x mean of y 
##  44.22826  62.97903

Das Ergebnis ist signifikant, die Gruppen unterscheiden sich also.

Standardmäßig wird ein zweiseitiger Test durchgeführt. Über den Parameter alternative kann dies geändert werden.

# t-test einseitig, kleiner
t.test(nw$age, epa$age, alternative = "less")
# t-test einseitig, größer
t.test(nw$age, epa$age, alternative = "greater")
# t-test zweiseitig
t.test(nw$age, epa$age, alternative = "two.sided")

Die Funktion nimmt auch die Syntax erklärt durch entgegen. Wir schauen, ob sich das Alter im Datensatz epa in den Geschlechtergruppen unterscheidet, also Alter erklärt durch Geschlecht.

# t-test Alter zwischen Geschlechtern im Datensatz "epa"
t.test(age ~ sex, data=epa)
## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = -4.7535, df = 595.03, p-value = 2.51e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -9.923268 -4.120822
## sample estimates:
## mean in group m mean in group w 
##        59.79646        66.81851

Das Ergebnis ist signifikant. Auch hier lässt sich wieder die Effektstärke Cohen’s D berechnen.

# Effektsärke 
lsr::cohensD(age ~ sex, data=epa)
## [1] 0.3837919

Der Effekt ist mit 0.3837919 eher klein.

Ein gutes Anleitungsvideo zum unabhängigen t-Test findet sich bei Youtube: https://www.youtube.com/watch?v=3Yz-XNRasGM

8.2.3.3.3 abhängiger t-Test

Bei gepaarten Werten muss ein abhängiger t-Test berechnet werden. Dies erfolgt über den Parameter paired, der auf TRUE gesetzt wird. Auch hier können wir wieder die Richtung (größer, kleiner, zweiseitig) über den Parameter alternative verändern.

# erzeuge Testwerte
t0 <- rnorm(100, mean=12, sd=1)
t1 <- rnorm(100, mean=11, sd=3)

# gepaarter t-Test
t.test(t0, t1, paired=TRUE)
## 
##  Paired t-test
## 
## data:  t0 and t1
## t = 2.5679, df = 99, p-value = 0.01173
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.189944 1.481379
## sample estimates:
## mean of the differences 
##               0.8356614
# gepaarter  t-test einseitig, kleiner
t.test(t0, t1, paired=TRUE, alternative = "less")
# gepaarter t-test einseitig, größer
t.test(t0, t1, paired=TRUE, alternative = "greater")
# gepaarter t-test zweiseitig
t.test(t0, t1, paired=TRUE, alternative = "two.sided")

Die Effektstärke Cohen’s D lässt sich entsprechend berechnen, indem der Parameter method auf paired gesetzt wird:

# Effektsärke 
lsr::cohensD(t0, t1, method="paired")
## [1] 0.2567894

Der Effekt ist mit 0.2567894 eher gering.

Ein gutes Anleitungsvideo zum verbundenen t-Test findet sich bei Youtube: https://www.youtube.com/watch?v=H15M_8gh1Ok

8.2.3.4 F-Test

Der F-Test wird mit der Funktion var.test() berechnet. So können wir prüfen, ob die Varianz des Alters im Datensatz epa in den Geschlechtergruppen unterschiedlich ist.

# lade Testdaten
load(url("https://www.produnis.de/R/epa.RData"))

# F-Test 
var.test(age ~ sex , data=epa)
## 
##  F test to compare two variances
## 
## data:  age by sex
## F = 0.98352, num df = 338, denom df = 280, p-value = 0.8816
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.7848842 1.2294811
## sample estimates:
## ratio of variances 
##           0.983521

Der Test ist nicht signifikant, es liegt kein Unterschied in den Varianzen vor.

8.2.3.5 Levene-Test

Ebenso kann mit dem Levene-Test auf Varianzunterschied geprüft werden. Hierfür rufen wir die Funktion leveneTest() aus dem Zusatzpaket car auf.

# lade Testdaten
load(url("https://www.produnis.de/R/epa.RData"))

# Levene-Test 
car::leveneTest(epa$age, epa$sex)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.0232 0.8789
##       618

Der Test ist nicht signifikant, es liegt kein Unterschied in den Varianzen vor. Ein gutes Anleitungsvideo zum verbundenen Levene-Test findet sich bei Youtube: https://www.youtube.com/watch?v=717UWGnyQGA

8.2.3.6 Mann-Whitney-U-Test

Dies ist die nicht-parametrische alternative zum t-Test für unabhängige Stichproben. Er wird mit der Funktion wilcox.test() durchgeführt, wobei der Parameter paired auf FALSE gesetzt wird (ist die Standardeinstellung).

# lade Testdaten
load(url("https://www.produnis.de/R/epa.RData"))

# Man-Whitney-U-Test
wilcox.test(epa$age ~ epa$sex)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  epa$age by epa$sex
## W = 35862, p-value = 1.151e-07
## alternative hypothesis: true location shift is not equal to 0

Das Ergebnis ist signifikant, das Alter im Datensatz epa unterscheidet sich zwischen den Geschlechtern.

Mit dem Parameter conf.int kann das Konfidenzintervall des Schätzwertes mit ausgegeben werden.

# Man-Whitney-U-Test mit Konfidenzintervall
wilcox.test(epa$age ~ epa$sex, conf.int=TRUE)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  epa$age by epa$sex
## W = 35862, p-value = 1.151e-07
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -9.999986 -4.999952
## sample estimates:
## difference in location 
##              -7.000029

Ein gutes Anleitungsvideo zum Mann-Whitney-U-Test finden Sei bei Youtube: https://www.youtube.com/watch?v=4lFyRcXoJB8

8.2.3.7 Wilcoxon-Test

Dies ist die nicht-parametrische alternative zum abhängigen t-Test. Er wird mit der Funktion wilcox.test() durchfegührt, wobei der Parameter paired auf TRUE gesetzt werden muss. Andernfalls rechnet R einen Man-Whitney-U-Test.

# erzeuge Testwerte
t0 <- sample(70:140, 100, replace=T)
t1 <- sample(74:170, 100, replace=T)

wilcox.test(t0, t1, paired=TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  t0 and t1
## V = 1436.5, p-value = 0.0002909
## alternative hypothesis: true location shift is not equal to 0

Mit dem Parameter conf.int kann das Konfidenzintervall des Schätzwertes mit ausgegeben werden.

# erzeuge Testwerte
t0 <- sample(70:140, 100, replace=T)
t1 <- sample(74:170, 100, replace=T)

wilcox.test(t0, t1, paired=TRUE, conf.int=TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  t0 and t1
## V = 1073, p-value = 6.001e-07
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -25.50005 -11.99994
## sample estimates:
## (pseudo)median 
##      -18.99996

Ein gutes Anleitungsvideo zum Wilcoxon-Test finden Sie bei Youtube: https://www.youtube.com/watch?v=4lFyRcXoJB8

8.2.3.8 Shapiro-Wilk-Normalverteilungstest

Zum Test auf Normalverteilung kann der Shapiro-Wilk-Test über die Funktion shapiro.test() verwendet werden. Die Nullhypothese lautet, dass eine Normalverteilung vorliegt. Testen wir, ob im Datensatz nw das Alter (age) normalverteilt ist.

# lade Testdaten
load(url("https://www.produnis.de/R/nw.RData"))

# Teste "age" auf Normalverteilung
shapiro.test(nw$age)
## 
##  Shapiro-Wilk normality test
## 
## data:  nw$age
## W = 0.96636, p-value = 4.677e-06

Das Ergebnis ist signifikant, die Daten sind nicht normalverteilt.

Erzeugen wir normalverteilte Werte zur Gegenprobe.

# Test mit 100 zufälligen und normalverteilten Werten
shapiro.test(rnorm(100))
## 
##  Shapiro-Wilk normality test
## 
## data:  rnorm(100)
## W = 0.99514, p-value = 0.9787

Das Ergebnis ist (wie erwartet) nicht signifikant, das heisst, unsere Zufallszahlen sind wirklich normalverteilt.

8.2.4 Fallzahlkalkulation

8.2.4.1 klinische Studien

Zur Fallzahlkalkulation bei klinischen Studien (z.B. RCTs) kann das Zusatzpaket pwr installiert werden. Mit der Funktion pwr.t.test() können alle Werte berechnet werden. Ihr müssen die entsprechenden Werte für das Signifikanzlevel, Power und Cohen’s d (Effektstärke) übergeben werden.

# erechne Fallzahl mit
# Cohens d = 0.1
# alpha = 0.05
# Power = 0.8
pwr::pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample", alternative="two.sided")
## 
##      Two-sample t test power calculation 
## 
##               n = 1570.733
##               d = 0.1
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Über die Parameter type und alternative lassen sich die Testbedingungen festlegen.

# Lade Paket
library(pwr)

# unabhängige Stichprobe, zweiseitig
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample", alternative="two.sided")

# Einstichproben-Test, kleiner
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="one.sample", alternative="less")

# gepaarte Stichprobe, größer
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="paired", alternative="greater")

Möchten Sie nicht immer die gesamte Testausgabe angezeigt bekommen, sondern nur den Wert für n, können Sie dies mit einem Dollarzeichen $ hinter dem Funktionsaufruf spezifizieren.

# zeige nur "n"
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample")$n
## [1] 1570.733
# runde auf nächste ganze Zahl
ceiling(pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample")$n)
## [1] 1571

Sie benötigen 1571 Probanden pro Untersuchungsgruppe.

8.2.4.2 Umfragen

Soll eine Befragung (ein Survey) durchgeführt werden, kann die benötigte Fallzahl über das Zusatzpaket samplingbook errechnet werden.

Geht es um numerische Werte, wird die Funktion sample.size.mean() verwendet. Ihr übergibt man Fehlerspanne e (ein Präzisionswert, der besagt, wie breit das Konfidenzintervall sein darf), das Konfidenzniveau, Standardabweichung der Grundgesamtheit sowie die Populationsgröße. Standardmäßig ist die Populationsgröße auf N=Inf (für unendlich) gesetzt.

# Berechne Fallzahl für einen Survey
# Konfidenzniveau auf 95%
# In der Grundgesamtheit hat der Wert eine
# Standardabweichung von S=10
# e = 0,5
samplingbook::sample.size.mean(e=0.5, S=10, level=0.95 )
## 
## sample.size.mean object: Sample size for mean estimate
## Without finite population correction: N=Inf, precision e=0.5 and standard deviation S=10
## 
## Sample size needed: 1537

Ist die Grundgesamtheit bekannt, kann dies mit N angegeben werden.

# Die Grundgesamtheit ist 1.500 Personen groß
samplingbook::sample.size.mean(e=0.5, S=10, level=0.95, N=1500 )
## 
## sample.size.mean object: Sample size for mean estimate
## With finite population correction: N=1500, precision e=0.5 and standard deviation S=10
## 
## Sample size needed: 760

Sollen Prozentwerte kategorialer Variablen erfragt werden (z.B. der Anteil weiblicher Pflegepersonen), wird die Funktion sample.size.prop() verwendet. Ihr werden ebenso e, N und level (Konfidenzniveau) übergeben. Statt der Standardabweichung wird über den Parameter P die prozentuale Verteilung in der Grundgesamtheit angegeben.

# Berechne Fallzahl für einen Survey
# unsere Werte düren maximal 5% abweichen
# Konfidenzniveau auf 95%
# In der Grundgesamtheit hat der Wert eine
# Verteilung von 65%
# Die Grundgesamtheit ist 1.500 Personen groß
samplingbook::sample.size.prop(e=0.05, P=0.65, N=1500)
## 
## sample.size.prop object: Sample size for proportion estimate
## With finite population correction: N=1500, precision e=0.05 and expected proportion P=0.65
## 
## Sample size needed: 284

8.3 Diagramme plotten

In diesem Kapitel werden die Plotfunktionen des klassischen R vorgestellt. Sie erzeugen immernoch brauchbare Ergebnisse, und man kann viele Einstellungen vornehmen. Wenn Sie allerdings ein etwas komplexeres Diagramm aufbauen möchten, lohnt es sich nicht, die Spezialeinstellungen der klassischen Plotfunktionen zu durchstöbern. Schwenken Sie in diesem Falle lieber gleich zum moderneren ggplot() um, welches in Kapitel @ref(Kapitel_ggplot) vorgestellt wird. So kommen Sie viel “einfacher” an Ihre Ziel.

Die universelle Funktion zum erstellen von Diagrammen heisst plot(). Mit ihr werden ein Großteil der Diagramme erstellt, und viele Funktionen rufen letztendlich auch nur plot() auf.

Da plot() so universal ist, nimmt sie viele Parameter entgegen. Die nachfolgende Tabelle zeigt diese Parameter in einer Übersicht. Viele der Parameter sind auch in weiteren Plotfunktionen implementiert. Es kann dennoch nie schaden, sich mit ?barplot in der R-Konsole die explizieten Parameter zu vergegenwärtigen.


Parameter der Funktion plot()
Parameter Argumente / Beschreibung Beispiel
type legt den Diagrammtyp fest type="p"
(mittlerweile veraltet, l = Liniendiagramm
funktioniert aber noch) p = Punktdiagramm
b = (both) Kombination aus l und p
s = Stufendiagramm
h = Histogramm-artig(!)
n = keine Darstellung
col Füllfarbe col="blue"
(colors() zeigt alle Farben) col="#112233" (RGB)
gray = schwarz bis gray99 = weiß col="gray77"
main Überschrift main="Mein Plot"
sub Untertitel main="Abbildung 12"
xlab Beschriftung der X-Achse xlab="Jahrgänge"
ylab Beschriftung der Y-Achse ylab="Häufigkeit"
xlim Start- und Endwert der X-Achse xlim=c(2,10)
ylim Start- und Endwert der Y-Achse ylim=c(0,100)
xaxt X-Achse komplett unbeschriftet xaxt="n"
yaxt Y-Achse komplett unbeschriftet yaxt="n"
pch Darstellungssymbol pch=4
(
cex Symbolgröße cex=2
lty Linientyp lty=1
(
lwd Strichstärke der Linien lwd=3
bty Boxtyp, setzt Rahmen ums Plot bty="o"
n = kein Rahmen
o = komplett
C = oben, unten und links
L = unten und links
u = unten, rechts und links
7 = oben und rechts
(die Buchstaben sehen so aus wie der Rahmen)
font Fontformat font=2
1 = normal
2 = fett
3 = kursiv
5 = griechische Symbole


8.3.1 Punktwolke

Die Punktwolke wird mit besagter Funktion plot() erstellt. Ihr werden zwei Vektoren mit den entsprechenden x und y Werten übergeben.

# erzeuge Testwerte
x <- sort(rnorm(100, mean=12, sd=2))
y <- sort(rnorm(100, mean=140, sd=15))

# Plotten
plot(x,y) 

Über die Parameter aus der Tabelle kann das Plot aufgehübscht werden.

plot(x,y, col="blue", 
          ylab="Körpergröße", xlab="Alter", 
          main="Ein Plot mit R base",
          sub="Dieses Plot wurde aus Zufallswerten erstellt",
          pch=13) 

Mit der Funktion abline() kann dem Plot eine Gerade hinzugefügt werden, z.B. eine Regressionsgerade. Die Funktion benötigt als Angaben den Schnitt durch die Y-Achse sowie den Steigungswert. Die Funktion lm() liefert für Regressionsgeraden genau diese Werte.

# Regressionswerte von Körpergröße ~ Alter
lm(y~x)
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      37.157        8.323

Das heisst, wir rufen die Funktion innerhalb von abline() auf:

plot(x,y, col="blue", 
          ylab="Körpergröße", xlab="Alter", 
          main="Ein Plot mit R base",
          sub="Dieses Plot wurde aus Zufallswerten erstellt",
          pch=13) 
abline(lm(y~x))

In rot sieht sie evtl. besser aus:

plot(x,y, col="blue", 
          ylab="Körpergröße", xlab="Alter", 
          main="Ein Plot mit R base",
          sub="Dieses Plot wurde aus Zufallswerten erstellt",
          pch=13) 
abline(lm(y~x), col="red")

Mit der Funktion lines() können Linien hinzugefügt werden. Dafür übergeben wir die entsprechenden Koordinaten des Start- und Endpunkts. So können wir z.B. den Mittelwert für Alter markieren.

# Mittelwer als x-Koordinate
mx <- mean(x)

# Berechne Schnittpunkt mit Regressionsgeraden
ymx <-  coef(lm(y~x))[1] + mx*coef(lm(y~x))[2]

plot(x,y, col="blue", 
          ylab="Körpergröße", xlab="Alter", 
          main="Ein Plot mit R base",
          sub="Dieses Plot wurde aus Zufallswerten erstellt",
          pch=13) 
abline(lm(y~x), col="red")

# Füge Mittelwertlinie hinzu
lines(x=c(mx, mx), y=c(0, ymx), lty=3)

Mit der Funktion polygon() können Polygone hinzugefügt werden.

# Mittelwer als x-Koordinate
mx <- mean(x)

# Berechne Schnittpunkt mit Regressionsgeraden
ymx <-  coef(lm(y~x))[1] + mx*coef(lm(y~x))[2]

plot(x,y, col="blue", 
          ylab="Körpergröße", xlab="Alter", 
          main="Ein Plot mit R base",
          sub="Dieses Plot wurde aus Zufallswerten erstellt",
          pch=13) 
abline(lm(y~x), col="red")

# Füge Mittelwertlinie hinzu
lines(x=c(mx, mx), y=c(0, ymx), lty=3)

# füge Polygon hinzu
polygon(x=c(x, max(x), min(x)), y=c(y, 0, 0), col="snow2")

So kann ein Diagramm der Standardnormalverteilung erstellt werden.

# erzeuge 10.000 Zufallswerte der Standardnormalverteilung
x <- sort(rnorm(10000))

# erzeuge die passenden Dichterwerte
y <- dnorm(x)

# beginne mit Plot
plot(x, y, ylab="Dichtefunktion", xlab="x", main="Die Standardnormalverteilung",
     sub="Bereich von 1 Standardabweichung") 

# füge gefülltes Polygon hinzu
polygon(x=c(x, max(x), min(x)), y=c(y, 0,0 ), col="snow2") 

# füge 1 Standardabweichung hinzu
lines(x=c(-1, -1), y=c(0, dnorm(-1)), lty=3)
lines(x=c(1, 1), y=c(0, dnorm(1)), lty=3)

# füge Polygon für 1. Standardbweichung hinzu
xsd1 <- x[x<1 & x>-1]
ysd1 <- dnorm(xsd1)
polygon(x=c(xsd1, max(xsd1), min(xsd1)), y=c(ysd1, 0, 0 ), col="green") 

Der Bereich von zwei Standardabweichungen kann wie folgt hinzugefügt werden. Dies “übermalt” jedoch die grüne Fläche.

plot(x, y, ylab="Dichtefunktion", xlab="x", main="Die Standardnormalverteilung",
     sub="Bereich von 2 Standardabweichungen") 

# füge 2. Standardabweichung hinzu
lines(x=c(-2, -2), y=c(0, dnorm(-2)), lty=6)
lines(x=c(2, 2), y=c(0, dnorm(2)), lty=6)

# füge Polygon für 2. Standardabweichung
# dies übermalt das 1. Polygon
xsd2 <- x[x<2 & x>-2]
ysd2 <- dnorm(xsd2)
polygon(x=c(xsd2, max(xsd2), min(xsd2)), y=c(ysd2, 0, 0 ), col="blue") 

Liegen mehrere Variablen in einem Datenframe vor, stellt plot() nach Kombinationen getrennt dar.

x <- data.frame(Lesetest = sort(rnorm(100, mean=20, sd=6)),
                Rechentest = rnorm(100, mean=15, sd=12),
                Rechtschreibung = sort(rnorm(100, mean=21, sd=4))
)

plot(x)

Je mehr Variablen vorliegen, desto komplexer wird der Plot.

x <- data.frame(Lesetest = sort(rnorm(100, mean=20, sd=6)),
                Rechentest = rnorm(100, mean=15, sd=12),
                Rechtschreibung = sort(rnorm(100, mean=21, sd=4)),
                Grammatikfehler = sort(rnorm(100, mean=21, sd=4), decreasing=T)
)

plot(x)

Mit Farben kann man es etwas aufhübschen.

x <- data.frame(Lesetest = sort(rnorm(100, mean=20, sd=6)),
                Rechentest = rnorm(100, mean=15, sd=12),
                Rechtschreibung = sort(rnorm(100, mean=21, sd=4)),
                Grammatikfehler = sort(rnorm(100, mean=21, sd=4), decreasing=T)
)

plot(x, col=rainbow(3))

8.3.2 verfügbare Farben

Die Funktion colors() zeigt die verfügbaren Farben an.

# zeige verfügbare Farben an
colors()
##   [1] "white"                "aliceblue"            "antiquewhite"        
##   [4] "antiquewhite1"        "antiquewhite2"        "antiquewhite3"       
##   [7] "antiquewhite4"        "aquamarine"           "aquamarine1"         
##  [10] "aquamarine2"          "aquamarine3"          "aquamarine4"         
##  [13] "azure"                "azure1"               "azure2"              
##  [16] "azure3"               "azure4"               "beige"               
##  [19] "bisque"               "bisque1"              "bisque2"             
##  [22] "bisque3"              "bisque4"              "black"               
##  [25] "blanchedalmond"       "blue"                 "blue1"               
##  [28] "blue2"                "blue3"                "blue4"               
##  [31] "blueviolet"           "brown"                "brown1"              
##  [34] "brown2"               "brown3"               "brown4"              
##  [37] "burlywood"            "burlywood1"           "burlywood2"          
##  [40] "burlywood3"           "burlywood4"           "cadetblue"           
##  [43] "cadetblue1"           "cadetblue2"           "cadetblue3"          
##  [46] "cadetblue4"           "chartreuse"           "chartreuse1"         
##  [49] "chartreuse2"          "chartreuse3"          "chartreuse4"         
##  [52] "chocolate"            "chocolate1"           "chocolate2"          
##  [55] "chocolate3"           "chocolate4"           "coral"               
##  [58] "coral1"               "coral2"               "coral3"              
##  [61] "coral4"               "cornflowerblue"       "cornsilk"            
##  [64] "cornsilk1"            "cornsilk2"            "cornsilk3"           
##  [67] "cornsilk4"            "cyan"                 "cyan1"               
##  [70] "cyan2"                "cyan3"                "cyan4"               
##  [73] "darkblue"             "darkcyan"             "darkgoldenrod"       
##  [76] "darkgoldenrod1"       "darkgoldenrod2"       "darkgoldenrod3"      
##  [79] "darkgoldenrod4"       "darkgray"             "darkgreen"           
##  [82] "darkgrey"             "darkkhaki"            "darkmagenta"         
##  [85] "darkolivegreen"       "darkolivegreen1"      "darkolivegreen2"     
##  [88] "darkolivegreen3"      "darkolivegreen4"      "darkorange"          
##  [91] "darkorange1"          "darkorange2"          "darkorange3"         
##  [94] "darkorange4"          "darkorchid"           "darkorchid1"         
##  [97] "darkorchid2"          "darkorchid3"          "darkorchid4"         
## [100] "darkred"              "darksalmon"           "darkseagreen"        
## [103] "darkseagreen1"        "darkseagreen2"        "darkseagreen3"       
## [106] "darkseagreen4"        "darkslateblue"        "darkslategray"       
## [109] "darkslategray1"       "darkslategray2"       "darkslategray3"      
## [112] "darkslategray4"       "darkslategrey"        "darkturquoise"       
## [115] "darkviolet"           "deeppink"             "deeppink1"           
## [118] "deeppink2"            "deeppink3"            "deeppink4"           
## [121] "deepskyblue"          "deepskyblue1"         "deepskyblue2"        
## [124] "deepskyblue3"         "deepskyblue4"         "dimgray"             
## [127] "dimgrey"              "dodgerblue"           "dodgerblue1"         
## [130] "dodgerblue2"          "dodgerblue3"          "dodgerblue4"         
## [133] "firebrick"            "firebrick1"           "firebrick2"          
## [136] "firebrick3"           "firebrick4"           "floralwhite"         
## [139] "forestgreen"          "gainsboro"            "ghostwhite"          
## [142] "gold"                 "gold1"                "gold2"               
## [145] "gold3"                "gold4"                "goldenrod"           
## [148] "goldenrod1"           "goldenrod2"           "goldenrod3"          
## [151] "goldenrod4"           "gray"                 "gray0"               
## [154] "gray1"                "gray2"                "gray3"               
## [157] "gray4"                "gray5"                "gray6"               
## [160] "gray7"                "gray8"                "gray9"               
## [163] "gray10"               "gray11"               "gray12"              
## [166] "gray13"               "gray14"               "gray15"              
## [169] "gray16"               "gray17"               "gray18"              
## [172] "gray19"               "gray20"               "gray21"              
## [175] "gray22"               "gray23"               "gray24"              
## [178] "gray25"               "gray26"               "gray27"              
## [181] "gray28"               "gray29"               "gray30"              
## [184] "gray31"               "gray32"               "gray33"              
## [187] "gray34"               "gray35"               "gray36"              
## [190] "gray37"               "gray38"               "gray39"              
## [193] "gray40"               "gray41"               "gray42"              
## [196] "gray43"               "gray44"               "gray45"              
## [199] "gray46"               "gray47"               "gray48"              
## [202] "gray49"               "gray50"               "gray51"              
## [205] "gray52"               "gray53"               "gray54"              
## [208] "gray55"               "gray56"               "gray57"              
## [211] "gray58"               "gray59"               "gray60"              
## [214] "gray61"               "gray62"               "gray63"              
## [217] "gray64"               "gray65"               "gray66"              
## [220] "gray67"               "gray68"               "gray69"              
## [223] "gray70"               "gray71"               "gray72"              
## [226] "gray73"               "gray74"               "gray75"              
## [229] "gray76"               "gray77"               "gray78"              
## [232] "gray79"               "gray80"               "gray81"              
## [235] "gray82"               "gray83"               "gray84"              
## [238] "gray85"               "gray86"               "gray87"              
## [241] "gray88"               "gray89"               "gray90"              
## [244] "gray91"               "gray92"               "gray93"              
## [247] "gray94"               "gray95"               "gray96"              
## [250] "gray97"               "gray98"               "gray99"              
## [253] "gray100"              "green"                "green1"              
## [256] "green2"               "green3"               "green4"              
## [259] "greenyellow"          "grey"                 "grey0"               
## [262] "grey1"                "grey2"                "grey3"               
## [265] "grey4"                "grey5"                "grey6"               
## [268] "grey7"                "grey8"                "grey9"               
## [271] "grey10"               "grey11"               "grey12"              
## [274] "grey13"               "grey14"               "grey15"              
## [277] "grey16"               "grey17"               "grey18"              
## [280] "grey19"               "grey20"               "grey21"              
## [283] "grey22"               "grey23"               "grey24"              
## [286] "grey25"               "grey26"               "grey27"              
## [289] "grey28"               "grey29"               "grey30"              
## [292] "grey31"               "grey32"               "grey33"              
## [295] "grey34"               "grey35"               "grey36"              
## [298] "grey37"               "grey38"               "grey39"              
## [301] "grey40"               "grey41"               "grey42"              
## [304] "grey43"               "grey44"               "grey45"              
## [307] "grey46"               "grey47"               "grey48"              
## [310] "grey49"               "grey50"               "grey51"              
## [313] "grey52"               "grey53"               "grey54"              
## [316] "grey55"               "grey56"               "grey57"              
## [319] "grey58"               "grey59"               "grey60"              
## [322] "grey61"               "grey62"               "grey63"              
## [325] "grey64"               "grey65"               "grey66"              
## [328] "grey67"               "grey68"               "grey69"              
## [331] "grey70"               "grey71"               "grey72"              
## [334] "grey73"               "grey74"               "grey75"              
## [337] "grey76"               "grey77"               "grey78"              
## [340] "grey79"               "grey80"               "grey81"              
## [343] "grey82"               "grey83"               "grey84"              
## [346] "grey85"               "grey86"               "grey87"              
## [349] "grey88"               "grey89"               "grey90"              
## [352] "grey91"               "grey92"               "grey93"              
## [355] "grey94"               "grey95"               "grey96"              
## [358] "grey97"               "grey98"               "grey99"              
## [361] "grey100"              "honeydew"             "honeydew1"           
## [364] "honeydew2"            "honeydew3"            "honeydew4"           
## [367] "hotpink"              "hotpink1"             "hotpink2"            
## [370] "hotpink3"             "hotpink4"             "indianred"           
## [373] "indianred1"           "indianred2"           "indianred3"          
## [376] "indianred4"           "ivory"                "ivory1"              
## [379] "ivory2"               "ivory3"               "ivory4"              
## [382] "khaki"                "khaki1"               "khaki2"              
## [385] "khaki3"               "khaki4"               "lavender"            
## [388] "lavenderblush"        "lavenderblush1"       "lavenderblush2"      
## [391] "lavenderblush3"       "lavenderblush4"       "lawngreen"           
## [394] "lemonchiffon"         "lemonchiffon1"        "lemonchiffon2"       
## [397] "lemonchiffon3"        "lemonchiffon4"        "lightblue"           
## [400] "lightblue1"           "lightblue2"           "lightblue3"          
## [403] "lightblue4"           "lightcoral"           "lightcyan"           
## [406] "lightcyan1"           "lightcyan2"           "lightcyan3"          
## [409] "lightcyan4"           "lightgoldenrod"       "lightgoldenrod1"     
## [412] "lightgoldenrod2"      "lightgoldenrod3"      "lightgoldenrod4"     
## [415] "lightgoldenrodyellow" "lightgray"            "lightgreen"          
## [418] "lightgrey"            "lightpink"            "lightpink1"          
## [421] "lightpink2"           "lightpink3"           "lightpink4"          
## [424] "lightsalmon"          "lightsalmon1"         "lightsalmon2"        
## [427] "lightsalmon3"         "lightsalmon4"         "lightseagreen"       
## [430] "lightskyblue"         "lightskyblue1"        "lightskyblue2"       
## [433] "lightskyblue3"        "lightskyblue4"        "lightslateblue"      
## [436] "lightslategray"       "lightslategrey"       "lightsteelblue"      
## [439] "lightsteelblue1"      "lightsteelblue2"      "lightsteelblue3"     
## [442] "lightsteelblue4"      "lightyellow"          "lightyellow1"        
## [445] "lightyellow2"         "lightyellow3"         "lightyellow4"        
## [448] "limegreen"            "linen"                "magenta"             
## [451] "magenta1"             "magenta2"             "magenta3"            
## [454] "magenta4"             "maroon"               "maroon1"             
## [457] "maroon2"              "maroon3"              "maroon4"             
## [460] "mediumaquamarine"     "mediumblue"           "mediumorchid"        
## [463] "mediumorchid1"        "mediumorchid2"        "mediumorchid3"       
## [466] "mediumorchid4"        "mediumpurple"         "mediumpurple1"       
## [469] "mediumpurple2"        "mediumpurple3"        "mediumpurple4"       
## [472] "mediumseagreen"       "mediumslateblue"      "mediumspringgreen"   
## [475] "mediumturquoise"      "mediumvioletred"      "midnightblue"        
## [478] "mintcream"            "mistyrose"            "mistyrose1"          
## [481] "mistyrose2"           "mistyrose3"           "mistyrose4"          
## [484] "moccasin"             "navajowhite"          "navajowhite1"        
## [487] "navajowhite2"         "navajowhite3"         "navajowhite4"        
## [490] "navy"                 "navyblue"             "oldlace"             
## [493] "olivedrab"            "olivedrab1"           "olivedrab2"          
## [496] "olivedrab3"           "olivedrab4"           "orange"              
## [499] "orange1"              "orange2"              "orange3"             
## [502] "orange4"              "orangered"            "orangered1"          
## [505] "orangered2"           "orangered3"           "orangered4"          
## [508] "orchid"               "orchid1"              "orchid2"             
## [511] "orchid3"              "orchid4"              "palegoldenrod"       
## [514] "palegreen"            "palegreen1"           "palegreen2"          
## [517] "palegreen3"           "palegreen4"           "paleturquoise"       
## [520] "paleturquoise1"       "paleturquoise2"       "paleturquoise3"      
## [523] "paleturquoise4"       "palevioletred"        "palevioletred1"      
## [526] "palevioletred2"       "palevioletred3"       "palevioletred4"      
## [529] "papayawhip"           "peachpuff"            "peachpuff1"          
## [532] "peachpuff2"           "peachpuff3"           "peachpuff4"          
## [535] "peru"                 "pink"                 "pink1"               
## [538] "pink2"                "pink3"                "pink4"               
## [541] "plum"                 "plum1"                "plum2"               
## [544] "plum3"                "plum4"                "powderblue"          
## [547] "purple"               "purple1"              "purple2"             
## [550] "purple3"              "purple4"              "red"                 
## [553] "red1"                 "red2"                 "red3"                
## [556] "red4"                 "rosybrown"            "rosybrown1"          
## [559] "rosybrown2"           "rosybrown3"           "rosybrown4"          
## [562] "royalblue"            "royalblue1"           "royalblue2"          
## [565] "royalblue3"           "royalblue4"           "saddlebrown"         
## [568] "salmon"               "salmon1"              "salmon2"             
## [571] "salmon3"              "salmon4"              "sandybrown"          
## [574] "seagreen"             "seagreen1"            "seagreen2"           
## [577] "seagreen3"            "seagreen4"            "seashell"            
## [580] "seashell1"            "seashell2"            "seashell3"           
## [583] "seashell4"            "sienna"               "sienna1"             
## [586] "sienna2"              "sienna3"              "sienna4"             
## [589] "skyblue"              "skyblue1"             "skyblue2"            
## [592] "skyblue3"             "skyblue4"             "slateblue"           
## [595] "slateblue1"           "slateblue2"           "slateblue3"          
## [598] "slateblue4"           "slategray"            "slategray1"          
## [601] "slategray2"           "slategray3"           "slategray4"          
## [604] "slategrey"            "snow"                 "snow1"               
## [607] "snow2"                "snow3"                "snow4"               
## [610] "springgreen"          "springgreen1"         "springgreen2"        
## [613] "springgreen3"         "springgreen4"         "steelblue"           
## [616] "steelblue1"           "steelblue2"           "steelblue3"          
## [619] "steelblue4"           "tan"                  "tan1"                
## [622] "tan2"                 "tan3"                 "tan4"                
## [625] "thistle"              "thistle1"             "thistle2"            
## [628] "thistle3"             "thistle4"             "tomato"              
## [631] "tomato1"              "tomato2"              "tomato3"             
## [634] "tomato4"              "turquoise"            "turquoise1"          
## [637] "turquoise2"           "turquoise3"           "turquoise4"          
## [640] "violet"               "violetred"            "violetred1"          
## [643] "violetred2"           "violetred3"           "violetred4"          
## [646] "wheat"                "wheat1"               "wheat2"              
## [649] "wheat3"               "wheat4"               "whitesmoke"          
## [652] "yellow"               "yellow1"              "yellow2"             
## [655] "yellow3"              "yellow4"              "yellowgreen"

Der Output ist 657 Farben lang.

# wieviele Farben gibt es?
length(colors())
## [1] 657

Zu erwähnen ist noch die Farbe gray, die in den Abstufungen gray0 (=schwarz) bis gray99 (= weiß) knapp 99 Grauschattierungen ermöglicht. Dies ist bei schwarz-weiß Publikationen sehr hilfreich.

Soll es möglichst bunt sein, erzeugt die Funktion rainbow() die gewünschte Anzahl an Farben.

8.3.3 Histogram

Ein Histogram wird mit der Funktion hist() erstellt. Plotten wir ein Histogramm der Variable Alter im Beispieldatensatz pf8.

# lade Datensatz "pf8"
load(url("https://www.produnis.de/R/pf8.RData"))

# erstelle Histogramm von Variable "Alter" aus Datensatz "pf8"
hist(pf8$Alter)

Zum Vergleich:

# zum Vergleich plot()
plot(pf8$Alter, type="h")

Wie Sie sehen ist das plot()-Resultat nur entfernt dem Histogram ähnlich.

Die Klassengröße des Histograms kann mit dem Parameter breaks() übergeben werden. Entweder schreiben wir die Anzahl der “Balken,” dann versucht R das “möglichst gut” auszuführen, oder wir übergeben einen Vektor mit den konkreten Breakpunkten.

# Histogramm mit "möglichst" 4 Balken
hist(pf8$Alter, breaks=4)

Wie Sie sehen stellt R dennoch 5 Balken dar, weil es so “besser” aussieht als mit nur 4 Balken. Die Anzahl der Balken ist also kein konkreter Befehl, sondern eher ein Serviervorschlag.

Über die Parameter kann der Plot nun aufgehübscht werden.

# Histogramm mit "möglichst" 4 Balken
# in schön
hist(pf8$Alter, breaks=4, 
     # Überschrift
     main="Mein Histogramm",
     # Unterschrift
     sub="Abbildung x.y: Ein Histogramm",
               # X-achse
               xlab="Alter in Jahren",
               # Y-Achse
               ylab="Häufigkeit",
               # blaue Balken
               col="blue")

8.3.4 Kreisdiagramm

Ein Kreisdiagramm wird mit der Funktion pie() erzeugt. Die Funktion benötigt einen Vektor mit den Anteilsgrößen.

# Kreisdiagramm 60 / 40
pie(c(60,40))

Soll z.B. ein Kreisdiagramm aus einem Faktor erzeugt werden, gelingt dies über die Integration der Funktion table().

# Kreisdiagramm von "Geschlecht" in "pf8"
pie(table(pf8$Geschlecht))

Mit den Parametern col, clockwise und labels können Sie weitere Anpassungen vornehmen.

pie(c(60,40), labels=c("ja", "nein"), col=c("azure", "coral"))

pie(table(pf8$Geschlecht), clockwise=TRUE, col=c("burlywood", "chartreuse","black"))

8.3.5 Säulendiagramm

Säulendiagramme werden mit der Funktion barplot() erstellt. Sie benötigt einen Vektor mit den jeweiligen Ausprägungen der Säulen sowie optional die Namen der Säulen über den Parameter names.arg:

# Säulendiagramm 
barplot(c(10,40,20), names.arg=c("Dritter", "Erster", "Zweiter"))

Liegen die Daten als Faktor vor, können sie mit Hilfe der Funktion table() übergeben werden.

# Säulendiagramm von "Standort" in "pf8"
barplot(table(pf8$Standort))

Liegen die Daten als wide table vor (z.B als Matrix mit bedeutsamen rownames), so wie unsere Matrix der Pflegeberufe, erstellt barplot() gruppierte Säulen. Über das Argument legend.text können die Reihennamen übergeben werden.

# Säulendiagramm der Matrix Pflegeberufe
barplot(Pflegeberufe, beside=T, col=rainbow(5), legend.text = rownames(Pflegeberufe))

Wird der Parameter beside weggelassen (entspricht FALSE), werden die Reihen in einer gemeinsamen Säule zusammengefasst.

# Säulendiagramm der Matrix Pflegeberufe
barplot(Pflegeberufe, col=rainbow(5), legend.text = rownames(Pflegeberufe))

Das entspricht vom Prinzip her der Abbildung 4.5.

8.3.6 Balkendiagramm

Balkendiagramme werden genau so erstellt wie Säulendiagramme, nämlich über die Funktion barplot() mit dem Parameter horiz=TRUE.

# Balkendiagramm der Matrix Pflegeberufe
# diesmal Schwarz-weiss
barplot(Pflegeberufe, 
            # schalte um auf "Balkendiagramm"
            horiz=TRUE, 
            # Graustufen   führende 0 kann weggelassen werden
            col=gray(c(.1, .2, .3, .4, .5)), 
            legend.text = rownames(Pflegeberufe))

# Balkendiagramm "category" aus "mma"
# diesmal Schwarz-weiss
barplot(table(mma$category), 
         horiz=TRUE, 
         col=rainbow(7), 
         # Legende aus Faktorenlevels
         legend.text = levels(mma$category))

8.3.7 Boxplot

Boxplots können mit der Funktion boxplot() erstellt werden.

# lade Datensatz
load(url("https://www.produnis.de/R/nw.RData"))

# Boxplot des Alters in Datensatz "nw"
boxplot(nw$age)

Die Funktion erkennt auch die Syntax “erklärt durch.”

# Boxplot des Alters in Datensatz "nw" erklärt durch "sex"
boxplot(age ~ sex, data=nw)

Über die Parameter frame, border, col und notch können weitere Anpassungen erfolgen.

boxplot(age ~ sex, data=nw, frame = FALSE, border="blue", horizontal = TRUE)

boxplot(age ~ sex, data=nw, frame = FALSE, col=c("green", "red"), notch = TRUE)

Weitere innere Gruppierungen können mit einem Sternchen * übergeben werden.

boxplot(age~sex*timework, data=nw, col=c("blue", "skyblue"), notch = TRUE)
## Warning in bxp(list(stats = structure(c(21, 35, 46, 52, 61, 22, 33, 42.5, : some
## notches went outside hinges ('box'): maybe set notch=FALSE

8.3.8 Diagramme speichern

Zum Speichern der Diagramme bietet die Formate jpeg, png, bmp und tiff an. Und genau so heissen auch die jeweiligen Speicherfunktionen. Dieser Befehl speichert das Diagramm im Format png:

# Speichere nachfolgendes als PNG
png("Testplot.png", width = 500, height = 500, units = "px", pointsize = 12)

# erzeuge Diagramm
pie(c(60,40))

# speichern abschließen 
dev.off()

Wie Sie sehen wird zunächst die Funktion png() aufgerufen. In ihr wird der Dateiname des Diagramms sowie dessen Dimensionen angegeben. Anschließend folgt der Aufruf der Plotbefehle, in diesem Falle pie(). Mit dem Befehl dev.off() (für device off ) wird das Speichern abgeschlossen. Die eigentlichen Diagrammbefehle müssen also zwischen der Speicherfunktion (png()) und dev.off() stehen.

Soll das Diagramm als jpeg gespeichert werden, lautet die Befehlskette entsprechend:

# Speichere nachfolgendes als JPEG
jpeg("Testplot2.jpg", width = 500, height = 500, units = "px", pointsize = 12)

# erzeuge Diagramm
barplot(c(10,40,20), names.arg=c("Dritter", "Erster", "Zweiter"))

# speichern abschließen
dev.off()

Genau so funktionieren auch die Funktionen bmp() und tiff(), aber diese Formate verwendet heute eigentlich niemand mehr.

8.4 ggplot

Das Zusatzpaket ggplot (für grammar of graphics plot) ist Hadley Wickhams8 R-Implementation der Grammar of Graphics von Leland Wilkinson9. Es ist eines der ersten Pakete des Tidyverse (welches damals noch nicht so hieß). Die Funktion zum plotten heisst “ggplot(),” das Installationspaket in R lautet jedoch ggplot2.

# ggplot2 installieren
install.packages("ggplot2", dependencies=T)

# ggplot2 aktivieren
library(ggplot2)

Da ggplot aus dem Tidyverse stammt (siehe Kapitel 6), ist es von Vorteil, wenn die zu verabreitenden Datensätze dem Prinzip “ein Fall pro Zeile” (tidy data bzw. long table) folgen. Das bedeutet, dass jede Beobachtung (auch Wiederholungen) in einer eigenen Zeile steht, und die jeweiligen Variablen durch die Spalten repräsentiert werden.

Die Idee einer Grammatik für Diagramme ist, dass man jedes Diagramm aus den selben Komponenten aufbauen kann:

  1. einem tidy Datensatz

  2. einem Koordinatensystem

  3. und Geomen (Punkte, Balken, Linien, Flächen)

In der Praxis ist ja auch so: wir haben einen Datensatz, und wir haben eine Vorstellung darüber, was wir gerne plotten möchten, z.B. ein Boxplot des Alters der Probanden nach Geschlecht gruppiert. Das Plott soll so aussehen, dass die Boxplots nebeneinander stehen.

Übersetzt bedeutet das, dass auf der X-Achse die Boxplots ausgerichtet sind, und dass die Y-Achse das Alter darstellt.

Diese Angaben werden an ggplot übergeben, wobei wir (ähnlich der Pipe) unsere Angaben mit einem + Zeichen aneinanderreihen können.

Die Zuweisung von Variablen auf das Koordinatensystem erfolgt über die Funktion aes() (für aesthetics). In userem Beispiel würden wir so beginnen:

# Lade Testdatensatz
load(url("https://www.produnis.de/R/nw.RData"))

# Grundlegende Zuweisung an gglot()
ggplot(nw, aes(x=sex, y=age))

Dies produziert noch kein Diagramm, da noch wichtige Angaben fehlen, aber wir haben die grundlegenden Zuweisungen vorgenommen. Aus dem Datensatz nw wurde die Variable sex der X-Achse zugewiesen, und Variable age wurde der Y-Achse zugewiesen. Jetzt muss ggplot() noch wissen, was überhaupt gezeichnet werden soll, also welches “Geom” zum Einsatz kommen soll.

In unserem Beispiel ist das geom_boxplot(), welches wir mit einem + Zeichen an ggplot() übergeben:

ggplot(nw, aes(x=sex, y=age)) +
   geom_boxplot()

Weitere Beispiele zu Boxplots finden Sie in Abschnitt 8.4.5.

Über das + Zeichen können beliebig viele Layer (Schichten) dem Plot hinzugefügt werden. Dabei werden die aesthetics übernommen, sie können aber in jedem Geom per aes() neu definiert werden.

ggplot biete hierfür viele nützliche Komponenten, unter anderem:

  • Geome (Flächen, Linien, Punkte, usw), beginnend mit geom_, z.B.:

    • gemo_area(), gemo_density(), gemo_smooth(), gemo_histogram(), gemo_point(), gemo_jitter(), gemo_line(), gemo_boxplot(), gemo_col(), gemo_text()
  • statistische Funktionen (Häufigkeiten, Mittelwerte, Dichte, uws), beginnend mit stat_, z.B.:

    • stat_count(), stat_density(), stat_function(), stat_summary()
  • Koordinatensystem, beginnend mit coor_, z.B.:

    • coord_flip(), coord_trans()
  • Aussehen (Theme), beginnend mit theme_, z.B.:

    • theme(), theme_bw(), theme_classic()
nützliche ggplot() Stylingbefehle
Eigenschaft Befehl
Plot Überschrift ggtitle("Überschrift", subtitle="Unterschrift")
theme(plot.title = element_text(size=8))
theme(plot.subtitle = element_text(size=8), face="bold"))
Plot Untertitel ggtitle(caption="Untertitel")
theme(plot.caption = element_text(size=8), face="italic")
Legende Titel theme(legend.title.x = element_text(size=8))
Legende Text theme(legend.text.x = element_text(size=8))
X-Achse Titel xlab("Titel der X-Achse")
theme(axis.title.x = element_text(size=8))
X-Achse Text theme(axis.text.x = element_text(size=8))
X-Achse Ticks theme(axis.ticks.x = element_text(size=8))
scale_x_continuous(breaks = c(0, 2, 6, 8, 15, 25))
X-Achse Länge xlim(-3,3)
Y-Achse Titel ylab("Titel der Y-Achse")
theme(axis.title.y = element_text(size=8))
Y-Achse Text theme(axis.text.y = element_text(size=8))
Y-Achse Ticks theme(axis.ticks.y = element_text(size=8))
scale_y_continuous(breaks = c(0, 2, 6, 8, 15, 25))
Y-Achse Länge ylim(-3,3)
Text annotate(geom="text", x=0, y=1, label="Huhu", color="black")
Plot aufteilen facet_grid(. ~ Gruppe, scales="free", space="free")
Plot drehen coord_flip()
Theme blank theme_void()
Theme klassisch theme_classic()
vertikale Linie geom_vline(xintercept=0, linetype="dotted")
horizontale Linie geom_hline(yintercept=0, linetype="dashed")

In RStudio ist ein “Cheatsheet” (siehe Abschnitt 7.1) zu ggplot verlinkt. Hier finden Sie (fast) alle verfügbaren Komponenten zusammengefasst. Klicken Sie in der Menüzeile von RStudio auf Help \(\rightarrow\) Cheatsheets \(\rightarrow\) Data Visualization with ggplot2.

8.4.1 Punktwolke

Schauen wir uns die Arbeitsjahre der Pflegenden im Nachtdienst an, indem wir eine Punktwolke erstellen. Auf der X-Achse sollen die Jahre als Pflegefachkraft abgebildet werden, auf der Y-Achse die Arbeitsjahre im Nachtdienst.

  ggplot(nw, aes(x=workyear, y=worknight)) +
    geom_point()

Als weitere aesthetic können wir festlegen, dass die Farben der Punkte nach Geschlecht unterschiedlich sein sollen. Im Geom geom_point() können wir zudem das Aussehen der Punkte über den Parameter shape ändern.

  ggplot(nw, aes(x=workyear, y=worknight, col=sex)) +
    geom_point(shape=8)

Beachten Sie, dass zwischen col (Rahmen) und fill (Füllung) unterschieden wird, wenn Sie Farben zuweisen.

Die Achsen möchten wir anders beschriften:

  ggplot(nw, aes(x=workyear, y=worknight, col=sex)) +
    geom_point(shape=8) +
    # Beschriftung der X-Achse
    xlab("Arbeitsjahre in der Pflege") +
    # Beschriftung der Y-Achse
    ylab("Arbeitsjahre als Nachtwache") +
    # Beschriftung der Legendenbox
    labs(color="Geschlecht")

Die Textgröße der Achsen könne angepasst werden:

  ggplot(nw, aes(x=workyear, y=worknight, col=sex)) +
    geom_point(shape=8) +
    # Beschriftung der X-Achse
    xlab("Arbeitsjahre in der Pflege") +
    # Beschriftung der Y-Achse
    ylab("Arbeitsjahre als Nachtwache") +
    # Beschriftung der Legendenbox
    labs(color="Geschlecht") +
    # Textgröße X-Beschriftung
    theme(axis.title.x = element_text(size=18)) + 
    # Textgröße der X-Ticks
    theme(axis.text.x = element_text(size=10)) + 
    # Textgröße Y-Beschriftung
    theme(axis.title.y = element_text(size=18)) + 
    # Textgröße der Y-Ticks
    theme(axis.text.y = element_text(size=10))

Nun können wir noch eine Regressionsgerade mit geom_smooth() hinzufügen:

  ggplot(nw, aes(x=workyear, y=worknight, col=sex)) +
    geom_point(shape=8) +
    # Beschriftung der X-Achse
    xlab("Arbeitsjahre in der Pflege") +
    # Beschriftung der Y-Achse
    ylab("Arbeitsjahre als Nachtwache") +
    # Beschriftung der Legendenbox
    labs(color="Geschlecht") +
    # Textgröße X-Beschriftung
    theme(axis.title.x = element_text(size=18)) + 
    # Textgröße der X-Ticks
    theme(axis.text.x = element_text(size=10)) + 
    # Textgröße Y-Beschriftung
    theme(axis.title.y = element_text(size=18)) + 
    # Textgröße der Y-Ticks
    theme(axis.text.y = element_text(size=10))  +
    # Regressionsgerade
    geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'

Das Plot könne wir zudem in verschiedene Bereiche aufteilen. Hierzu nehmen wir die Variable timework, die angibt, ob die Pflegefachpersonen in Voll- oder Teilzeit arbeiten. Die Variable übergeben wir der Komponente facet_grid().

  ggplot(nw, aes(x=workyear, y=worknight, col=sex)) +
    geom_point(shape=8) +
    # Beschriftung der X-Achse
    xlab("Arbeitsjahre in der Pflege") +
    # Beschriftung der Y-Achse
    ylab("Arbeitsjahre als Nachtwache") +
    # Beschriftung der Legendenbox
    labs(color="Geschlecht") +
    # Textgröße X-Beschriftung
    theme(axis.title.x = element_text(size=18)) + 
    # Textgröße der X-Ticks
    theme(axis.text.x = element_text(size=10)) + 
    # Textgröße Y-Beschriftung
    theme(axis.title.y = element_text(size=18)) + 
    # Textgröße der Y-Ticks
    theme(axis.text.y = element_text(size=10))  +
    # Regressionsgerade
    geom_smooth(method="lm") +
    # Plot aufteilen
    facet_grid(. ~ timework, scales = "free", space = "free")
## `geom_smooth()` using formula 'y ~ x'

Das Plot selbst kann jederzeit als Objekt gespeichert werden. So können Sie anschließend darauf zugreifen. Das ist sehr hilfreich, wenn man am Aussehen herumfeilt.

# Speichere als Objekt "p"
p <- ggplot(nw, aes(x=workyear, y=worknight, col=sex)) +
       # Beschriftung der X-Achse
       xlab("Arbeitsjahre in der Pflege") +
       # Beschriftung der Y-Achse
       ylab("Arbeitsjahre als Nachtwache") +
       # Beschriftung der Legendenbox
       labs(color="Geschlecht") +
       # Textgröße X-Beschriftung
       theme(axis.title.x = element_text(size=18)) + 
       # Textgröße der X-Ticks
       theme(axis.text.x = element_text(size=10)) + 
       # Textgröße Y-Beschriftung
       theme(axis.title.y = element_text(size=18)) + 
       # Textgröße der Y-Ticks
       theme(axis.text.y = element_text(size=10))  +
       geom_smooth(method="lm") +
       # Plot aufteilen
       facet_grid(. ~ timework, scales = "free", space = "free")

# was wurde gespeichert?
p
## `geom_smooth()` using formula 'y ~ x'

Der gespeicherte Plot enthält bereits die Regressionsgeraden.

Von hier aus können Sie leichter “herumspielen”:

# spiele mit "p" herum
p + geom_point(shape=2)  
## `geom_smooth()` using formula 'y ~ x'

# spiele mit "p" herum
p + geom_point(shape=10)  
## `geom_smooth()` using formula 'y ~ x'

Hier eine Auflistung, welche shape-Formen es gibt:

Abbildung 8.1: shape nach Zahlen

Abbildung 8.2: shape nach Formen

Natürlich funktioniert ggplot() auch in Kombination mit der Pipe.

nw %>% 
  drop_na() %>% 
    ggplot(aes(x=bed)) + 
    geom_bar(fill=rainbow(6)) + 
    coord_flip()

8.4.2 Diagramme speichern

Zum Speichern der Plots in eine Datei wird die Funktion ggsave() verwendet. Sie speichert das zuletzt erzeugte Plot.

# Plot speichern in Datei "MeinPlot.png"
ggsave("MeinPlot.png", units="mm", width=400, height=200, dpi=300, pointsize=7) 

Haben Sie ein Plot als Objekt zugewiesen, und möchten dieses nun in einer Datei speichern, übergeben Sie den Objektnamen mit dem Parameter plot.

# Plot liegt in Objekt "p"
# speichern in Datei "MeinPlot.png"
ggsave("MeinPlot.png", plot=p, units="mm", width=400, height=200, dpi=300, pointsize=7) 

8.4.3 Histrogramm

Histogramme können mit dem geom_histogram() hinzugefügt werden.

# Histogramm des Alters aus dem Datensatz "epa"
# lade Datensatz
load(url("https://www.produnis.de/R/epa.RData"))

epa %>% 
  ggplot(aes(x=age)) +
  geom_histogram(aes(y=..density..), col="white", fill="plum4") 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Mit der Funktion stat_density() können wir die Dichtefunktion des Alters hinzufügen:

epa %>% 
  ggplot(aes(x=age)) +
  geom_histogram(aes(y=..density..), col="white", fill="plum4") + 
  stat_density(geom="line", colour="blue", linetype = "dotted")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Die Dichteverteilung kann auch als Fläche angezeigt werden. Über den Parameter alpha kann eingestellt werden, wie durchsichtig die Fläche sein soll.

epa %>% 
  ggplot(aes(x=age)) +
  geom_histogram(aes(y=..density..), col="white", fill="plum4") + 
  stat_density(geom="area", colour="blue", linetype = "dotted", alpha=0.5)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Die Dichteverteilung der Normalverteilung kann mit der Funktion stat_funtion() ebenso darübergelegt werden.

epa %>% 
  ggplot(aes(x=age)) +
  geom_histogram(aes(y=..density..), col="white", fill="plum4") + 
  stat_density(geom="area", colour="blue", linetype = "dotted", alpha=0.5) +
  stat_function(fun=dnorm, args=(c(mean=62,sd=18)), colour="red", linetype = "dashed")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

8.4.4 Säulen- und Balkendiagramm

Säulendiagramme werden mit geom_bar() erzeugt. Die Werte der Variable müssen als Faktor vorliegen.

# Lade Testdatensatz
load(url("https://www.produnis.de/R/nw.RData"))

p <-nw %>% 
      drop_na() %>% 
         ggplot(aes(x=bed))
p + geom_bar()

Ein Balkendiagramm wird erzeugt, indem die Kompontent coord_flip() hinzugefügt wird.

p + geom_bar() + coord_flip()

Sie können gruppierte Plots über die Aesthetic fill erzeugen.

  nw %>% 
    drop_na() %>% 
    ggplot(aes(x=bed, fill=timework)) +
    geom_bar()

Über den Parameter position können Sie die Balken nebeneinander gruppieren…

  nw %>% 
    drop_na() %>% 
    ggplot(aes(x=bed, fill=timework)) +
    geom_bar(position="dodge")

… oder als Prozentwerte plotten.

  nw %>% 
    drop_na() %>% 
    ggplot(aes(x=bed, fill=timework)) +
    geom_bar(position="fill")

Möchten Sie ein Barplot aus bereits berechneten Häufigkeiten erstellen, müssen Sie dies dem Geom über den Parameter stat=identity mitteilen. Die Daten müssen als Datenframe übergeben werden.

 # erzeuge Testdaten
 kategorie     <- c("Erdbeere", "Schokolade", "Vanille")
 haeufigkeiten <- c(10,20,30)
 df <- data.frame(kategorie, haeufigkeiten)

 ggplot(df, aes(x=kategorie, y=haeufigkeiten, fill=kategorie))+
    geom_bar(stat="identity", col="black")

8.4.5 Boxlots

Boxplots werden über das geom_boxplot() erstellt. Aus dem Nachtwachendatensatz erstellen wir ein Boxplot über die “Anzahl der Dienste am Stück” (nights).

load(url("https://www.produnis.de/R/nw.RData"))
ggplot(nw, aes(y=nights))+
    geom_boxplot()

Das Aussehen kann beliebig angepasst werden.

ggplot(nw, aes(y=nights))+
    # Die Ausreisser stylen
    geom_boxplot(colour="darkblue", fill="skyblue", outlier.colour="blue", outlier.shape=8, outlier.size=3) +
    # zeichne Whisker-Linien
    stat_boxplot(geom ='errorbar') +
    # beschrifte Y-Achse
    ylab("Dienste hintereinander") + 
    # Textgröße Y-Achsenskala
    theme(axis.title.y = element_text(size=18)) + 
    # individuelle Ticks Y-Achse
    scale_y_continuous(breaks = c(0, 2, 4, 6, 8, 10 , 15, 25)) +
    # beschrifte X-Achse
    xlab(NULL) +
    # X-Achsenbeschriftung löschen
    theme(axis.ticks.x = element_blank(), axis.text.x = element_blank())

So können auch mehere Gruppen abgebildet werden, zum Beispiel nach Geschlecht gruppiert:

ggplot(nw, aes(x=sex, y=nights, fill=sex)) +
    # Die Ausreisser stylen
    geom_boxplot(colour="darkblue", outlier.colour="blue", outlier.shape=8, outlier.size=3) +
    # zeichne Whisker-Linien
    stat_boxplot(geom ='errorbar') +
    # beschrifte Y-Achse
    ylab("Dienste hintereinander") + 
    # Textgröße Y-Achsenskala
    theme(axis.title.y = element_text(size=18)) + 
    # individuelle Ticks Y-Achse
    scale_y_continuous(breaks = c(0, 2, 4, 6, 8, 10 , 15, 25)) +
    # beschrifte X-Achse
    xlab(NULL) +
    # X-Achsenbeschriftung löschen
    theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) +
    # Zeichne den Mittelwert als roten Punkt
    stat_summary(fun=mean, colour="darkred", geom="point", shape=18, size=5, show.legend = F)

Hier ein etwas komplexeres Beispiel aus dem Nachtwachen-Datensatz. Zuerst selektieren wir die Variablen zu Pflegestufen, Demenz und freiheitsentziehenden Maßnahmen (FEM) und bringen sie ins long table Format. Dann führen wir eine neue gruppierende Variable “Gruppe” ein. Diese dient als Trenner zwischen den Boxplotgruppen.

# Lade Nachtwachendatensatz
load(url("https://www.produnis.de/R/nw.RData"))
nwbp <- nw %>% 
  # suche die richtigen Spalten heraus
  select(Stufe0, Stufe1, Stufe2, Stufe3, Demenz, Bettgitter, Bettgurt, Schlafmittel) %>% 
  # erzeuge daraus eine "long table"
  pivot_longer(cols=c(Stufe0, Stufe1, Stufe2, Stufe3, Demenz, Bettgitter, Bettgurt, Schlafmittel), names_to="Variable", values_to="Anzahl") %>% 
  # Füge manuell die "Trenner"-Variable hinzu
  bind_cols(., Gruppe=factor(rep(c("Pflegestufe", "Pflegestufe", "Pflegestufe", "Pflegestufe", "Demenz", "FEM", "FEM", "FEM"), 276)))

ggplot(nwbp, aes(Variable, Anzahl))+
  geom_boxplot(outlier.size = 3,outlier.shape=1) +
  stat_boxplot(geom ='errorbar') +
  ylab(NULL) + xlab(NULL)+
  theme(axis.text.x = element_text(size=12))+
  theme(axis.text.y = element_text(size=11))+
  scale_y_continuous(limits=c(0, 70))+
  scale_y_continuous(breaks = c(0,10,20,30,40,50,60,70,100))+
  facet_grid(. ~ Gruppe, scales = "free", space = "free") # "Gruppe" ist der Spaltenname
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Warning: Removed 321 rows containing non-finite values (stat_boxplot).

## Warning: Removed 321 rows containing non-finite values (stat_boxplot).

8.4.6 Likert-Plots

Likkert-Plots können über das Zusatzpaket likert geplottet werden. Die Daten müssen in einem Datenframe vorliegen und an die Funktion likert() übergeben werden.

# Lade Nachtwachendatensatz
load(url("https://www.produnis.de/R/nw.RData"))

# Erzeuge Datenframe "Pause"
pause <- data.frame(nw$pause)

# übergebe Datenframe an "likert()"
pause <- likert::likert(pause)

# plotte
plot(pause)

Sollen mehrere Variablen gemeinsam ausgewertet werden, erzeugt man einfach ein Datenframe, wobei wie gewohnt jede Spalte eine Variable darstellt.

# Erzeuge neues Datenframe "n"
n <- data.frame(nw[, 29:34])

# übergebe an "likert()"
n <- likert::likert(n)

# plot
plot(n)

8.4.7 Kreisdiagramm

In ggplot() gibt es kein Geom, mit dem Kreisdiagramme erstellt werden können. Der Trick besteht darin, einfach ein geom_bar() auf einem circulären Koordinationsystem zu plotten. Dies erreichen wir, indem coord_polar() hinzugefügt wird.

Wir übergeben ggplot() ein Datenframe, mit einer Variable für die Gruppennamen, und einer Variable mit den Anteilswerten.

Dass die Anteile bereits berechnet sind teilen wir über den Parameter stat=“identity” mit.

# Erzeuge Datenframe
df <- data.frame(Gruppe=c("A", "B", "C"), Anteil=c(60, 25, 15))

# pfusche Pie-Chart zurecht
ggplot(df, aes(x="", y=Anteil, fill=Gruppe)) +
    geom_bar(stat="identity", width=1) + 
    coord_polar("y", start=0) +
    # entferne Achsen und Ticks
    theme_void()


  1. große Schlarmann, J; Galatsch, M (2014): Regressionsmodelle für ordinale Zielvariablen, GMS Medizinische Informatik, Biometrie und Epidemiologie, 10(1):2-10, doi: 10.3205/mibe000154, https://www.egms.de/static/en/journals/mibe/2014-10/mibe000154.shtml↩︎

  2. Cohen, J (1992): A Power Primer, Psychological Bulletin, 112(7):155-159↩︎

  3. Wickham, H (2010): A layered grammar of graphics, Journal of Computational and Graphical Statistics, 19(1):3-28, doi: 10.1198/jcgs.2009.07098, https://vita.had.co.nz/papers/layered-grammar.html↩︎

  4. Wilkinson, L (2005): The Grammar of Graphics, Springer, ISBN 978-0-387-28695-2↩︎