Kapitel 8 Statistik mit R und RStudio
Ebenso wie die Statistik, ist dieses Kapitel in folgende Bereiche unterteilt:
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 1Fü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.125Kumulierte 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 8Und 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.000Eine 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.58.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.00Die Funktion describe() aus dem psych-Zusatzpaket liefert noch detailiertere Übersichten:
psych::describe(x)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 8 8.75 10.43 6.5 8.75 1.48 1 34 33 1.69 1.35 3.69Das Funktioniert auch bei kompletten Datensätzen.
# Zusammenfassung von pf8
psych::describe(pf8)## vars n mean sd median trimmed mad min max range skew
## Standort* 1 731 2.85 1.61 2 2.81 1.48 1 5 4 0.28
## Alter 2 730 37.00 18.00 28 34.98 12.60 15 86 71 0.74
## Geschlecht* 3 725 1.61 0.49 2 1.63 0.00 1 3 2 -0.40
## Größe 4 730 174.04 9.62 173 173.84 10.38 151 206 55 0.26
## Gewicht 5 720 75.53 16.88 73 73.93 16.31 45 175 130 1.12
## Bildung* 6 506 4.88 1.77 5 4.98 2.97 1 7 6 -0.39
## Beruf* 7 731 40.42 37.04 34 38.58 48.93 1 104 103 0.26
## Familienstand* 8 729 2.11 1.12 2 1.99 1.48 1 6 5 0.83
## Kinder 9 727 0.38 0.79 0 0.18 0.00 0 5 5 2.19
## Wohnort* 10 729 1.41 0.49 1 1.39 0.00 1 2 1 0.36
## Rauchen* 11 729 1.21 0.40 1 1.13 0.00 1 2 1 1.45
## SportHäufig 12 677 2.56 1.88 2 2.41 1.48 0 21 21 1.91
## SportMinuten 13 667 62.42 40.14 60 60.47 44.48 0 400 400 1.53
## SportWie* 14 656 1.64 0.69 2 1.55 1.48 1 3 2 0.62
## SportWarum* 15 454 3.75 1.64 4 3.49 1.48 1 8 7 1.06
## LebenZufrieden 16 724 7.91 1.59 8 8.07 1.48 1 11 10 -1.10
## kurtosis se
## Standort* -1.53 0.06
## Alter -0.70 0.67
## Geschlecht* -1.75 0.02
## Größe -0.26 0.36
## Gewicht 2.40 0.63
## Bildung* -1.26 0.08
## Beruf* -1.55 1.37
## Familienstand* 0.48 0.04
## Kinder 4.45 0.03
## Wohnort* -1.87 0.02
## Rauchen* 0.11 0.01
## SportHäufig 12.75 0.07
## SportMinuten 8.64 1.55
## SportWie* -0.77 0.03
## SportWarum* 0.79 0.08
## LebenZufrieden 1.86 0.06Die Lagewerte lassen sich auch einzeln bestimmen.
Mit der Funktion min() wird das Minimum der Verteilung ausgegeben.
# Minimum bestimmen
min(x)## [1] 1Dementsprechend liefert max() das Maximum.
# Maximum bestimmen
max(x)## [1] 34Die Spannweite wird mit range() abgefragt.
# Spannweite bestimmen
range(x)## [1] 1 34Das arithmetische Mittel wird mit der Funktion mean() bestimmt
# x-quer bestimmen
mean(x)## [1] 8.75Sind fehlende Werte in der Reihe enthalten, gibt R nur NA zurück.
# mean für "Gewicht" im Datensatz "pf8"
mean(pf8$Gewicht)## [1] NABeheben 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.5325Der Median wird mit der Funktion median() ermittelt.
# Median bestimmen
median(x)## [1] 6.5Auch hier müssen enthaltene NAs ausgeblendet werden:
# Median bestimmen, OHNE NAs
median(pf8$Gewicht, na.rm=TRUE)## [1] 73In 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
## 3Es 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] 7Die 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.0Es lassen sich beliebige Quantile berechnen.
# Quartile bestimmen, ohne NAs
quantile(x, probs = c(.333, .666))## 33.3% 66.6%
## 5.331 7.000Und auch die Art der Berechnung lässt sich über den Parameter type festlegen. Beispielsweise verwendet type=6 die Methode von SPSS.
# Quartile bestimmen, ohne NAs
# so wie SPSS
quantile(x, probs = c(.333, .666), type=6)## 33.3% 66.6%
## 4.994 7.0008.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.43004Die Varianz der Stichprobe wird mit var() berechnet.
# Varianz von x bestimmen
var(x, na.rm=T)## [1] 108.7857Der Interquartilsabstand (oder einfach nur Quartilsabstand) wird mit der Funktion IQR() bestimmt.
# Quartilsabstand bestimmen
IQR(x, na.rm=T)## [1] 2.5Ebenso wie bei quantile() lässt sich über den Parameter type die Berechnungsart ändern.
Mit type=6 nutzt R die Methode wie in SPSS:
# Quartilsabstand bestimmen
# Methode wie in SPSS
IQR(x, na.rm=T, type=6)## [1] 3.58.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 1Ebenso 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 18.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.3758728.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.9884616Mit 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.9884616Sollten 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] NAMit 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.1067484Fü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.9969651Auch 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.05920716Auch 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.9969651Die 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.99696518.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.055Hierbei 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.8291Es 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-16Auß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-16Ebenfalls 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: 5Die Odds Ratio (OR) für \(\beta\)Temp= -0,2322 errechnet sich mit der exp()-Funktion:
exp(-0.2322)## [1] 0.7927875Dies 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.9984265Das 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 Galatsch7.
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 maessigDie 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.9129275Mit 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-11Die 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.5493pom.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.2928789Das 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.1065Der 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.106526Ist 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.2464parallel=T~Konflikt-1bedeutet, dass equal slopes nur fürKonfliktangenommen wird.parallel=F~Zufriedenhbedeutet, dass equal slopes nur fürZufriedenhnicht 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-03Mit 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.769318.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 estimatesDies 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 estimatesHat 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.2252Der 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.4627Mit 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-12get.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.571979get.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.28470388.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.968434Beachten 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.0316Die untere Grenze des 95%-Konfidenzintervalls liegt bei 396.0316 Euro.
# obere Grenze des 95%KI
xquer + fehler95## [1] 403.9684Die 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.936868Das 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.7472Die untere Grenze des 99%-Konfidenzintervalls liegt bei 394.7472 Euro.
# obere Grenze des 99%KI
xquer + fehler99## [1] 405.2528Die 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.50562Das 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.4125813Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.4125813 mg/dl.
# obere Grenze des 95%KI
d + fehler95## [1] 0.5474187Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.5474187 mg/dl.
# Länge des Intervalls
2*fehler95## [1] 0.1348373Das 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.3883663Die untere Grenze des 99%-Konfidenzintervalls liegt bei 0.3883663 mg/dl.
# obere Grenze des 99%KI
d + fehler99## [1] 0.5716337Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.5716337 mg/dl.
# Länge des Intervalls
2*fehler99## [1] 0.1832675Das 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.07632963Beachten 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.2736704Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.2736704 = 27,36704%.
# obere Grenze des KI
k + fehler95## [1] 0.4263296Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.4263296 = 42,63296%.
# Länge des KI
2*fehler95## [1] 0.1526593Das 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.2496859Die untere Grenze des 99%-Konfidenzintervalls liegt bei 0.2496859 = 24,96859%.
# obere Grenze des KI
k + fehler99## [1] 0.4503141Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.4503141 = 45,03141%.
# Länge des KI
2*fehler99## [1] 0.2006283Das 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 dˆ 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-p1Beachten 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.03446183Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.03446183 = 3,446183%.
# obere Grenze des 95%KI
d + fehler95## [1] 0.2655382Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.2655382 = 26,55382%.
# Länge des KI
2*fehler95## [1] 0.2310763Die 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.001842898Die untere Grenze des 99%-Konfidenzintervalls liegt bei -0.0018429 = -0,18429%.
# obere Grenze des 99%KI
d + fehler99## [1] 0.3018429Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.3018429 = 30,18429%.
# Länge des KI
2*fehler99## [1] 0.3036858Die 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 missingnessDas 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: holmDie 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: bonferroniDie 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.4765Der 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.22826Der 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.4109627Dieser Wert kann nun nach Cohen’s Angaben8 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.97903Das 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.81851Das 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.3837919Der 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 = 4.4894, df = 99, p-value = 1.933e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.730891 1.888698
## sample estimates:
## mean of the differences
## 1.309794# 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.4489378Der Effekt ist mit 0.4489378 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.983521Der 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
## 618Der 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 0Das 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.000029Ein 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 = 1282, p-value = 1.934e-05
## alternative hypothesis: true location shift is not equal to 0Mit 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 = 1348.5, p-value = 8.483e-05
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
## -22.999958 -8.499945
## sample estimates:
## (pseudo)median
## -15.99997Ein 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-06Das 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.98387, p-value = 0.2623Das 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] 1571Sie 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: 1537Ist 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: 760Sollen 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: 2848.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 | 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
## 57.593 6.808Das 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] 657Zu 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.
Eine Übersicht der Farben in “echten Farbtönen” findet sich bei Björn Walther: https://www.bjoernwalther.com/farben-in-r-der-col-befehl/
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"))
Und auch lustige Dinge kann man plotten9:
par(mar = c(0, 1, 0, 1))
pie(
c(280, 60, 20),
c('Sky', 'Sunny side of pyramid', 'Shady side of pyramid'),
col = c('#0292D8', '#F7EA39', '#C4B632'),
init.angle = -50, border = NA
)
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 Wickhams10
R-Implementation der Grammar of Graphics von Leland Wilkinson11. 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:
einem tidy Datensatz
einem Koordinatensystem
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()
…
| 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()
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↩︎
Cohen, J (1992): A Power Primer, Psychological Bulletin, 112(7):155-159↩︎
https://www.shirin-glander.de/2017/09/moving-my-blog-to-blogdown/↩︎
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↩︎
Wilkinson, L (2005): The Grammar of Graphics, Springer, ISBN 978-0-387-28695-2↩︎