31  Schließende Statistik

31.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.

31.1.1 Normalverteilung

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

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

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

  • \(n\) = Stichprobenumfang

  • \(s\) = Standardabweichung

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

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

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

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

# Fehlerwert ausgeben
fehler95
## [1] 3.968434

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

31.1.1.1 Unterschied von Mittelwerten

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

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

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

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

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

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

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

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

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

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

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

  • \(n_{A}\) = 13

  • \(n_{B}\) = 10

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

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

# Differenz der Mittelwerte
d   <- xb - xa

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

31.1.2 Binomialverteilung

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

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

  • \(n\) = Umfang der Stichprobe

  • \(k\) = geschätzter Anteil

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

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

# Daten erzeugen
n <- 150
k <- 0.35

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

# anzeigen
fehler95
## [1] 0.07632963

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

31.1.2.1 Unterschied von Anteilswerten

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

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

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

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

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

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

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

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

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

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

# Differenz der Anteile
d <- p2-p1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

31.2 Signifikanztests

31.2.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-Quadrat-Approximation kann inkorrekt sein
## 
##  Pearson's Chi-squared test
## 
## data:  xtabelle
## X-squared = 43.886, df = 44, p-value = 0.4765

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

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

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

31.2.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.

31.2.3 t-Test

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

31.2.3.1 Einstichproben t-Test

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

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

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

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

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

# aktiviere Zusatzpaket "lsr"
library(lsr)

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

Dieser Wert kann nun nach Cohen’s Angaben (1992) 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

31.2.4 unabhängiger t-Test

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

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

Das Ergebnis ist signifikant, die Gruppen unterscheiden sich also.

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

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

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

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

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

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

Der Effekt ist mit 0.3837919 eher klein.

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

31.2.4.1 abhängiger t-Test

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

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

# gepaarter t-Test
t.test(t0, t1, paired=TRUE)
## 
##  Paired t-test
## 
## data:  t0 and t1
## t = 2.2375, df = 99, p-value = 0.0275
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.08220316 1.37023444
## sample estimates:
## mean difference 
##       0.7262188
# 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.2237486

Der Effekt ist mit 0.2237486 eher gering.

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

31.2.5 F-Test

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

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

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

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

31.2.6 Levene-Test

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

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

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

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

31.2.7 Mann-Whitney-U-Test

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

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

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

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

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

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

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

31.2.8 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 = 1474.5, p-value = 0.0003056
## alternative hypothesis: true location shift is not equal to 0

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

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

wilcox.test(t0, t1, paired=TRUE, conf.int=TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  t0 and t1
## V = 809.5, p-value = 6.183e-09
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -33.49995 -18.49997
## sample estimates:
## (pseudo)median 
##      -26.49996

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

31.2.9 Shapiro-Wilk-Normalverteilungstest

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

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

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

Das Ergebnis ist signifikant, die Daten sind nicht normalverteilt.

Erzeugen wir normalverteilte Werte zur Gegenprobe.

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

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

31.3 Varianzanalysen

31.3.1 ANOVA mit Messwiederholung

Eine ANOVA mit Messwiederholung kann wie folgt in R durchgeführt werden. Wir nutzen den Testdatensatz Messwiederholung.rda.

load(url("https://www.produnis.de/R/Messwiederholung.rda"))

Der Datensatz besteht aus 200 Probanden, bei denen zu drei Zeitpunkten ein Wert erhoben wurde.

head(Messwiederholung)
##      Name        t1       t2       t3
## 1  Brycen 11.728110 8.190290 6.202859
## 2   Waail  9.781605 8.646142 6.311455
## 3 Abigail 11.937967 9.386623 4.577478
## 4   Holly  4.597950 8.204762 7.250849
## 5  Sheila  8.277670 8.322449 6.952457
## 6 Rebekah  7.891258 9.104856 7.191639

Zunächst müssen wir die Daten ins long table-Format überführen.

# Überführe die Daten ins "long table"-Format
df <- Messwiederholung %>%
  pivot_longer(cols=t1:t3,
               names_to = "Zeitpunkt",
               values_to = "Wert") %>%
  mutate(Zeitpunkt = factor(Zeitpunkt),
         Name = factor(Name))

# anschauen
head(df)
## # A tibble: 6 × 3
##   Name   Zeitpunkt  Wert
##   <fct>  <fct>     <dbl>
## 1 Brycen t1        11.7 
## 2 Brycen t2         8.19
## 3 Brycen t3         6.20
## 4 Waail  t1         9.78
## 5 Waail  t2         8.65
## 6 Waail  t3         6.31

Zunächst verschaffen wir uns einen deskriptiven Überblick über die Daten.

# mal deskriptiv anschauen
df %>%
  group_by(Zeitpunkt) %>%
  summarise(M = mean(Wert),
            SD = sd(Wert),
            Q1 = quantile(Wert, probs = 0.25),
            Q2 = quantile(Wert, probs = 0.5),
            Q3 = quantile(Wert, probs = 0.75),
            IQR = IQR(Wert),
            N = n())
## # A tibble: 3 × 8
##   Zeitpunkt     M    SD    Q1    Q2    Q3   IQR     N
##   <fct>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 t1         9.88 2.14   8.41  9.77 11.3   2.91   200
## 2 t2         7.99 0.935  7.24  8.09  8.60  1.36   200
## 3 t3         6.92 0.868  6.31  6.89  7.51  1.21   200

Das geht auch graphisch:

# mal graphisch anschauen
boxplot(df$Wert ~ df$Zeitpunkt)

ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt)

# mit Regressionsgerade
ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt) +
  geom_smooth(aes(x = as.numeric(Zeitpunkt)), method = 'lm', se = F)
## `geom_smooth()` using formula = 'y ~ x'

Die Plots legen nahe, dass es einen Unterschied in den Zeitpunktgruppen gibt.

Wir überprüfen mittels Shapiro-Test und QQ-Plot, ob die Daten normalverteilt sind.

## qq-Plot (Normalverteilung)
ggpubr::ggqqplot(df, x="Wert", 
                 facet.by = "Zeitpunkt", 
                 color = "Zeitpunkt",
                 scales = "free")

# Teste "Wert" auf Normalverteilung
shapiro.test(df$Wert[df$Zeitpunkt=="t1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t1"]
## W = 0.99364, p-value = 0.5481
shapiro.test(df$Wert[df$Zeitpunkt=="t2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t2"]
## W = 0.98766, p-value = 0.08025
shapiro.test(df$Wert[df$Zeitpunkt=="t3"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t3"]
## W = 0.99116, p-value = 0.2625

Die Testergebnisse sind nicht signifikant, also liegt Normalverteilung vor.

Wir dürfen eine ANOVA mit Messwiederholung rechnen. Mit R-Hausmittel könnten wir so vorgehen:

# ANOVA mit Messwiederholung 
fit <- aov(Wert ~ Zeitpunkt + Error(Name), data=df)
summary(fit)
## 
## Error: Name
##            Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 199  395.8   1.989               
## 
## Error: Within
##            Df Sum Sq Mean Sq F value Pr(>F)    
## Zeitpunkt   2  898.0   449.0   212.3 <2e-16 ***
## Residuals 398  841.6     2.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Das Ergebnis ist signifikant. Mit der Funktion pairwise.t.test() können wir nun schauen, welche Gruppenunterschiede es konkret gibt:

# geht auch mit R-Hausmitteln so:
pairwise.t.test(df$Wert, df$Zeitpunkt, 
                p.adjust="bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$Wert and df$Zeitpunkt 
## 
##    t1      t2     
## t2 < 2e-16 -      
## t3 < 2e-16 1.1e-12
## 
## P value adjustment method: bonferroni

Die Unterschiede sind zwischen allen Gruppen signifikant.

Das Zusatzpaket rstatix bietet mit der Funktion anova_test() eine alternative Vorgehensweise.

# ANOVA mit Messwiederholung 
fit <- rstatix::anova_test(data = df, dv = Wert, wid = Name,
                           within = Zeitpunkt, effect.size = "pes")

# auf Sphärizität prüfen
fit$`Mauchly's Test for Sphericity`
##      Effect     W        p p<.05
## 1 Zeitpunkt 0.668 4.29e-18     *

Mauchly’s Test ist signifikant, das heisst, es liegt keine Sphärizität vor. R hat in diesem Falle bereits die Freiheitsgrade automatisch korrigiert (mittels Greenhouse-Geisser-Verfahren).

Schauen wir uns die ANOVA an:

rstatix::get_anova_table(fit)
## ANOVA Table (type III tests)
## 
##      Effect DFn    DFd      F        p p<.05   pes
## 1 Zeitpunkt 1.5 298.73 212.34 2.18e-48     * 0.516

Die ANOVA ist signifikant. Jetzt können wir analysieren, welche konkreten Gruppenunterschiede vorliegen.

# post-hoc Analyse: paarweise Gruppenvergleiche
df %>%
  rstatix::pairwise_t_test(Wert ~ Zeitpunkt, paired = TRUE,
                           p.adjust.method = "bonferroni") %>%
  as.data.frame()
##    .y. group1 group2  n1  n2 statistic  df        p    p.adj p.adj.signif
## 1 Wert     t1     t2 200 200  11.24384 199 5.02e-23 1.51e-22         ****
## 2 Wert     t1     t3 200 200  18.28475 199 1.79e-44 5.37e-44         ****
## 3 Wert     t2     t3 200 200  11.26870 199 4.23e-23 1.27e-22         ****

Die Unterschiede sind zwischen allen Gruppen signifikant.

Die Effektstärke kann wie folgt berechnet werden.

# Effektstärke (nur solche, die post-hoc signifikant sind, interpretieren)
df %>%
  rstatix::cohens_d(Wert ~ Zeitpunkt, paired = TRUE) %>%
  as.data.frame()
##    .y. group1 group2   effsize  n1  n2 magnitude
## 1 Wert     t1     t2 0.7950598 200 200  moderate
## 2 Wert     t1     t3 1.2929274 200 200     large
## 3 Wert     t2     t3 0.7968171 200 200  moderate

Achtung, nur solche Werte interpretieren, die oben signifkant waren. In unserem Beispiel sind alle Werte signifikant, darum dürfen wir alle Effekte interpretieren.

31.3.2 Friedman ANOVA

Sind die Werte nicht normalverteilt, kann ein Friedman-ANOVA gerechnet werden. Wir nutzen den Testdatensatz MarioANOVA.rda. Er enthält Informationen über 47 Super Mario Charaktere, deren Fehlerquote beim Spielen eines schweren Levels zu 3 Zeitpunkten erfasst wurde.

load(url("https://www.produnis.de/R/MarioANOVA.rda"))

# anschauen
head(MarioANOVA)
## # A tibble: 6 × 8
##   Name             Alter Kingdom          Geschlecht BadGuy    t1    t2    t3
##   <fct>            <dbl> <fct>            <fct>      <lgl>  <dbl> <dbl> <dbl>
## 1 Mario               35 Mushroom Kingdom männlich   FALSE  0.372 0.472 0.257
## 2 Luigi               35 Mushroom Kingdom männlich   FALSE  1.27  1.34  1.14 
## 3 Prinzessin Peach    25 Mushroom Kingdom weiblich   FALSE  6.15  1.40  1.98 
## 4 Bowser              45 Bowser's Castle  männlich   TRUE   4.85  0.648 1.27 
## 5 Toad                20 Mushroom Kingdom männlich   FALSE  1.34  0.174 1.77 
## 6 Toadette            19 Mushroom Kingdom weiblich   FALSE  0.727 0.413 2.09

In den Variablen t1 bis t3 liegen die Werte der 3 Messzeitpunkte. Auch hier überführen wir die Daten ins long table-Format

df <- MarioANOVA %>%
  select(Name, t1:t3) %>%
  pivot_longer(cols=t1:t3,
               names_to = "Zeitpunkt",
               values_to = "Wert")

Schauen wir uns die Daten zunächst deskriptiv an.

# mal deskriptiv anschauen
df %>%
  group_by(Zeitpunkt) %>%
  summarise(M = mean(Wert),
            SD = sd(Wert),
            Q1 = quantile(Wert, probs = 0.25),
            Q2 = quantile(Wert, probs = 0.5),
            Q3 = quantile(Wert, probs = 0.75),
            IQR = IQR(Wert),
            N = n())
## # A tibble: 3 × 8
##   Zeitpunkt     M    SD    Q1    Q2    Q3   IQR     N
##   <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 t1        2.75  4.61  0.750 1.34   2.46 1.71     47
## 2 t2        0.762 0.585 0.356 0.509  1.15 0.791    47
## 3 t3        0.901 0.701 0.266 0.650  1.34 1.07     47

Das geht auch graphisch wie oben beschrieben:

# mal graphisch anschauen
boxplot(df$Wert ~ df$Zeitpunkt)

ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt)

# mit Regressionsgerade
ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt) +
  geom_smooth(aes(x = as.numeric(Zeitpunkt)), method = 'lm', se = F)

Nun prüfen wir mittels QQ-Plot und Shapiro-Test auf Normalverteilung.

## qq-Plot (Normalverteilung)
ggpubr::ggqqplot(df, x="Wert", 
                 facet.by = "Zeitpunkt", 
                 color = "Zeitpunkt",
                 scales = "free")

# Teste "Wert" auf Normalverteilung
shapiro.test(df$Wert[df$Zeitpunkt=="t1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t1"]
## W = 0.52189, p-value = 3.729e-11
shapiro.test(df$Wert[df$Zeitpunkt=="t2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t2"]
## W = 0.89427, p-value = 0.0004751
shapiro.test(df$Wert[df$Zeitpunkt=="t3"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t3"]
## W = 0.90101, p-value = 0.0007756

Die Werte für t1 bis t3 sind signifikant, also nicht normalverteilt.

Wir führen daher eine Friedman-ANOVA mit Messwiederholung durch.

# Friedman-Test
rstatix::friedman_test(Wert ~ Zeitpunkt | Name, data=df)
## # A tibble: 1 × 6
##   .y.       n statistic    df       p method       
## * <chr> <int>     <dbl> <dbl>   <dbl> <chr>        
## 1 Wert     47      12.8     2 0.00165 Friedman test

Das Ergebnis ist signifikant, es gibt also einen Unterschied.

# post-hoc-test
df %>%
  rstatix::wilcox_test(Wert ~ Zeitpunkt, paired=T,
                       p.adjust.method = "bonferroni")
## # A tibble: 3 × 9
##   .y.   group1 group2    n1    n2 statistic        p    p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl> <chr>       
## 1 Wert  t1     t2        47    47       917 0.000102 0.000306 ***         
## 2 Wert  t1     t3        47    47       882 0.000536 0.002    **          
## 3 Wert  t2     t3        47    47       506 0.546    1        ns

Der Unterschied existiert zwischen t1 und t2 sowie zwischen t1 und t3.

Die Effktstärke wird wie folgt berechnet:

# Effektstärke
rstatix::friedman_effsize(Wert ~ Zeitpunkt | Name, data = df)
## # A tibble: 1 × 5
##   .y.       n effsize method    magnitude
## * <chr> <int>   <dbl> <chr>     <ord>    
## 1 Wert     47   0.136 Kendall W small

Der Effekt des Gesamtmodells ist eher gering. Schauen wir uns die einzelnen Gruppeneffekte an.

Achtung, nur solche Werte interpretieren, die oben signifkant waren. In unserem Beispiel ist das der Unterschied von t1 zu t2 sowie von t1 zu t3..

# nur die signifikanten anschauen (t1 -> t2 und t1 -> t3)
df %>%
  rstatix::wilcox_effsize(Wert ~ Zeitpunkt, paired = TRUE)
## # A tibble: 3 × 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 Wert  t1     t2      0.545     47    47 large    
## 2 Wert  t1     t3      0.491     47    47 moderate 
## 3 Wert  t2     t3      0.0895    47    47 small

Der Effekt von t1 nach t2 ist “groß”, der von t1 nach t3 “moderat”.

31.3.3 einfaktorielle 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 Beobachtung als fehlend gelöscht

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

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

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

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

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

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

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

31.4 Fallzahlkalkulation

31.4.1 klinische Studien

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

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

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

# Lade Paket
library(pwr)

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

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

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

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

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

Sie benötigen 1571 Probanden pro Untersuchungsgruppe.

31.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 - nach oben sowie nach unten - sein darf, für gewöhnlich 5% bzw. 0.05), das Konfidenzniveau, die Standardabweichung S innerhalb der Grundgesamtheit (füre gewöhnlich 50% bzw. 0.5) 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%, Fehlerbereich 5%.
# In der Grundgesamtheit hat der Wert eine
# Standardabweichung von S=0.5
samplingbook::sample.size.mean(e=0.05, S=0.5, level=0.95 )
## 
## sample.size.mean object: Sample size for mean estimate
## Without finite population correction: N=Inf, precision e=0.05 and standard deviation S=0.5
## 
## Sample size needed: 385

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

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

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

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