Zu den Daten

Die folgende \(2 \times 2\)-Tabelle zeigt die Daten einer Studie, welche durchgeführt wurde, um die Wirksamkeit von Aspirin zur Prävention eines frühzeitigen Todes nach einem Herzinfarkt zu evaluieren. Hier findest du weitere Informationen zu den Daten.

##          outcome
## treatment survived death
##   placebo      354    52
##   aspirin      725    85

Du kannst dir diese Tabelle in R selbst erstellen:

# Matrix anhand der gegebenen Werte erstellen
m <- matrix(c(354, 52, 725, 85), byrow = TRUE, ncol = 2)

# Zeilen und Spalten benennen
rownames(m) <- c("placebo", "aspirin")
colnames(m) <- c("survived", "died")

# aus einer Matrix eine Tabelle machen
tab <- as.table(m)
tab
##         survived died
## placebo      354   52
## aspirin      725   85

Manchmal möchte man auf Grundlage einer \(2 \times 2\)-Tabelle Einzeldaten mit denselben Zellhäufigkeiten erzeugen (z. B. um Grafiken zu erstellen). Dazu gibt es eine praktische Funktion aus dem Paket rstatix.

library(rstatix)
df.aspirin <- counts_to_cases(tab)

# Variablen des erstellten Dataframes benennen
names(df.aspirin) <- c("gruppe", "outcome")
head(df.aspirin)
##    gruppe  outcome
## 1 placebo survived
## 2 placebo survived
## 3 placebo survived
## 4 placebo survived
## 5 placebo survived
## 6 placebo survived

Proportionen pro Gruppe

Aufgaben

Verwende für die Aufgaben die oben gezeigte Tabelle tab.

  1. Berechne anhand der Normalapproximation für beide Gruppen jeweils die Proportion \(\hat{\pi}\) der Personen, die gestorben sind, sowie ein zugehöriges 95%-Konfidenzintervall. Der Standardfehler ist

    \[ SE_{\hat{\pi}} = \sqrt{\frac{\hat{\pi}(1-\hat{\pi})}{n}}, \]

  2. Vergleiche die manuellen Berechnungen mit dem Resultat der in R implementierten Funktion prop.test(). Interpretiere!

Lösungen

Wir stützen die Berechnungen auf die folgenden Angaben:

tab
##         survived died
## placebo      354   52
## aspirin      725   85

Anhand der Tabelle berechnen wir zuerst die Proportion, den Standardfehler und schliesslich das 95%-Konfidenzintervall. Dabei gilt:

\[ SE_{\hat{\pi}} = \sqrt{\frac{\hat{\pi}(1-\hat{\pi})}{n}}. \]

cases <- tab[, 2] # Todesfälle pro Gruppe
n <- tab[, 1] + tab[, 2] # n pro Gruppe
prop <- cases/n # Proportion pro Gruppe
se.prop <- sqrt((prop * (1-prop))/n) # Standardfehler pro Gruppe
CI.lower <- prop-1.96*se.prop
CI.upper <- prop+1.96*se.prop
cbind(CI.lower, CI.upper)
##           CI.lower  CI.upper
## placebo 0.09557231 0.1605853
## aspirin 0.08383222 0.1260443

Mit prop.test() erhalten wir die folgenden Resultate.

Placebo:

prop.test(x = cases["placebo"], n = n["placebo"])
## 
##  1-sample proportions test with continuity correction
## 
## data:  cases["placebo"] out of n["placebo"], null probability 0.5
## X-squared = 223.16, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.09793556 0.16545553
## sample estimates:
##         p 
## 0.1280788

Aspirin:

prop.test(x = cases["aspirin"], n = n["aspirin"])
## 
##  1-sample proportions test with continuity correction
## 
## data:  cases["aspirin"] out of n["aspirin"], null probability 0.5
## X-squared = 504.1, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.08510476 0.12860969
## sample estimates:
##         p 
## 0.1049383

Bemerkung: Für einen Binomialanteil gibt es verschiedene Konfidenzintervalle. Das oben mithilfe der Normalapproximation berechnete Intervall heisst Wald-Intervall und wird in der folgenden Tabelle als asymptotic bezeichnet. Für einen einzelnen Anteil berechnet prop.test() ein Score-Intervall (Wilson-Intervall), standardmässig mit Stetigkeitskorrektur. Die Methoden werden später genauer behandelt.

library(binom)
binom.confint(cases["aspirin"],n["aspirin"],methods=c("asymptotic","prop.test"))
##       method  x   n      mean      lower     upper
## 1 asymptotic 85 810 0.1049383 0.08383261 0.1260439
## 2  prop.test 85 810 0.1049383 0.08510476 0.1286097

Die Proportionen betragen also 12.81 % (Placebo) bzw. 10.49 % (Aspirin).

Da die Gruppengrössen und die Anzahlen der Todesfälle und Überlebenden ausreichend gross sind, liefern das Wald-Intervall und die Funktion prop.test() hier ähnliche Ergebnisse. Wir sehen, dass sich die 95%-Konfidenzintervalle der beiden Gruppen überlappen. Das bedeutet jedoch nicht zwingend, dass die Differenz der Proportionen statistisch nicht signifikant ist, da der Standardfehler der Differenz ein anderer ist als die Standardfehler der Einzelproportionen.

Für einen einzelnen Anteil testet prop.test() standardmässig, ob sich die Proportion von 0.5 unterscheidet (\(H_0: \pi = 0.5\)). In diesem Fall interessiert uns dieser Test (und folglich sein P-Wert) nicht. Mit dem Argument p = … könnte man \(H_0\) individuell anpassen. Es kommt bei Proportionen jedoch häufig vor, dass man primär an der Präzision der Schätzung interessiert ist, die sich anhand eines 95%-Konfidenzintervalls beurteilen lässt.

Proportionsdifferenzen

Aufgaben

Es wird weiterhin mit der \(2 \times 2\)-Tabelle von oben gerechnet.

  1. Berechne die Proportionsdifferenz zwischen der Placebo- und der Aspiringruppe und berechne mithilfe der Normalapproximation ein 95%-Konfidenzintervall und einen P-Wert für diese Differenz.

  2. Wiederhole die Berechnungen mithilfe der Funktion prop.test(). Vergleiche und interpretiere die Ergebnisse.

Lösungen

Wir berechnen die Differenz

\[ \hat{\pi}_1-\hat{\pi}_2 \]

und den Standardfehler für diese Differenz

\[ SE_{\hat{\pi_1}-\hat{\pi_2}} = \sqrt{SE_{\hat{\pi_1}}^2 + SE_{\hat{\pi_2}}^2}. \]

Die Reihenfolge der Gruppen bestimmt das Vorzeichen der Differenz und die Orientierung des Konfidenzintervalls. Interpretiere die Ergebnisse entsprechend der gewählten Reihenfolge. Der zweiseitige P-Wert bleibt bei einer Umkehrung der Reihenfolge unverändert.

prop.diff <- prop["placebo"] - prop["aspirin"]
se.prop.diff <- sqrt(se.prop["placebo"]^2 + se.prop["aspirin"]^2)
CI.diff.lower <- prop.diff-1.96*se.prop.diff
CI.diff.upper <- prop.diff+1.96*se.prop.diff
cbind(CI.diff.lower, CI.diff.upper)
##         CI.diff.lower CI.diff.upper
## placebo   -0.01561689    0.06189798

Den zweiseitigen P-Wert für \(H_0: \pi_1 - \pi_2 = 0\) berechnen wir mit einem Wald-Test anhand der standardnormalapproximierten Teststatistik:

z <- prop.diff/se.prop.diff
p <- 2*pnorm(-abs(z))
p
##   placebo 
## 0.2419047

Die Funktion prop.test() liefert ähnliche Ergebnisse. Für die Differenz zweier Anteile berechnet sie ein Wald-Intervall mit Stetigkeitskorrektur. Ihr Test verwendet unter der Nullhypothese einen gemeinsam geschätzten Anteil; der obige Wald-Test verwendet dagegen die separat geschätzten Anteile zur Berechnung des Standardfehlers.

Ohne Stetigkeitskorrektur entspricht der Test von prop.test() hier einem Score-Test, dessen quadratische Teststatistik exakt mit der Pearson-Chi-Quadrat-Teststatistik übereinstimmt. Das Konfidenzintervall der Anteilsdifferenz wird dagegen nach dem Wald-Verfahren berechnet; nur für einen einzelnen Anteil liefert prop.test() ein Score- beziehungsweise Wilson-Intervall. Mit correct = FALSE entfällt die Stetigkeitskorrektur, der Unterschied zwischen den Standardfehlern des Score- und des manuellen Wald-Tests bleibt jedoch bestehen. Die Herleitung erläutert der folgende Exkurs; die Grundlagen der Likelihood-Inferenz behandeln wir später.

prop.test(x = cases, n = n)
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  cases out of n
## X-squared = 1.2264, df = 1, p-value = 0.2681
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.01746498  0.06374608
## sample estimates:
##    prop 1    prop 2 
## 0.1280788 0.1049383

Der beobachtete Anteil der verstorbenen Personen ist in der Aspiringruppe um etwa 2.3 Prozentpunkte niedriger. Für die Differenz Placebo minus Aspirin liefert prop.test() ein 95%-Konfidenzintervall von etwa [-0.017; 0.064], also von −1.7 bis 6.4 Prozentpunkten. Unter den Modellannahmen sind damit sowohl negative als auch positive Differenzen mit den Daten vereinbar. Das Intervall enthält 0; der zweiseitige Test ist auf dem 5%-Niveau nicht statistisch signifikant. Der P-Wert von 26.81 % ist unter \(H_0\) die Wahrscheinlichkeit einer mindestens so grossen Teststatistik wie der beobachteten. Er ist keine Wahrscheinlichkeit dafür, dass \(H_0\) zutrifft. Ein 95%-Konfidenzverfahren überdeckt bei wiederholter Stichprobenziehung unter denselben Modellannahmen in ungefähr 95 % der Fälle den wahren Parameter.