pacman::p_load(tidyverse, magrittr, effectsize, rstatix,
scales, parameters, conflicted)
conflicts_prefer(stats::chisq.test)
conflicts_prefer(stats::fisher.test)35 Der \(\mathcal{X}^2\)-Test
Letzte Änderung am 02. April 2024 um 09:52:43
“You’re so square, baby, I don’t care.” — Elvis Presley
Was macht der \(\mathcal{X}^2\)-Test?
Ein \(\mathcal{X}^2\)-Test vergleicht die Anteile zweier oder mehrerer Gruppen. Da Anteile Wahrscheinlichkeiten sind, vergleicht der \(\mathcal{X}^2\)-Test damit auch Wahrscheinlichkeiten.
Der \(\mathcal{X}^2\)-Test wird häufig verwendet, wenn wir zwei Faktoren mit jeweils zwei Leveln miteinander vergleichen wollen. Das heißt wir haben zum Beispiel unseren Faktor animal mit den beiden Leveln cat und dog. Wir schauen uns jetzt den Infektionsstatus mit Flöhen auf den Tieren an. Wir erhalten wiederum einen Faktor infected mit zwei Leveln yes und no. Wir sind bei dem \(\mathcal{X}^2\)-Test nicht auf nur Faktoren mit zwei Leveln eingeschränkt. Traditionell wird aber versucht ein 2x2 Setting zu erreichen.
Im Bereich der Agrarwissenschaften kommt der \(\mathcal{X}^2\)-Test eher selten vor. Im Bereich der Humanwissenschaften und vor allem der Epidemiologie ist der \(\mathcal{X}^2\)-Test weit verbreitet.
Das eigentlich besondere an dem \(\mathcal{X}^2\)-Test ist gra nicht mal der Test selber sondern die Datenstruktur die der \(\mathcal{X}^2\)-Test zugrunde liegt: der Vierfeldertafel oder 2x2 Kreuztabelle. Wir werden diese Form von Tabelle noch später im maschinellen Lernen und in der Testdiagnostik wiederfinden.
Ein Wort zur Klausur
Wir rechnen den \(\mathcal{X}^2\)-Test per Hand in der Klausur. Teilweise füllen wir die \(2x2\) Tabelle aus, teilweise berechnen wir nur die Randsummen.
Bitte schau dir die Aufgaben in den gesammelten Klausurfragen auf GitHub an um eine Idee zu haben, welche Fragen zum \(\mathcal{X}^2\)-Test drankommen.
Wenn kein \(\chi^2_{\alpha=5\%}\) in der Klausur gegeben ist, setzen wir \(\chi^2_{\alpha=5\%} = 3.84\).
35.1 Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
Am Ende des Kapitels findest du nochmal den gesamten R Code in einem Rutsch zum selber durchführen oder aber kopieren.
35.2 Daten für den \(\mathcal{X}^2\)-Test
Wie eben schon benannt schauen wir uns für den \(\mathcal{X}^2\)-Test eine Vierfeldertafel oder aber 2x2 Kreuztabelle an. In Tabelle 35.1 sehen wir eine solche 2x2 Kreuztabelle. Da wir eine Mindestanzahl an Zellbelegung brauchen um überhaupt mit dem \(\mathcal{X}^2\)-Test rechnen zu können, nutzen wir hier gleich aggrigierte Beispieldaten. Wir brauchen mindestens fünf Beobachtungen je Zelle, dass heißt mindestens 20 Tiere. Da wir dann aber immer noch sehr wenig haben, ist die Daumenregel, dass wir etwa 30 bis 40 Beobachtungen brauchen. In unserem Beispiel schauen wir uns 65 Tiere an.
| Infected | ||||
| Yes (1) | No (0) | |||
| Animal | Dog | \(23_{\;\Large a}\) | \(10_{\;\Large b}\) | \(\mathbf{a+b = 33}\) |
| Cat | \(18_{\;\Large c}\) | \(14_{\;\Large d}\) | \(\mathbf{c+d = 32}\) | |
| \(\mathbf{a+c = 41}\) | \(\mathbf{b+d = 24}\) | \(n = 65\) |
In der Tabelle sehen wir, dass in den zeieln die Level des Faktors animal angegeben sind und in den Spalten die Level des Faktors infected. Wir haben somit \(23\) Hunde, die mit Flöhen infiziert sind, dann \(10\) Hunde, die nicht mit Flöhen infiziert sind. Auf der Seite der Katzen haben wir \(18\) Katzen, die infiziert sind und \(14\) Katzen, die keine Flöhe haben. An den Rändern stehen die Randsummen. Wir haben \(33\) Hunde und \(32\) Katzen sowie \(41\) infizierte Tiere und \(24\) nicht infizierte Tiere. Somit haben wir dann in Summe \(n = 65\) Tiere. Diese Form der Tabelle wird uns immer wieder begegnen.
Bevor wir jetzt diese 2x2 Kreuztabelle verwenden, müssen wir uns nochmal überlegen, welchen Schluss wir eigentlich über die Nullhypothese machen. Wie immer können wir nur die Nullhypothese ablehnen. Daher überlegen wir uns im Folgenden wie die Nullhypothese in dem \(\mathcal{X}^2\)-Test aussieht. Dann bilden wir anhand der Nullhypothese noch die Alternativehypothese.
35.3 Hypothesen für den \(\mathcal{X}^2\)-Test
Der \(\mathcal{X}^2\)-Test betrachtet die Zellbelegung gegeben den Randsummen um einen Unterschied nachzuweisen. Daher haben wir die Nullhypothese als Gleichheitshypothese. In unserem Beispiel lautet die Nullhypothese, dass die Zahlen in den Zellen gegeben der Randsummen gleich sind. Wir betrachten hier nur die Hypothesen in Prosa und die mathematischen Hypothesen. Es ist vollkommen ausreichend, wenn du die Nullhypothese des \(\mathcal{X}^2\)-Test nur in Prosa kennst.
\[ H_0: \; \mbox{Zellbelegung sind gleichverteilt gegeben der Randsummen} \]
Die Alternative lautet, dass sich die Zahlen in den Zellen gegeben der Randsummen unterscheiden.
\[ H_A: \; \mbox{Zellbelegung sind nicht gleichverteilt gegeben der Randsummen} \]
Wir schauen uns jetzt einmal den \(\mathcal{X}^2\)-Test theoretisch an bevor wir uns mit der Anwendung des \(\mathcal{X}^2\)-Test in R beschäftigen.
35.4 \(\mathcal{X}^2\)-Test theoretisch
In Tabelle 35.1 von oben hatten wir die beobachteten Werte. Das sind die Zahlen, die wir in unserem Experiment erhoben und gemessen haben. Der \(\mathcal{X}^2\)-Test vergleicht nun die beobachteten Werte mit den anhand der Randsummen zu erwartenden Werte. Daher ist die Formel für den \(\mathcal{X}^2\)-Test wie folgt.
\[ \chi^2_{calc} = \sum\cfrac{(O - E)^2}{E} \]
mit
- \(O\) für die beobachteten Werte
- \(E\) für die nach den Randsummen zu erwartenden Werte
In Tabelle 35.2 kannst du sehen wie wir anhand der Randsummen die erwartenden Zellbelegungen berechnen. Hierbei können auch krumme Zahlen rauskommen. Wir würden keinen Unterschied zwischen Hunde und Katzen gegeben deren Infektionsstatus erwarten, wenn die Abweichungen zwischen den beobachteten Werten und den zu erwartenden Werten klein wären. Wir berechnen nun die zu erwartenden Werte indem wir die Randsummen der entsprechenden Zelle multiplizieren und durch die Gesamtanzahl teilen.
| Infected | ||||
| Yes (1) | No (0) | |||
| Animal | Dog | \(\cfrac{41 \cdot 33}{65} = 20.82\) | \(\cfrac{24 \cdot 33}{65} = 12.18\) | \(\mathbf{33}\) |
| Cat | \(\cfrac{41 \cdot 32}{65} = 20.18\) | \(\cfrac{24 \cdot 32}{65} = 11.82\) | \(\mathbf{32}\) | |
| \(\mathbf{41}\) | \(\mathbf{24}\) | \(n = 65\) |
Wir können dann die Formel für den \(\mathcal{X}^2\)-Test entsprechend ausfüllen. Dabei ist wichtig, dass die Abstände quadriert werden. Das ist ein Kernkonzept der Statistik, Abstände bzw. Abweichungen werden immer quadriert.
\[\begin{aligned} \chi^2_{calc} &= \cfrac{(23 - 20.82)^2}{20.82} + \cfrac{(10 - 12.18)^2}{12.18} + \\ &\phantom{=}\;\; \cfrac{(18 - 20.18)^2}{20.18} + \cfrac{(14 - 11.82)^2}{11.82} = 1.25 \end{aligned}\]Es ergibt sich ein \(\chi^2_{calc}\) von \(1.25\) mit der Regel, dass wenn \(\chi^2_{calc} \geq \chi^2_{\alpha=5\%}\) die Nullhypothese abgelehnt werden kann. Mit einem \(\chi^2_{\alpha=5\%} = 3.84\) können wir die Nullhypothese nicht ablehnen. Es besteht kein Zusammenhang zwischen den Befall mit Flöhen und der Tierart. Oder anders herum, Hunde und Katzen werden gleich stark mit Flöhen infiziert.
Wie stark ist nun der beobachtete Effekt? Wir konnten zwar die Nullhypothese nicht ablehnen, aber es wäre auch von Interesse für die zukünftige Versuchsplanung, wie stark sich die Hunde und Katzen im Befall mit Flöhen unterscheiden. Wir haben nun die Wahl zwischen zwei statistischen Maßzahlen für die Beschreibung eines Effektes bei einem \(\mathcal{X}^2\)-Test.
Zum einen Cramers \(V\), das wir in etwa interpretieren wie eine Korrelation und somit auch einheitslos ist. Wenn Cramers \(V\) gleich 0 ist, dann haben wir keinen Unterschied zwischen den Hunden und Katzen gegeben dem Flohbefall. Bei einem Cramers \(V\) von 0.5 haben wir einen sehr starken Unterschied zwischen dem Flohbefall zwischen den beiden Tierarten. Der Vorteil von Cramers V ist, dass wir Cramers \(V\) auch auf einer beliebig großen Kreuztabelle berechnen können. Die Formel ist nicht sehr komplex.
\[ V = \sqrt{\cfrac{\mathcal{X}^2/n}{\min(c-1, r-1)}} \]
mit
- \(\mathcal{X}^2\) gleich der Chi-Quadrat-Statistik
- \(n\) gleich der Gesamtstichprobengröße
- \(r\): Anzahl der Reihen
- \(c\): Anzahl der Spalten
\[ V = \sqrt{\cfrac{1.26/65}{1}} = 0.14 \]
Wir setzen also den Wert \(\chi^2_{calc}\) direkt in die Formel ein. Wir wissen ja auch, dass wir \(n = 65\) Tiere untersucht haben. Da wir eine 2x2 Kreuztabelle vorliegen haben, haben wir \(r = 2\) und \(c = 2\) und somit ist das Minimum von \(r-1\) und \(c-1\) gleich \(1\). Wir erhalten ein Cramers \(V\) von \(0.14\) was für einen schwachen Effekt spricht. Wir können uns auch grob an der folgenden Tabelle der Effektstärken für Carmers \(V\) orientieren.
| schwach | mittel | stark | |
|---|---|---|---|
| Cramers \(V\) | 0.1 | 0.3 | 0.5 |
Zum anderen haben wir noch die Möglichkeit die Odds Ratios oder das Chancenverhältnis zu berechnen. Die Odds Ratios lassen sich direkter als Effekt interpretieren als Cramers \(V\) haben aber den Nachteil, dass wir die Odds Ratios nur auf einer 2x2 Kreuztabelle berechnen können. Wichtig bei der Berechnung der Odds Ratios und der anschließenden Interpretation ist die obere Zeile der 2x2 Kreuztabelle. Die obere Zeile ist der Bezug für die Interpretation. Wir nutzen folgende Formel für die Berechnung der Odds Ratios.
\[ \mbox{Odds Ratio} = OR = \cfrac{a\cdot d}{b \cdot c} = \cfrac{23\cdot 14}{10 \cdot 18} = 1.79 \]
Es ist zwingend notwendig für die folgenden Interpretation der Odds Ratios, dass in den Spalten links die \(ja\) Spalte steht und rechts die \(nein\) Spalte. Ebenso interpretieren wie die Odds Ratios im Bezug zur oberen Zeile. In unserem Fall ist also die Chance sich mit Flöhen zu infizieren bei Hunden 1.79 mal größer als bei Katzen. Diese Interpretation ist nur korrekt, wenn die 2x2 Kreuztabelle wie beschrieben erstellt ist!
35.5 \(\mathcal{X}^2\)-Test in R
Der \(\mathcal{X}^2\)-Test wird meist in anderen Funktionen in R verwendet und nicht direkt. Wenn du Fragen dazu hast, schreib mir einfach eine Mail.
Wenn wir den \(\mathcal{X}^2\)-Test in R rechnen wollen nutzen wir die Funktion chisq.test(), die eine Matrix von Zahlen verlangt. Dies ist etwas umständlich. Wir müssen nur beachten, dass wir die Matrix so bauen, wie wir die Matrix auch brauchen. Deshalb immer mal doppelt schauen, ob deine Matrix auch deinen beobachteten Werten entspricht.
mat <- matrix(c(23, 10, 18, 14),
byrow = TRUE, nrow = 2,
dimnames = list(animal = c("dog", "cat"),
infected = c("yes", "no")))
chisq.test(mat, correct = FALSE)
Pearson's Chi-squared test
data: mat
X-squared = 1.2613, df = 1, p-value = 0.2614Als ein mögliches Effektmaß können wir Cramers \(V\) berechnen. Wir nutzen hierzu die Funktion cramers_v(). Auf einer reinen 2x2 Kreuztabelle wäre aber Pearsons \(\phi\) durch die Funktion phi() vorzuziehen. Siehe dazu auch \(\phi\) and Other Contingency Tables Correlations auf der Hilfeseite des R Paketes {effectsize}. Wir bleiben hier dann aber bei Cramers \(V\).
cramers_v(mat) Cramer's V (adj.) | 95% CI
--------------------------------
0.06 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Wir sehen, dass der Effekt mit einem \(V = 0.06\) schwach ist. Ein Wert von 0 bedeutet keine Assoziation und ein Wert von 1 einen maximalen Zusammenhang. Wir können die Werte von \(V\) wie eine Korrelation interpretieren. Der in R berechnete Wert unterscheidet sich von unseren händsichen berechneten Wert, da wir hier mit dem adjustierten Cramers \(V\) rechnen. Wir würden auch in der Anwendung den adjustierten Wert verwenden, aber für das Verständnis reicht der händisch berechnete Wert.
Haben wir eine geringe Zellbelegung von unter 5 in einer der Zellen der 2x2 Kreuztabelle, dann verwenden wir den Fisher Exakt Test. Der Fisher Exakt Test hat auch den Vorteil, dass wir direkt die Odds Ratios wiedergegeben bekommen. Wir können auch den Fisher Exakt Test rechnen, wenn wir viele Beobachtungen pro Zelle haben und einfach an die Odds Ratios rankommen wollen. Der Unterschied zwischen dem klassischen \(\mathcal{X}^2\)-Test und dem Fisher Exakt Test ist in der praktischen Anwendung nicht so groß.
fisher.test(mat)
Fisher's Exact Test for Count Data
data: mat
p-value = 0.3102
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.5756546 5.6378438
sample estimates:
odds ratio
1.77279 Wir sehen auch hier den nicht signifikanten \(p\)-Wert sowie eine Odds Ratio von 1.77. Hunde haben aso eine um 1.77 höhere Chance sich mit Flöhen zu infizieren.
35.6 Test auf gleiche oder gegebene Anteile
Neben dem Test auf absolute Anteile, wie der \(\mathcal{X}^2\)-Test es rechnet, wollen wir manchmal auch relative Anteile testen. Also haben wir nicht 8 kranke Erdbeeren gezählt sondern 8 kranke von 12 Erdbeeren. Damit sind dann 66% bzw. 0.66 kranke Erdbeeren vorhanden. Wir rechnen also mit Wahrscheinlichkeiten, also eng. proportions. Deshalb können wir hier auch die R Funktion prop.test() nutzen. Wichtig ist zu wissen, dass wir trotz allem erst die Anzahl an beschädigten Erdbeeren \(x\) sowie die absolute Anzahl an Erdbeeren \(n\) brauchen. Das wird jetzt aber gleich in dem Beispiel etwas klarer. Im Prinzip ist der prop.test() also eine andere Art den \(\mathcal{X}^2\)-Test zu rechnen.
35.6.1 Vergleich eines Anteils \(p_1\) gegen einen gegebenen Anteil \(p_0\)
Nehmen wir ein etwas konstruiertes Beispiel zur Erdbeerernte. Wir haben einen neuen Roboter entwickelt, der Erdbeeren erntet. Nun stellen wir fest, dass von 100 Erdbeeren 76 heile sind. Jetzt lesen wir im Handbuch, dass der Ernteroboter eigentlich 84% der Erdbeeren heile ernten sollte. Sind jetzt 76 von 100, also 76%, signifikant unterschiedlich von 84%? Oder können wir die Nullhypothese der Gleichheit zwischen den beiden Wahrscheinlichkeiten nicht ablehnen? Damit können wir den Test auf Anteile nutzen, wenn wir eine beobachtete Wahrscheinlichkeit oder Anteil \(p_1\) gegen eine gegebene Wahrscheinlichkeit oder Anteil \(p_0\) vergleichen wollen.
In unserem Fall haben wir mit \(p_1\) die Wahrscheinlichkeit vorliegen, dass die Erdbeeren in unserem Experiment heile sind. Wir wissen also, dass \(p_1 = \cfrac{x}{n} = \cfrac{76}{100} = 0.76 = 76\%\) ist. Damit ist die Wahrscheinlichkeit \(p_1\) auch die beobachte Wahrscheinlichkeit. Wir erwarten auf der anderen Seite die Wahrscheinlichkeit \(p_0 = 0.84 = 84\%\). Die Roboterfirma hat uns aj zugesagt, dass 84% der Erdbeeren heile bleiben sollen.
Wir berechnen jetzt einmal die beobachten Werte. Zum einen haben wir \(x=76\) heile Erdbeeren von \(n=100\) Erdbeeren gezählt. Damit ergeben sich dann \(x = 76\) heile Erdbeeren und \(24\) beschädigte Erdbeeren. Die Summe muss ja am Ende wieder 100 Erdbeeren ergeben.
\[ \begin{aligned} O &= \begin{pmatrix} x &|& n - x \end{pmatrix} = \begin{pmatrix} 76 &|& 100 - 76 \end{pmatrix} \\ & = \begin{pmatrix} 76 &|& 24 \end{pmatrix} \end{aligned} \]
Dann müssen wir uns noch die Frage stellen, welche Anzahl an heilen Erdbeeren hätten wir erwartet? In diesem Fall ja 84% heile Erdbeeren. Das macht dann bei 100 Erdbeeren \(0.84 \cdot 100 = 84\) heile Erdbeeren und \((1 - 0.84) \cdot 100 = 16\) beschädigte Erdbeeren.
\[ \begin{aligned} E &= \begin{pmatrix} n \cdot p &|& n \cdot (1 - p) \end{pmatrix} = \begin{pmatrix} 100 \cdot 0.84 &|& 100 \cdot (1 - 0.84) \end{pmatrix} \\ & = \begin{pmatrix} 84 &|& 16 \end{pmatrix} \end{aligned} \]
Jetzt müssen wir nur noch die beobachteten Anteile mit den zu erwarteten Anteilen durch die \(\mathcal{X}^2\)-Formel in Kontext setzten. Wir subtrahieren von jedem beobachten Anteil den zu erwartenden Anteil, quadrieren und addieren auf. Dann erhalten wir die \(\mathcal{X}^2\)-Statistik.
\[ \chi^2_{calc} = \sum\cfrac{(O-E)^2}{E} = \cfrac{(76 - 84)^2}{84} + \cfrac{(24 - 16)^2}{16} = 4.76 \]
Wir können diese einfache Berechnung dann auch schnell nochmal in R durchführen. Wir setzten dafür einfach x, n und p fest und berechnen dann die beobachten Anteile \(O\) sowwie die zu erwartenden Anteile \(E\). Wir berechnen dann hier auch den gleichen Wert für \(\mathcal{X}^2\)-Statistik.
x <- 76
n <- 100
p <- 0.84
O <- cbind(x, n - x)
E <- cbind(n * p, n * (1 - p))
sum((abs(O - E))^2/E)[1] 4.761905Und wie immer gibt es auch eine Funktion prop.test(), die uns ermöglicht einen beobachteten Anteil x/n zu einem erwarteten Anteil p zu vergleichen. Auch hier sehen wir, dass sich die \(\mathcal{X}^2\)-Statistik aus der R Funktion nicht von unser berechneten \(\mathcal{X}^2\)-Statistik unterscheidet.
prop.test(x = 76, n = 100, p = 0.84, correct = FALSE)
1-sample proportions test without continuity correction
data: 76 out of 100, null probability 0.84
X-squared = 4.7619, df = 1, p-value = 0.0291
alternative hypothesis: true p is not equal to 0.84
95 percent confidence interval:
0.6676766 0.8330867
sample estimates:
p
0.76 Was ist die Yates Korrektur, die wir mit correct = TRUE auswählen? Die Korrektur von Yates subtrahiert von jedem Summanden des Chi-Quadrat-Tests 0.5, so dass die Chi-Quadrat-Statistik geringer ausfällt.
35.6.2 Vergleich zweier Anteile \(p_1\) und \(p_2\)
Nun haben wir nicht einen Ernteroboter sondern zwei brandneue Robotertypen. Einmal die Marke Picky und einmal den Roboter Colly. Wir wollen jetzt wieder bestimmen, wie viel Erdbeeren bei der Ernte beschädigt werden. Erdbeeren sind ja auch ein sehr weiches Obst. Wir vergleichen hier also wieder zwei Anteile miteinander. Der Roboter Picky beschädigt 76 von 100 Erdbeeren und der Roboter Colly detscht 91 von 100 Erdbeeren an. Das sind ganz schön meise Werte, aber was will man machen, wir haben jetzt nur die beiden Roboter vorliegen. Die Frage ist nun, ob sich die beiden Roboter in der Häufigkeit der beschädigten Erdbeeren unterscheiden. Wir können hier eine 2x2 Kreuztabelle in der Tabelle 35.5 aufmachen und die jeweiligen Anteile berechnen. Picky beschädigt 76% der Erdbeeren und Colly ganze 91%.
| Damaged | ||||
| Yes (1) | No (0) | |||
| Robot | Picky | \(76\) | \(24\) | \(\mathbf{100}\) |
| Colly | \(91\) | \(9\) | \(\mathbf{100}\) | |
| \(\mathbf{167}\) | \(\mathbf{33}\) | \(n = 200\) |
Wir können die erwarteten Anteile jetzt wie schon bekannt berechnen oder aber wir nutzen folgende Formel um \(\hat{p}\), die Wahrscheinlichkeit für ein Ereignis, zu berechnen. Wir nutzen dann \(\hat{p}\) um die erwarteten Werte \(E\) aus zurechnen. Dafür addieren wir alle beobachteten \(x\) zusammen und teilen diese Summe durch die gesammte Anzahl an Beobachtungen.
\[ \hat{p} = \cfrac{\sum x}{\sum n} = \cfrac{79 + 91}{100 + 100} = \cfrac{170}{200} = 0.85 \]
Im Folgenden die Rechenschritte nochmal in R aufgedröselt zum besseren nachvollziehen. Wie auch schon im obigen Beispiel berechnen wir erst die beobachten Anteil \(O\) sowie die erwartenden Anteile \(E\). Dann nutzen wir die Formel des \(\mathcal{X}^2\)-Test um die \(\mathcal{X}^2\)-Statistik zu berechnen.
x <- c(76, 91)
n <- c(100, 100)
p <- sum(x)/sum(n)
O <- cbind(x, n - x)
E <- cbind(n * p, n * (1 - p))
sum((abs(O - E))^2/E)[1] 8.165487Auch hier vergleichen wir nochmal unser händisches Ergebnis mit dem Ergebnis der R Funktion prop.test(). Der Funktion übergeben wir dann einmal die beobachteten Anteile \(x\) sowie dann die jeweils Gesamtanzahlen \(n\). Wichtig ist hier, dass wir als 95% Konfidenzintervall die Differenz der beiden Wahrscheinlichkeiten erhalten.
prop.test(x = c(76, 91), n = c(100, 100), correct = FALSE)
2-sample test for equality of proportions without continuity
correction
data: c(76, 91) out of c(100, 100)
X-squared = 8.1655, df = 1, p-value = 0.00427
alternative hypothesis: two.sided
95 percent confidence interval:
-0.25076198 -0.04923802
sample estimates:
prop 1 prop 2
0.76 0.91 Wir wir sehen unterscheiden sich die beiden Anteil von \(76/100\) gleich 76% von \(91/100\) gleich 91%. Damit sollten wir den Roboter Picky nehmen, denn da werden prozentual weniger Erdbeeren zermanscht. Ob das jetzt gut oder schlecht ist 76% Erdbeeren zu zerstören ist aber wieder eine Frage, die die Statistik an dieser Stelle nicht beantworten kann.
35.6.3 Vergleich mehrerer Anteile \(p_1\), \(p_2\) bis \(p_k\)
Im Folgenden haben wir jetzt nicht mehr zwei Gruppen, die wir miteinander vergleichen wollen, sondern mehrere Behandlungen als Gruppen. Wir haben uns in einem Experiment die beschädigten Erdbeeren nach vier neue Erntearten, A bis D, angeschaut. Beide Vektoren können wir dann in die Funktion prop.test() stecken.
damaged <- c(A = 105, B = 100, C = 139, D = 96)
berries <- c(106, 113, 156, 102)
prop.test(damaged, berries)
4-sample test for equality of proportions without continuity
correction
data: damaged out of berries
X-squared = 11.747, df = 3, p-value = 0.008303
alternative hypothesis: two.sided
sample estimates:
prop 1 prop 2 prop 3 prop 4
0.9905660 0.8849558 0.8910256 0.9411765 Wir erhalten dann einen \(p\)-Wert für die Signifikanz von 0.0056 wieder. Was testen wir hier eigentlich? Unsere Nullhypothese ist, dass alle paarweisen Wahrscheinlichkeiten zwischen den Gruppen gleich sind.
\[ H_0: \; \mbox{Der Anteil der beschädigten Erdbeeren ist in den vier Gruppen ähnlich hoch} \]
\[ H_A: \; \mbox{Der Anteil der beschädigten Erdbeeren in mindestens einer der Behandlungen unterschiedlich ist.} \]
Wie wir sehen ist das nicht der Fall. Wir haben hier vier Wahrscheinlichkeiten vorliegen und mindestens zwei unterscheiden sich. Welche das sind, ist wieder die Frage. Hierzu nutzen wir dann gleich die Funktion pairwise.prop.test(). Wir immer geht die Ausgabe auch schöner und aufgeräumter.
prop.test(damaged, berries) |>
model_parameters()4-sample test for equality of proportions without continuity correction
Proportion | Chi2(3) | p
---------------------------------------------------
99.06% / 88.50% / 89.10% / 94.12% | 11.75 | 0.008
Alternative hypothesis: two.sidedWarum ist ein Test auf Anteile ein \(\mathcal{X}^2\)-Test? Hierfür brauchen wir noch die Informationen zu den nicht beschädigten Erdbeeren.
non_damaged <- berries - damaged
non_damaged A B C D
1 13 17 6 Nun können wir uns erstmal eine Tabelle bauen auf der wir dann den \(\mathcal{X}^2\)-Test und den prop.test() rechnen können. Der \(\mathcal{X}^2\)-Test ist nicht nur auf eine 2x2 Kreuztabelle beschränkt. Wir können in einem \(\mathcal{X}^2\)-Test auch andere \(n \times m\) Tabellen testen. Auf der anderen Seite ist der prop.test() auf eine \(n \times 2\) Tabelle beschränkt. Es müssen also immer zwei Spalten sein.
damaged_tab <- cbind(damaged, non_damaged) |>
as.table()
damaged_tab damaged non_damaged
A 105 1
B 100 13
C 139 17
D 96 6Wir erhalten die Tabelle 35.4 mit den beobachteten Werten sowie die Tabelle 35.5 mit den erwarteten Werten. Die Berechnung der erwarteten Werte kennen wir schon aus dem klassischen \(\mathcal{X}^2\)-Test. Hier machen wir die Berechnungen nur auf einer größeren Tabelle.
| Damaged | ||||
| Yes (1) | No (0) | |||
| damaged | non_damaged | berries | ||
| Group | A | \(105\) | \(1\) | \(\mathbf{106}\) |
| B | \(100\) | \(13\) | \(\mathbf{113}\) | |
| C | \(139\) | \(17\) | \(\mathbf{156}\) | |
| D | \(96\) | \(6\) | \(\mathbf{102}\) | |
| \(\mathbf{440}\) | \(\mathbf{37}\) | \(n = 477\) |
Jetzt können wir auch die Anteile der beschädigten Erdbeeren von allen Erdbeeren berechnen pro Ernteart berechnen.
\[ \begin{aligned} Pr(A) &= \cfrac{105}{106} = 0.9906 = 99.06\%\\ Pr(B) &= \cfrac{100}{113} = 0.8850 = 88.50\%\\ Pr(C) &= \cfrac{139}{156} = 0.8910 = 89.10\%\\ Pr(D) &= \cfrac{96}{102} = 0.9411 = 94.12\% \end{aligned} \]
| Damaged | ||||
| Yes (1) | No (0) | |||
| damaged | non_damaged | berries | ||
| Group | A | \(\cfrac{440 \cdot 106}{477} = 97.78\) | \(\cfrac{37 \cdot 106}{477} = 8.22\) | \(\mathbf{106}\) |
| B | \(\cfrac{440 \cdot 113}{477} = 104.23\) | \(\cfrac{37 \cdot 113}{477} = 8.77\) | \(\mathbf{113}\) | |
| C | \(\cfrac{440 \cdot 156}{477} = 143.90\) | \(\cfrac{37 \cdot 156}{477} = 12.10\) | \(\mathbf{156}\) | |
| D | \(\cfrac{440 \cdot 102}{477} = 94.09\) | \(\cfrac{37 \cdot 102}{477} = 7.91\) | \(\mathbf{102}\) | |
| \(\mathbf{440}\) | \(\mathbf{37}\) | \(n = 477\) |
Schauen wir uns nun an, ob es einen Unterschied zwischen den vier Erntearten A bis D für die Erdbeeren gibt. Einmal nutzen wir hierfür die Funktion chisq.test() und einmal die Funktion prop.test().
chisq.test(damaged_tab) |>
model_parameters()Pearson's Chi-squared test
Chi2(3) | p
---------------
11.75 | 0.008prop.test(damaged_tab) |>
model_parameters()4-sample test for equality of proportions without continuity correction
Proportion | Chi2(3) | p
---------------------------------------------------
99.06% / 88.50% / 89.10% / 94.12% | 11.75 | 0.008
Alternative hypothesis: two.sidedNachdem wir beide Funktionen gerechnet haben, sehen wir, dass beide Tests auf der \(\mathcal{X}^2\) Statistik basieren. Das macht ja auch Sinn, denn wir rechnen ja die Proportions indem wir die beobachteten Werte durch die Gesamtzahl berechnen. Hier haben wir im Prinzip die gleiche Idee wie schon in den beiden obigen Beispielen umgesetzt. Wir können daher den \(\mathcal{X}^2\)-Test auch einmal per Hand rechnen und kommen auf fast die gleiche \(\mathcal{X}^2\)-Statistik. Wir haben eine leichte andere Statistik, da wir hier mehr runden.
\[ \begin{aligned} \chi^2_{calc} &= \cfrac{(105 - 97.78)^2}{97.78} + \cfrac{(1 - 8.22)^2}{8.22} + \\ &\phantom{=}\;\; \cfrac{(100 - 104.23)^2}{104.23} + \cfrac{(13 - 8.77)^2}{8.77} + \\ &\phantom{=}\;\; \cfrac{(139 - 143.90)^2}{143.90} + \cfrac{(17 - 12.10)^2}{12.10} + \\ &\phantom{=}\;\; \cfrac{(96 - 94.09)^2}{94.09} + \cfrac{(6 - 7.91)^2}{7.91} = 11.74 \approx 10.04 \end{aligned} \]
Wir wissen nun, dass es mindestens einen paarweisen Unterschied zwischen den Wahrscheinlichkeiten für eine Beschädigung der Erdbeeren der vier Behandlungen gibt.
Führen wir den Test nochmal detaillierter mit der Funktion prop_test() aus dem R Paket {rstatix} durch. Es handelt sich hier um eine Alternative zu der Standardfunktion prop.test(). Das Ergebnis ist das Gleiche, aber die Aufarbeitung der Ausgabe ist anders und manchmal etwas besser weiterverarbeiteten.
Der Zugang ist etwas anders, deshalb bauen wir uns erstmal eine Tabelle mit den beschädigten und nicht beschädigten Erdbeeren. Dann benennen wir uns noch die Tabelle etwas um, damit haben wir dann einen besseren Überblick. Eigentlich unnötig, aber wir wollen uns hier ja auch mal mit der Programmierung beschäftigen.
damaged_tab <- rbind(damaged, non_damaged) |>
as.table()
dimnames(damaged_tab) <- list(
Infected = c("yes", "no"),
Groups = c("A", "B", "C", "D")
)
damaged_tab Groups
Infected A B C D
yes 105 100 139 96
no 1 13 17 6Dann können wir die Funktion prop_test() nutzen um den Test zu rechnen. Wir erhalten hier viele Informationen und müssen dann schauen, was wir dann brauchen. Dafür nutzen wir dann die Funktion select() und mutieren dann die Variablen in der Form, wie wir die Variablen haben wollen.
prop_test(damaged_tab, detailed = TRUE) |>
select(matches("estimate"), p) |>
mutate(p = pvalue(p)) |>
mutate_if(is.numeric, round, 2)# A tibble: 1 x 5
estimate1 estimate2 estimate3 estimate4 p
<dbl> <dbl> <dbl> <dbl> <chr>
1 0.99 0.88 0.89 0.94 0.008Wir erhalten auch hier das gleiche Ergebnis. War auch zu erwarten, denn im Kern sind die beiden Funktionen prop.test() und prop_test() gleich. Nun können wir uns einmal den paarweisen Vergleich anschauen. Wir wollen dann im Anschluß noch das Compact letter display für die Darstellung der paarweisen Vergleiche berechnen. Dafür brauchen wir dann noch folgende R Pakete.
Du kannst mehr über das Compact letter display in dem Kapitel Multiple Vergleiche oder Post-hoc Tests erfahren.
pacman::p_load(rcompanion, multcompView, multcomp)Die Standardfunktion für die paarweisen Vergleiche von Anteilen in R ist die Funktion pairwise.prop.test(). Wir wollen hier einmal die \(p\)-Werte nicht adjustieren, deshalb schreiben wir auch p.adjust.method = "none".
pairwise.prop.test(damaged, berries,
p.adjust.method = "none")
Pairwise comparisons using Pairwise comparison of proportions
data: damaged out of berries
A B C
B 0.0035 - -
C 0.0040 1.0000 -
D 0.1118 0.2264 0.2466
P value adjustment method: none Wir können auch die Funktion pairwise_prop_test() aus dem R Paket {rstatix} nutzen. Hier haben wir dann eine andere Ausgabe der Vergleiche. Manchmal ist diese Art der Ausgabe der Ergebnisse etwas übersichlicher. Wir nutzen dann noch die Funktion pvalue() um die \(p\)-Werte einmal besser zu formatieren.
pairwise_prop_test(damaged_tab,
p.adjust.method = "none") |>
mutate(p = pvalue(p.adj))# A tibble: 6 x 5
group1 group2 p p.adj p.adj.signif
<chr> <chr> <chr> <dbl> <chr>
1 A B 0.004 0.00354 **
2 A C 0.004 0.00401 **
3 B C >0.999 1 ns
4 A D 0.112 0.112 ns
5 B D 0.226 0.226 ns
6 C D 0.247 0.247 ns Dann benutzen wir noch die Funktion multcompLetters() um uns das Compact letter display wiedergeben zu lassen.
pairwise.prop.test(damaged, berries,
p.adjust.method = "none") |>
pluck("p.value") |>
fullPTable() |>
multcompLetters() |>
pluck("Letters") A B C D
"a" "b" "b" "ab" Hier sehen wir dann auch den Unterschied zwischen den beiden Funktionen. Wir können für die Funktion pairwise_prop_test() etwas einfacher das Compact letter display berechnen lassen. Wir müssen uns nur eine neue Spalte contrast mit den Vergleichen bauen.
pairwise_prop_test(damaged_tab,
p.adjust.method = "none") |>
mutate(contrast = str_c(group1, "-", group2)) |>
pull(p, contrast) |>
multcompLetters() A B C D
"a" "b" "b" "ab" Am Ende sehen wir, dass sich die Behandlung \(A\) von den Behandlungen \(B\) und \(C\) unterscheidet, da sich die Behandlungen nicht den gleichen Buchstaben teilen. Die Behandlung \(A\) unterscheidet sich aber nicht von der Behandlung \(D\). In dieser Art und Weise können wir dann alle Behandlungen durchgehen. Du kannst mehr über das Compact letter display in dem Kapitel Multiple Vergleiche oder Post-hoc Tests erfahren.