stat-tests-chi-test.qmd

```{r echo = FALSE}
pacman::p_load(tidyverse, readxl, knitr, kableExtra, Hmisc)
```

# Der $\mathcal{X}^2$-Test {#sec-chi-test}

*Letzte Änderung am `r format(fs::file_info("stat-tests-chi-test.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`*

> *"You're so square, baby, I don't care." --- Elvis Presley*

{{< video https://youtu.be/B6YNtxBsK3c >}}

::: callout-note
## 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.

::: {.callout-caution collapse="true"}
## 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](https://github.com/jkruppa/teaching/tree/main/Klausur) 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$.
:::

## Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

```{r echo = TRUE}
#| warning: false
#| message: false
pacman::p_load(tidyverse, magrittr, effectsize, rstatix,
               scales, parameters, conflicted)
conflicts_prefer(stats::chisq.test)
conflicts_prefer(stats::fisher.test)
```

Am Ende des Kapitels findest du nochmal den gesamten R Code in einem Rutsch zum selber durchführen oder aber kopieren.

## 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 @tbl-chi-square-obs 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$       |

: Eine 2x2 Tabelle als Beispiel für unterschiedliche Flohinfektionen bei Hunden und Katzen. Dargestellt sind die *beobachteten* Werte. {#tbl-chi-square-obs}

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.

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

## $\mathcal{X}^2$-Test theoretisch

In @tbl-chi-square-obs 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_{D} = \sum\cfrac{(O - E)^2}{E}
$$

mit

-   $O$ für die beobachteten Werte
-   $E$ für die nach den Randsummen zu erwartenden Werte

In @tbl-chi-square-exp 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$    |

: Eine 2x2 Tabelle als Beispiel für unterschiedliche Flohinfektionen bei Hunden und Katzen. Dargestellt sind die zu *erwartenden* Werte. {#tbl-chi-square-exp}

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.

```{=tex}
\begin{aligned} 
\chi^2_{D} &= \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_{D}$ von $1.25$ mit der Regel, dass wenn $\chi^2_{D} \geq \chi^2_{\alpha=5\%}$ die Nullhypothese abgelehnt werden kann. Mit einem $\chi^2_{\alpha=5\%} = `r round(qchisq(p = 0.05, df = 1, lower.tail = FALSE), 2)`$ 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_{D}$ 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!

## $\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.]{.aside}

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.

```{r}

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)

```

Als 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](https://easystats.github.io/effectsize/reference/phi.html) auf der Hilfeseite des R Paketes `{effectsize}`. Wir bleiben hier dann aber bei Cramers $V$.

```{r}
cramers_v(mat) 
```

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

```{r}
fisher.test(mat)
```

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.

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

### 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_{D} = \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.

```{r}
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)
```

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

```{r}
prop.test(x = 76, n = 100, p = 0.84, correct = FALSE) 
```

::: column-margin
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.
:::

### Vergleich zweier Anteile $p_1$ und $p_2$ {#sec-chisquare-prop-test}

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 @tbl-chi-square-exp-3 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$    |

: Eine 2x2 Tabelle für die zwei Ernteroboter und der beschädigten Erdbeeren. Dargestellt sind die *beobachteten* Werte. {#tbl-chi-square-exp-3}

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.

```{r}
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)
```

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

```{r}
prop.test(x = c(76, 91), n = c(100, 100), correct = FALSE) 
```

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.

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

```{r}
#| warning: false
#| message: false
damaged  <- c(A = 105, B = 100, C = 139, D = 96)
berries <- c(106, 113, 156, 102)
prop.test(damaged, berries)
```

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.

```{r}
#| warning: false
#| message: false
prop.test(damaged, berries) |> 
  model_parameters()
```

Warum ist ein Test auf Anteile ein $\mathcal{X}^2$-Test? Hierfür brauchen wir noch die Informationen zu den nicht beschädigten Erdbeeren.

```{r}
non_damaged <- berries - damaged 
non_damaged
```

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.

```{r}
#| warning: false
#| message: false
damaged_tab <- cbind(damaged, non_damaged) |> 
  as.table()
damaged_tab
```

Wir erhalten die @tbl-chi-square-exp-2 mit den beobachteten Werten sowie die @tbl-chi-square-exp-3 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$    |

: Eine 4x2 Tabelle für die vier Erntearten und der beschädigten Erdbeeren. Dargestellt sind die *beobachteten* Werte. {#tbl-chi-square-exp-2}

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$    |

: Eine 4x2 Tabelle für die vier Erntearten und der beschädigten Erdbeeren. Dargestellt sind die zu *erwartenden* Werte, die sich aus den Randsummen ergeben würden. {#tbl-chi-square-exp-3}

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()`.

```{r}
#| warning: false
#| message: false
chisq.test(damaged_tab) |> 
  model_parameters()
prop.test(damaged_tab) |> 
  model_parameters()
```

Nachdem 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_{D} &= \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.

```{r}
#| warning: false
#| message: false
damaged_tab <- rbind(damaged, non_damaged) |> 
  as.table()
dimnames(damaged_tab) <- list(
  Infected = c("yes", "no"),
  Groups = c("A", "B", "C", "D")
)
damaged_tab 
```

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

```{r}
#| warning: false
#| message: false
prop_test(damaged_tab, detailed = TRUE) |> 
  select(matches("estimate"), p) |> 
  mutate(p = pvalue(p)) |> 
  mutate_if(is.numeric, round, 2)
```

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

::: column-margin
Du kannst mehr über das *Compact letter display* in dem Kapitel [Multiple Vergleiche oder Post-hoc Tests](#sec-posthoc) erfahren.
:::

```{r}
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"`.

```{r}
#| warning: false
#| message: false
pairwise.prop.test(damaged, berries, 
                   p.adjust.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.

```{r}
#| warning: false
#| message: false
pairwise_prop_test(damaged_tab, 
                   p.adjust.method = "none") |> 
  mutate(p = pvalue(p.adj))
```

Dann benutzen wir noch die Funktion `multcompLetters()` um uns das *Compact letter display* wiedergeben zu lassen.

```{r}
#| warning: false
#| message: false
pairwise.prop.test(damaged, berries, 
                   p.adjust.method = "none") |> 
  pluck("p.value") |> 
  fullPTable() |> 
  multcompLetters() |> 
  pluck("Letters")
```

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.

```{r}
#| warning: false
#| message: false

pairwise_prop_test(damaged_tab, 
                   p.adjust.method = "none") |> 
  mutate(contrast = str_c(group1, "-", group2)) |> 
  pull(p, contrast) |> 
  multcompLetters() 
```

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](#sec-posthoc) erfahren.