stat-tests-effect.qmd

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

# Der Effektschätzer {#sec-effect}

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

> *"Alles, was wir hören, ist eine Meinung, keine Tatsache. Alles, was wir sehen, ist eine Perspektive, nicht die Wahrheit." --- Marcus Aurelius, Meditationen*

Der Effektschätzer. Ein seltsames Kapitel, denn ich tue mich sehr schwer, dieses Kapitel irgendwo in die Linearität dieses Buches hier einzuordnen. Deshalb ist dieses Kapitel eigentlich immer an der falschen Stelle. Entweder hast du schon die statistischen Tests gelesen und du wüsstest gerne was die Effektschätzer sind *oder* du suchst hier nochmal die Beschreibung der Effektschätzer zum Beispiel aus der multiple Regression heraus. Also steht jetzt dieses Kapitel hier im Raum und du musst schauen, was du wirklich brauchst. Oder ob du dieses Kapitel erst überspringst und dann später nochmal hier liest. Eine wunderbare Übersicht über den Begriff *Effektschätzer* liefert das englische Buch [Doing Meta-Analysis with R: A Hands-On Guide](https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/effects.html). Der Effektschätzer wird dabei auch gerne Theta $\boldsymbol{\theta}$ genannt. Da wir dann aber später noch mit anderen Konzepten in die Quere kommen, nutze ich das etwas intuitivere Delta $\boldsymbol{\Delta}$.

::: {.callout-note collapse="true"}
## Zerforschen: Science Cops

{{< include zerforschen/zerforschen-science-cops.qmd >}}
:::

Wenn wir einen der vielen Effektschätzer berechnen wollen, dann nutzen wir dafür die Effektschätzer aus dem [R Paket `{effectsize}`](https://easystats.github.io/effectsize/index.html). Das R Paket `{effectsize}` liefert Effektschätzer für fast alle statistischen Gelegenheiten. Wir werden hier wie immer nur den groben Überblick abdecken. Vermutlich wird das Kapitel dann noch Anwachsen. Streng genommen gehört das @sec-test-diag zu den diagnostischen Tests auf einer 2x2 Kreuztabelle auch irgendwie zu Effektschätzern. Wenn du Spezifität und Sensitivität suchst bist du in dem Kapitel zu diagnostischen Tests richtig.

Wir unterscheiden hier erstmal grob in vier Arten von Effektschätzern:

1)  Effektschätzer, die einen *Mittelwertsunterschied* zwischen Behandlungen beschreiben.
2)  Effektschätzer, die einen *Anteilsunterschied* zwischen Behandlungen beschreiben.
3)  Effektschätzer, die einen *Wirkungsgrad* einer Behandlung zu einer Kontrolle beschreiben.
4)  Effektschätzer, die eine *Übereinstimmung* zwischen Bewertern beschreibt.

Daneben gibt es wie noch die Korrelation wie in @sec-lin-reg-corr beschrieben. Die Korrelation wollen wir aber in diesem Kapitel nicht vorgreifen bzw. wiederholen.

Am Ende muss du immer den Effekt im Kontext der Fragestellung bzw. des Outcomes $y$ bewerten. Der numerische Unterschied von $0.1$ cm kann in einem Kontext viel sein. Das Wachstum von Bakterienkolonien kann ein Unterschied von $0.1$ cm viel sein. Oder aber sehr wenig, wenn wir uns das Wachstum von Bambus pro Tag anschauen. Hier bist du gefragt, den Effekt in den Kontext richtig einzuordnen. Ebenso stellt sich die Frage, ob ein Unterschied von 6% viel oder wenig ist.

::: callout-tip
## Effektschätzer zwischen Behandlungen

Wenn wir uns einen Unterschied eines **Mittelwerts** anschauen, dann haben wir *keinen* Effekt vorliegen, wenn der Mittlwertsunterschied $\Delta$ zwischen der Gruppe $A$ und der Gruppe $B$ gleich 0 ist. Die Nullhypothese gilt. Beide Gruppen $A$ und $B$ haben den gleichen Mittelwert.

$$
\Delta_{A-B} = A - B = 0
$$

Wenn wir uns einen Unterschied eines **Anteils** anschauen, dann haben wir *keinen* Effekt vorliegen, wenn der Anteilsunterschied $\Delta$ zwischen der Gruppe $A$ und der Gruppe $B$ gleich 1 ist. Die Nullhypothese gilt. Beide Gruppen $A$ und $B$ haben den gleichen Anteil.

$$
\Delta_{A/B} = \cfrac{A}{B} = 1
$$

Dieses Wissen brauchen wir um die Signifikanzschwelle bei einem 95% Konfidenzintervall richtig zu setzen und interpretieren zu können. Siehe dazu auch nochmal das @sec-ki.
:::

## Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

```{r echo = TRUE}
#| message: false
pacman::p_load(tidyverse, magrittr, see, scales,
               effectsize, parameters, broom, readxl,
               emmeans, vcd, irr, conflicted)
conflicts_prefer(dplyr::filter)
conflicts_prefer(effectsize::oddsratio)
conflicts_prefer(magrittr::set_names)

```

An der Seite des Kapitels findest du den Link *Quellcode anzeigen*, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.

## Unterschied zweier Mittelwerte

Wir berechnen zwei Mittelwerte $\bar{y}_1$ und $\bar{y}_2$. Wenn wir wissen wollen wie groß der Effekt zwischen den beiden Mittelwerten ist, dann bilden wir die Differenz. Wir berechnen das $\Delta_{y_1-y_2}$ für $y_1$ und $y_2$ indem wir die beiden Mittelwerte voneinander abziehen.

$$
\Delta_{y_1-y_2} = \bar{y}_1 - \bar{y}_2
$$

Wenn es keinen Unterschied zwischen den beiden Mittelwerten $\bar{y}_1$ und $\bar{y}_2$ gibt, dann ist die Differenz $\Delta_{y_1-y_2} = \bar{y}_1 - \bar{y}_2$ gleich 0. Wir sagen, die Nullhypothese *vermutlich* gilt, wenn die Differenz klein ist. Warum schreiben wir hier vermutlich? Ein statistischer Test ist eine Funktion von $\Delta$, $s$ und $n$. Wir können auch mit kleinem $\Delta$ die Nullhypothese ablehnen, wenn $s$ und $n$ eine passende Teststatistik generieren. Siehe dazu auch das @sec-delta-n-s. Was wir besser annehmen können ist, dass die Relevanz klein ist. Effekt mit einem geringen Mittelwertsunterschied sind meistens nicht relevant. Aber diese Einschätzung hängt stark von der Fragestellung ab.

$$
H_0: \Delta_{y_1-y_2} = \bar{y}_1 - \bar{y}_2 = 0
$$

In @tbl-dog-cat-small-delta ist nochmal ein sehr simples Datenbeispiel gegeben an dem wir den Zusammenhang nochmal nachvollziehen wollen.

```{r}
#| echo: false
#| label: tbl-dog-cat-small-delta
#| tbl-cap: Beispiel für die Berechnung von einem Mittelwertseffekt an der Sprunglänge [cm] von Hunde und Katzenflöhen.
data_tbl <- tibble(animal = gl(2, 4, labels = c("cat", "dog")),
                   jump_length = c(8.0, 7.9, 8.3, 9.1, 8.0, 7.8, 9.2, 7.7)) 

data_tbl %>% 
  kable(align = "c", "pipe")

mean_tbl <- data_tbl %>% 
  group_by(animal) %>% 
  summarise(mean_animal = round(mean(jump_length), 1),
            var_animal = round(var(jump_length), 1))

mean_cat <- mean_tbl$mean_animal[1]
mean_dog <- mean_tbl$mean_animal[2]

```

Nehmen wir an, wir berechnen für die Sprungweite \[cm\] der Hundeflöhe einen Mittelwert von $\bar{y}_{dog} = `r mean_dog`$ und für die Sprungweite \[cm\] der Katzenflöhe einen Mittelwert von $\bar{y}_{cat} =`r mean_cat`$. Wie große ist nun der Effekt? Oder anders gesprochen, welchen Unterschied in der Sprungweite macht es aus ein Hund oder eine Katze zu sein? Was ist also der Effekt von `animal`? Wir rechnen $\bar{y}_{dog} - \bar{y}_{cat} = 8.2 - 8.3 = -0.1$. Zum einen wissen wir jetzt "die Richtung". Da wir ein Minus vor dem Mittelwertsunterschied haben, müssen die Katzenflöhe weiter springen als die Hundeflöhe, nämlich 0.1 cm. Dennoch ist der Effekt sehr klein.

### Cohen's d

Da der Mittlwertsunterschied *alleine* nnur eine eingeschränkte Aussage über den Effekt erlaubt, gibt es noch Effektschätzer, die den Mittelwertsunterschied $\Delta_{y_1-y_2}$ mit der Streuung $s^2$ sowie der Fallzahl zusammenbringt. Der bekannteste Effektschätzer für einen Mittelwertsunterschied bei großer Fallzahl mit mehr als 20 Beobachtungen ist Cohen's d. Wir können Cohen's d wie folgt berechnen.

$$
|d| = \cfrac{\bar{y}_1-\bar{y}_2}{\sqrt{\cfrac{s_1^2+s_2^2}{2}}}
$$

Wenn wir die berechneten Mittelwerte und die Varianz der beiden Gruppen in die Formel einsetzten ergibt sich ein absolutes Cohen's d von 0.24 für den Gruppenvergleich.

$$
|d| = \cfrac{8.2 - 8.3}{\sqrt{(0.5^2+0.3^2) /2}} = \cfrac{-0.1}{0.41} = \lvert-0.24\rvert
$$

Was denn nun Cohen's d *exakt* aussagt, kann niemand sagen. Aber wir haben einen Wust an möglichen Grenzen. Hier soll die Grenzen von Cohen (1988) einmal angegeben werden. Cohen (1988) hat in seiner Arbeit folgende Grenzen in @tbl-cohen-d für die Interpretation von $d$ vorgeschlagen. Mehr Informationen zu Cohen's d gibt es auf der Hilfeseite von `effectsize`: [Interpret standardized differences](https://easystats.github.io/effectsize/reference/interpret_cohens_d.html).

|     Cohen's d      | Interpretation des Effekts |
|:------------------:|:--------------------------:|
|     $d < 0.2$      |         Sehr klein         |
| $0.2 \leq d < 0.5$ |           Klein            |
| $0.5 \leq d < 0.8$ |           Mittel           |
|    $d \geq 0.8$    |           Stark            |

: Interpretation der Effektstärke nach Cohen (1988). {#tbl-cohen-d}

Wir können auch über die Funktion `cohens_d()` Cohen's d einfach in R berechnen. Die Funktion `cohens_d()` akzeptiert die Formelschreibweise. Die 95% Konfidenzintervalle sind mit Vorsicht zu interpretieren. Denn die Nullhypothese ist hier nicht so klar formuliert. Wir lassen also die 95% Konfidenzintervalle erstmal hier so stehen.

```{r}
cohens_d(jump_length ~ animal, data = data_tbl, pooled_sd = TRUE)
```

Dankenswerterweise gibt es noch die Funktion `interpret_cohens_d`, die es uns erlaubt auszusuchen nach welche Literturquelle wir den Wert von Cohen's d interpretieren wollen. Ob dieser Effekt relevant zur Fragestellung ist musst du selber entscheiden.

```{r}
interpret_cohens_d(0.24, rules = "cohen1988")
```

### Hedges' g

Soweit haben wir uns mit sehr großen Fallzahlen beschäftigt. Cohen's d ist dafür auch sehr gut geeigent und wenn wir mehr als 20 Beobachtungen haben, können wir Cohen's d auch gut anwenden. Wenn wir weniger Fallzahl vorliegen haben, dann können wir Hedges' g nutzen. Hedges' g bietet eine Verzerrungskorrektur für kleine Stichprobengrößen ($N < 20$) sowie die Möglichkeit auch für *unbalanzierte* Gruppengrößen einen Effektschätzeer zu berechnen. Die Formel sieht mit dem Korrekturterm recht mächtig aus.

$$
g = \cfrac{\bar{y}_1 - \bar{y}_2}{s^*} \cdot \left(\cfrac{N-3}{N-2.25}\right) \cdot \sqrt{\cfrac{N-2}{N}}
$$

mit

$$
s^* = \sqrt{\cfrac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}}
$$

Wir können aber einfach die Mittelwerte und die Varianzen aus unserem beispiel einsetzen. Da unsere beiden Gruppen gleich groß sind $n_1 = n_2$ und damit ein balanziertes Design vorliegt, sind Cohen's d und Hedges' g numerisch gleich. Wir können dann noch für die geringe Fallzahl korrigieren und erhalten ein händisches $g = 0.18$.

$$
g = \cfrac{8.2 - 8.3}{0.41} \cdot \left(\cfrac{8-3}{8-2.25}\right) \cdot \sqrt{\cfrac{8-2}{8}} = \lvert-0.24\rvert \cdot 0.87 \cdot 0.87 \approx 0.18
$$

mit

$$
s^* = \sqrt{\cfrac{(4-1)\cdot0.5^2 + (4-1)\cdot0.3^2}{4+4-2}} = \sqrt{\cfrac{0.75 + 0.27}{6}} = 0.41
$$

In R gibt es die Funktion `hedges_g()` die uns erlaubt in der Formelschreibweise direkt Hedges' g zu berechnen. Wir sehen hier eine Abweichung von unserer händischen Rechnung. Das ist aber in soweit nicht ungewöhnlich, da es noch eine Menge Varianten der Anpassung für die geringe Fallzahl gibt. In der Anwendung nutzen wir die Funktion aus dem Paket `{effectsize}` wie hier durchgeführt.

Wir ignorieren wie auch bei Cohen's d das 95% Konfidenzintervall, da die Interpretation ohne die Nullhypothese nicht möglich ist. Die Nullhypothese ist in diesem Fall komplexer. Wir lassen daher das 95% Konfidenzintervall erstmal einfach hier so stehen.

```{r}
hedges_g(jump_length ~ animal, data = data_tbl, pooled_sd = TRUE)
```

Auch für Hedges' g gibt es die Möglichkeit sich über die Funktion `interpret_hedges_g()` den Wert von $g=0.21$ interpretieren zu lassen. Nach Sawilowsky (2009) haben wir hier einen kleinen Effekt vorliegen. Ob dieser Effekt relevant zur Fragestellung ist musst du selber entscheiden.

```{r}
interpret_hedges_g(0.21, rules = "sawilowsky2009")
```

[Die Hilfeseite zu dem Paket `effectsize`](https://easystats.github.io/effectsize/reference/interpret_cohens_d.html) bietet eine Liste an möglichen Referenzen für die Wahl der Interpretation der Effektstärke. Du musst dann im Zweifel schauen, welche der Quellen und damit Grenzen du nutzen willst.

## Unterschied zweier Anteile

Neben den Unterschied zweier Mittelwerte ist auch häufig das Interesse an dem Unterschied zwischen zwei Anteilen. Nun unterscheiden wir zwischen Wahrscheinlichkeiten und Chancen. Beide Maßzahlen, die Wahrscheinlichkeit wie auch die Chance, beschreiben einen Anteil. Hier tritt häufig Verwirrung auf, daher hier zuerst ein Beispiel.

Wir behandelt $n = 65$ Hunde mit dem Antiflohmittel *FleaEx*. Um die Wirkung von *FleaEx* auch bestimmen zu können haben wir uns zwei Gruppen von Hunden ausgesucht. Wir haben Hunde, die mit Flöhe infiziert sind und Hunde, die nicht mit Flöhen infiziert sind. Wir schauen nun in wie weit *FleaEx* gegen Flöhe hilft im Vergleich zu einer Kontrolle.

|              |           |                   |                   |
|:------------:|:---------:|:-----------------:|:-----------------:|
|              |           |     **Group**     |                   |
|              |           |     *FleaEx*      |     *Control*     |
| **Infected** | *Yes (1)* | $18_{\;\Large a}$ | $23_{\;\Large b}$ |
|              | *No (0)*  | $14_{\;\Large c}$ | $10_{\;\Large d}$ |

: Eine 2x2 Tabelle als Beispiel für unterschiedliche Flohinfektionen bei nach einer Behandlung mit *FleaEx* für die Berechnung von Effektschätzern eines Anteils. {#tbl-2x2-ratio-delta}

Aus der @tbl-2x2-ratio-delta können wir entnehmen, dass 18 behandelte Hunde mit Flöhen infiziert sind und 14 Hunde keine Infektion aufweisen. Bei den Hunden aus der Kontrolle haben wir 23 infizierte und 10 gesunde Tiere beobachtet.

Wir können nun zwei Arten von Anteilen berechnen um zu beschreiben, wie sich der Anteil an infizierten Hunden verhält. Das bekanntere ist die Frequenz oder Wahrscheinlichkeit oder Risk Ratio ($RR$). Das andere ist das Chancenverhältnis oder Odds Ratio ($OR$). Beide kommen in der Statistik vor und sind unterschiedlich zu interpretieren. Es gibt verschiedene Typen von klinischen Studien, also Untersuchungen an Menschen. Einige Studien liefern nur $OR$ wieder andere Studientypen liefern $RR$. @george2020s liefert eine gute Übersicht über *What's the risk: differentiating risk ratios, odds ratios, and hazard ratios?*

Um die die Odds Ratio und die Risk Ratios auch in R berechnen zu können müssen wir einmal die 2x2 Kreuzabelle in R nachbauen. Wir nutzen dafür die Funktion `matrix()` und müssen schauen, dass die Zahlen in der 2x2 Kreuztabelle in R dann auch so sind, wie in der Datentabelle. Das ist jetzt ein schöner Codeblock, ist aber dafür da um sicherzustellen, dass wir die Zahlen richtig eintragen.

```{r}
cross_mat <- matrix(c(18, 23, 14, 10),
  nrow = 2, byrow = TRUE,
  dimnames = list(
    Infected = c("Yes", "No"),
    Group = c("FleaEx", "Control")
  )
)

cross_mat
```

Später werden wir das $OR$ und $RR$ wieder treffen. Das $OR$ kommt in der logistsichen Regression als Effektschätzer vor. Wir nutzen das $RR$ als Effektschätzer in der Poissonregression.

### Wahrscheinlichkeitsverhältnis oder Risk Ratio (RR)

Wir berechnen wir nun das Wahrscheinlichkeitsverhältnis oder Risk Ratio (RR)? Das Risk Ratio ist das Verhältnis von den infizierten Hunden in der Behandlung ($a$) zu allen infizierten Hunden ($a+c$) zu dem Verhältnis der gesunden Hunde in der Behandlung ($b$) zu allen gesunden Hunden ($b+d$). Das klingt jetzt etwas wirr, deshalb helfen manchmal wirklich Formeln, den Zusammenhang besser zu verstehen. Damit ist $Pr(\mbox{FleaEx}|\mbox{infected})$ die Wahrscheinlichkeit infiziert zu sein, wenn der Hund mit *FleaEx* behandelt wurde.

$$
Pr(\mbox{FleaEx}|\mbox{infected}) = \cfrac{a}{a+c} = \cfrac{18}{18+14} \approx 0.56
$$

Die Wahrscheinlichkeit $Pr(\mbox{Control}|\mbox{infected})$ ist dann die Wahrscheinlichkeit infiziert zu sein, wenn der Hund mit in der Kontrolle war.\]{.aside}

$$
Pr(\mbox{Control}|\mbox{infected}) = \cfrac{b}{b+d} = \cfrac{23}{23 + 10} \approx 0.70
$$

Das Risk Ratio ist mehr oder minder das Verhältnis von der beiden Spalten der @tbl-2x2-ratio-delta für die Behandlung. Wir erhalten also ein $RR$ von $0.76$. Damit mindert die Gabe von *FleaEx* die Wahrscheinlichkeit sich mit Flöhen zu infizieren.

$$
\Delta_{y_1/y_2} = RR = \cfrac{Pr(\mbox{FleaEx}|\mbox{infected})}{Pr(\mbox{Control}|\mbox{infected})} =  \cfrac{0.56}{0.70} \approx 0.80 
$$

Wir überprüfen kurz mit der Funktion `riskratio()` ob wir richtig gerechnet haben. Das 95% Konfidenzintervall können wir interpretieren, dafür brauchen wir aber noch einmal eine Idee was "kein Effekt" bei einem Risk Ratio heist.

```{r}
riskratio(cross_mat)
```

Wann liegt nun kein Effekt bei einem Anteil wie dem RR vor? Wenn der Anteil in der einen Gruppe genauso groß ist wie der Anteil der anderen Gruppe.

$$
H_0: RR = \cfrac{Pr(\mbox{dog}|\mbox{infected})}{Pr(\mbox{cat}|\mbox{infected})} = 1
$$

Wir interpretieren das $RR$ nun wie folgt. Unter der Annahme, dass ein kausaler Effekt zwischen der Behandlung und dem Outcome besteht, können die Werte des relativen Risikos auf folgende Art und Weise interpretiert werden:

-   $RR = 1$ bedeutet, dass die Behandlung keinen Einfluss auf das Outcome hat
-   $RR < 1$ bedeutet, dass das Risiko für das Outcome durch die Behandlung verringert wird, was ein "Schutzfaktor" ist
-   $RR > 1$ bedeutet, dass das Risiko für das Outcome durch die Behandlung erhöht wird, was ein "Risikofaktor" ist.

Das heist in unserem Fall, dass wir mit einem RR von $0.80$ eine protektive Behandlung vorliegen haben. Die Gabe von *FleaEx* reduziert das Risiko mit Flöhen infiziert zu werden. Durch das 95% Konfidenzintervall wissen wir auch, dass das $RR$ nicht signifikant ist, da die 1 im 95% Konfidenzintervall enthalten ist.

### Chancenverhältnis oder Odds Ratio (OR)

Neben dem Risk Ratio gibt es noch das Odds Ratio. Das Odds Ratio ist ein Chancenverhältnis. Wenn der Mensch an sich schon Probleme hat für sich Wahrscheinlichkeiten richtig einzuordnen, scheitert man allgemein an der Chance vollkommen. Dennoch ist das Odds Ratio eine gute Maßzahl um abzuschätzen wie die Chancen stehen, einen infizierten Hund vorzufinden, wenn der Hund behandelt wurde.

Schauen wir uns einmal die Formeln an. Im Gegensatz zum Risk Ratio, welches die Spalten miteinander vergleicht, vergleicht das Odds Ratio die Zeilen. Als erstes berechnen wir die Chance unter der Gabe von *FleaEx* infiziert zu sein wie folgt.

$$
Odds(\mbox{FleaEx}|\mbox{infected}) = a:b = 18:23 = \cfrac{18}{23} = 0.78
$$

Dann berechnen wir die Chance in der Kontrollgruppe infiziert zu sein wie folgt.

$$
Odds(\mbox{Control}|\mbox{infected}) = c:d = 14:10 = \cfrac{14}{10} \approx 1.40 
$$

Abschließend bilden wir das Chancenverhältnis der Chance unter der Gabe von *FleaEx* infiziert zu sein zu der Chance in der Kontrollgruppe infiziert zu sein. Es ergbit sich das Odds Ratio wie folgt.

$$
\Delta_{y_1/y_2} = OR =  \cfrac{Odds(\mbox{Flea}|\mbox{infected})}{Odds(\mbox{Control}|\mbox{infected})} = \cfrac{a \cdot d}{b \cdot c} = \cfrac{0.78}{1.40} \approx 0.56
$$

Wir überprüfen kurz mit der Funktion `oddsratio()` ob wir richtig gerechnet haben. Das 95% Konfidenzintervall können wir interpretieren, dafür brauchen wir aber noch einmal eine Idee was "kein Effekt" bei einem Odds Ratio heist.

```{r}
oddsratio(cross_mat)
```

Wann liegt nun kein Effekt bei einem Anteil wie dem OR vor? Wenn der Anteil in der einen Gruppe genauso groß ist wie der Anteil der anderen Gruppe.

$$
H_0: OR =  \cfrac{Odds(\mbox{dog}|\mbox{infected})}{Odds(\mbox{cat}|\mbox{infected})} = 1
$$

Wir interpretieren das $OR$ nun wie folgt. Unter der Annahme, dass ein kausaler Effekt zwischen der Behandlung und dem Outcome besteht, können die Werte des Odds Ratio auf folgende Art und Weise interpretiert werden:

-   $OR = 1$ bedeutet, dass die Behandlung keinen Einfluss auf das Outcome hat
-   $OR < 1$ bedeutet, dass sich die Chance das Outcome zu bekommen durch die Behandlung verringert wird, was ein "Schutzfaktor" ist
-   $OR > 1$ bedeutet, dass sich die Chance das Outcome zu bekommen durch die Behandlung erhöht wird, was ein "Risikofaktor" ist.

Das heist in unserem Fall, dass wir mit einem OR von $0.56$ eine protektive Behandlung vorliegen haben. Die Gabe von *FleaEx* reduziert die Chance mit Flöhen infiziert zu werden. Durch das 95% Konfidenzintervall wissen wir auch, dass das $OR$ nicht signifikant ist, da die 1 im 95% Konfidenzintervall enthalten ist.

### Odds Ratio (OR) zu Risk Ratio (RR)

Wenn wir das OR berechnet haben, wollen wir eventuell das $OR$ in dem Sinne eines Riskoverhältnisses berichten. Leider ist es nun so, dass wir das nicht einfach mit einem $OR$ machen können. Ein $OR$ von 3.5 ist ein großes Chancenverhältnis. Aber ist es auch 3.5-mal so wahrscheinlich? Nein so einfach können wir das OR nicht interpretieren. Wir können aber das $OR$ in das $RR$ umrechnen. Dafür brauchen wir aber das $p_0$. Dabei ist das $p_0$ das Basisrisiko also die Wahrscheinlichkeit des Ereignisses ohne die Intervention. Wenn wir nichts tun würden, wie wahrscheinlich wäre dann das Auftreten des Ereignisses? Es ergibt sich dann die folgende Formel für die Umrechnung des $OR$ in das $RR$. @grant2014converting gibt nochmal eine wissenschaftliche Diskussion des Themas zur Konvertierung von OR zu RR.

$$
RR = \cfrac{OR}{(1 - p_0 + (p_0 \cdot OR))}
$$

Schauen wir uns das einmal in einem Beispiel an. Wir nutzen für die Umrechnung die Funktion `oddsratio_to_riskratio()` aus dem R Paket `{effectsize}`. Wenn wir ein $OR$ von 3.5 haben, so hängt das $RR$ von dem Basisriskio ab. Wenn das Basisirisko für die Erkrankung ohne die Behandlung sehr hoch ist mit $p_0 = 0.85$, dann ist das $RR$ sehr klein.

```{r}
OR <- 3.5
baserate <- 0.85

oddsratio_to_riskratio(OR, baserate) %>% round(2)
```

Auf der anderen Seite nähert sich das $OR$ dem $RR$ an, wenn das Basisriskio für die Erkrankung mit $p_0 = 0.04$ sehr klein ist.

```{r}
OR <- 3.5
baserate <- 0.04

oddsratio_to_riskratio(OR, baserate) %>% round(2)
```

Weil wir natürlich das Basisrisiko nur abschätzen können, verbleibt hier eine gewisse Unsicherheit, wie das $RR$ zu einem gegebenen $OR$ aussieht.

### Foldchange

Der Foldchange (deu. *Zunahmenänderung* ungebräuchlich, abk. *FC*) ist ein statistisches Maß welches uns beschreibt, wie sehr sich eine Menge von einem Anfangs- zu einem Endwert verändert. Dabei ist es wichtig, dass wir es hier gleich mit einem Anteil zu tun haben. Beispielsweise entspricht ein Anfangswert von 30 und ein Endwert von 60 einer 2-fachen Änderung oder eben eines *two-fold change*. Wir wollen dabei natürlich wissen, wie sich der Wert von Ende zu Beginn geändert hat.

$$
FC = \cfrac{Endwert}{Anfangswert} = \cfrac{60}{30} = 2
$$

Eine Verringerung von 80 auf 20 wäre ein Foldchange von -0.75, während eine Veränderung von 20 auf 80 ein Foldchange von 3 wäre. Damit sind wir natürlich sehr unsymmetrisch um die 1.

$$
FC = \cfrac{Endwert}{Anfangswert} = \cfrac{20}{80} = 0.75
$$

Alle Foldchanges mit einer Verringerung pressen sich in den Zahlenraum $[0,1]$ während alle Erhöhungen sich im Zahlenraum $[1, \infty]$ bewegen. Um dies zu umgehen, verwenden wir $log2$ für den Ausdruck der Log Foldchange. Damit du das besser verstehst, dann die Sachlage einmal als ein Beispiel. Wenn du einen FC von 2 vorliegen hast und dann $log_2(2)$ berechnest, dann erhälst du einen $log_2FC$ von 1.

```{r}
FC <- c(0.1, 0.5, 0.75, 1, 2, 3, 4)
log2(FC) |> 
  round(2) |> 
  set_names(FC)
```

Wenn du einen log2-transformierten FC-Wert wieder zurück haben willst, dann nutze einfach die 2-er Potenz mit $2^{logFC}$. Auch hier einmal das schnelle Beispiel.

```{r}
2^log2(FC)
```

Wenn dir also einmal ein Foldchange über den Weg läuft, dann geht es eigentlich nur darum eine Veränderung gleichmäßig um die Eins darstellen zu können.

## Wirkungsgrad von Pflanzenschutzmitteln {#sec-effect-wirkung}

Neben den klassischen Effektmaßzahlen, die sich aus einem Mittelwert oder einem Anteil direkt berechnen, gibt es noch andere Effektmaße. Einer dieser Effektmaße ist der Wirkungsgrad für zum Beispiel ein Pflanzenschutzmittel. Wir können hier aber auch weiter denken und uns überlegen in wie weit wir eine Population von Schaderregern durch eine Behandlung reduzieren können. Unabdingbar ist in diesem Fall eine positive Kontrolle in der nichts gemacht wird sondern nur der normale Befall gemessen wird. Merke also, ohne eine Kontrolle geht bei der Bestimmung des Wirkungsgrades gar nichts.

::: callout-tip
## Übersicht über Wirkungsgrade

-   Die Seite [LdP Line by Ehab Bakr](https://www.ehabsoft.com/ldpline/onlinecontrol.htm) liefert nochmal die Formeln und Berechungsmöglichkeiten für andere Wirkungsgrade mit einem praktischen Online-Rechner. Leider ist die Seite sehr alt und die Sicherheitsbestimmungen lassen moderne Browser zweifeln.
-   Der [Vortrag von Andreas Büchse](http://www.biometrische-gesellschaft.de/fileadmin/AG_Daten/Landwirtschaft/tagungsberichte_pdf/buechse.pdf) *"Nutzung der Wirkungsgradberechnung nach Abbott und Henderson-Tilton in der angewandten Agrarforschung"* gibt auch nochmal weitere Einblicke und ist die Grundlage für die Beispiele hier.
:::

In der @tbl-effect-wirk sehen wir nochmal eine Übersicht über vier mögliche Arten den Wirkungsgrad zu berechnen. Je nachdem was wir vorliegen haben, entweder einen Befall mit Zähldaten oder aber eine Sterblichkeit bzw. Mortalität, müssen wir eine andere Methode wählen. Auch müssen wir unterscheiden, ob wir eine homogene Population vorliegen haben oder eher eine heterogene Population. Ich kann hier aber schon mal etwas abkürzen, meistens können wir nur den Wirkungsgrad nach @abbott1925method berechnen, da unser experimentelles Design nicht mehr hergibt. Der Wirkungsgrad nach Abbott hat nämlich nur die räumliche Komponente und nicht die zeitliche. So lässt sich Abbott häufig einfacher durchführen und ausrechnen. Für eine zeitliche Komponente wie bei Henderson-Tilton müssten wir dann den Befall wiederholt messen.

|                                        |                                                     |                                                                        |
|----------------------------------------|:---------------------------------------------------:|:----------------------------------------------------------------------:|
|                                        |               **Homogene Population**               |          **Heterogene Population (Vorbefall differenziert)**           |
| **Befall oder lebende Individuen**     | [Abbott](#sec-effect-abbott) \[*räumlicher Bezug*\] | [Henderson-Tilton](#sec-effect-ht) \[*zeitlich und räumlicher Bezug*\] |
| **Sterblichkeit oder tote Individuen** |         [Schneider-Orelli](#sec-effect-so)          |                     [Sun-Shepard](#sec-effect-ss)                      |

: Übersicht der möglichen Wirkungsgrade nach Population und gemessen Werten an den Beobachtungen. {#tbl-effect-wirk}

Im Folgenden gehen wir einmal die vier Wirkungsgrade durch. Neben Abbott stelle ich dann hier nochmal den Wirkungsgrad nach @henderson1955tests vor. Die anderen Wirkungsgrade nur als Formeln, da mir aktuell noch keine Anwendungen hier über den Weg gelaufen sind. Der Hauptfokus liegt dann aber am Ende auf den Wirkungsgrad nach Abbott. Den Wirkungsgrad brauchen wir eben dann doch am meisten, da wir nach Abbott nur den gezählten Befall nach der Behandlung bestimmen müssen. Der Wirkungsgrad nach Abbott entspricht meistens eher der experimentellen Wirklichkeit. Das soll dich natürlich nicht davon abschrecken auch die anderen Wirkunsggrade zu nutzen, aber dafür musst du meistens mehr Aufwand betreiben.

### Wirkungsgrad nach @henderson1955tests {#sec-effect-ht}

Der Wirkungsgrad nach Henderson-Tilton hate den großen Vorteil, dass wir auch eine zzeitliche Komponente mit berücksichtigen können. Wie verhält sich der BEfall der Kontrolle und der Behandlung über die Zeit des Versuches? Wir müssen dafür den Befall vor (eng. *before*, abk. *b*) und den Befall nach (eng. *after*, abk. *a*) der Applikation bestimmen. Dann erhalten wir folgende Formel, die die Anzahlen der Behandlung und Kontrolle räumlich wie zeitlich zusamenfasst. Nämlich einmal zu Beginn des Versuches und einmal zum Ende. Henderson-Tilton erlaubt es damit zusätzlich Unterschiede in der Populationsdynamik zu bewerten.

$$
WG_{HT} = \left(1 - \cfrac{T_a}{T_b}\cdot\cfrac{C_b}{C_a}\right) 
$$

mit

-   $T_{a}$ Anzahl lebend in behandelter Parzelle *nach* Applikation (eng. *after*)
-   $T_{b}$ Anzahl lebend in behandelter Parzelle *vor* Applikation (eng. *before*)
-   $C_{a}$ Anzahl lebend in Kontrolle *nach* Applikation (eng. *after*)
-   $C_{b}$ Anzahl lebend in Kontrolle *vor* Applikation (eng. *before*)

Hier nochmal ein kleines numerisches Beispiel mit $T_b = 300$ sowie $T_a = 30$ und $C_b = 500$ sowie $C_a = 1000$. Wir fangen also in einer Population mit einem Befall von $500$ Trepsen in der Kontrolle an und haben am Ende des Versuches dann $1000$ Trepsen in unserem Versuch vorliegen. In der Behandlung haben wir zu Beginn des Versuches $300$ Trepsen im Feld gezählt und finden nach der Applikation unseres Pflanzenschutzmittels am Ende des Versuches noch $30$ Trepsen wieder. Dann können wir einmal die Zahlen in die Formel einsetzen.

$$
WG_{HT} = \left(1 - \cfrac{T_a}{T_b}\cdot\cfrac{C_b}{C_a}\right)  = \left(1 - \cfrac{30}{300}\cdot\cfrac{500}{1000}\right) = 1 - (0.1 \cdot 0.5) = 0.95 
$$

Wir haben hier nun einen räumlichen und zeitlichen Bezug in unserem Wirkungsgrad abgebildet. Leider haben wir nicht immer einen Befall in unserem Bestand. Der Befall tritt erst während der Applikation auf. Dann müssen wir die Formel nach Abbott nutzen, die sich nur auf die Entwicklung des Befalls am Ende des Versuches konzentriert.

### Wirkungsgrad nach @abbott1925method {#sec-effect-abbott}

Im Folgenden wollen wir uns auf die Berechnung des Wirkungsgrad nach Abbott konzentrieren. Der Wirkungsgrad wird sehr häufig genutzt, so dass sich hier eine tiefere Betrachtung lohnt. Wir berechnen hier den Wirkungsgrad nach @abbott1925method erstmal sehr simpel an einem Beispiel. Dann nutzen wir R um die ganze Berechnung nochmal etwas komplexer durchzuführen. Abschließend schauen wir uns mit der Anpassung von @finner1989berechnung noch eine Erweiterung des Wirkungsgrads nach Abbott an. Der Wirkungsgrad $WG$ eines Schutzmittels im Vergleich zur Kontrolle berechnet sich wie folgt.

$$
WG_{Ab} = \left(\cfrac{C_a - T_a}{C_a}\right) = \left(1-\cfrac{T_a}{C_a}\right)
$$

mit

-   $T_a$ Anzahl lebend in der Behandlung *nach* Applikation (eng. *after*)
-   $C_a$ Anzahl lebend in der Kontrolle *nach* Applikation (eng. *after*)

Wir können auch nur den räumlichen Bezug mit dem Wirkungsgrad nach Abbott abbilden. Hier nochmal ein kleines numerisches Beispiel von oben mit $T_b = 300$ sowie $T_a = 30$ und $C_b = 500$ sowie $C_a = 1000$. Wir fangen also in einer Population mit einem Befall von $500$ Trepsen in der Kontrolle an und haben am Ende des Versuches dann $1000$ Trepsen in unserem Versuch vorliegen. In der Behandlung haben wir zu Beginn des Versuches $300$ Trepsen im Feld gezählt und finden nach der Applikation unseres Pflanzenschutzmittels am Ende des Versuches noch $30$ Trepsen wieder. Hier ist es wichtig, dass wir uns nicht den Zusammenhang über die Zeit anschauen können. Also entweder betrachten wir den Zusammenhang von *nach* unser Applikation oder aber den Vergleich von vorher und nachher von unserer Behandlung.

$$
WG_{Ab} = \left(1-\cfrac{T_a}{C_a}\right) = \left(1-\cfrac{30}{1000}\right) = 0.97
$$

Oder wir konzentrieren uns eben nur auf den zeitlichen Bezug in dem wir die Formel nach Abbott etwas modidifzieren und nur auf die Behandlung $T$ schauen.

$$
WG_{Ab^*} = \left(1-\cfrac{T_a}{T_b}\right) = \left(1-\cfrac{30}{300}\right) = 0.90
$$

Natürlich ist die Formel wieder auf sehr einfachen Daten gut zu verstehen. Wir haben aber später dann in der Anwendung natürlich Datensätze, so dass wir uns hier einmal zwei Beispieldaten anschauen wollen. Zuerst schauen wir uns einen Datensatz zu dem Befall mit Trespe, einem Wildkraut, an. Wir haben also Parzellen in denen sich die Trespe ausbreitet und haben verschiedene Behandlungen durchgeführt. Wichtig hierbei, wir haben auch Parzellen wo wir nichts gemacht haben, das ist dann unsere positive Kontrolle (*ctrl*).

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-effect-trespe
#| tbl-cap: "Trespebefall in einem randomisierten Blockdesign mit vier Behandlungsvarianten."

trespe_tbl <- read_excel("data/raubmilben_data.xlsx", sheet = "trespe") 

trespe_tbl %>% 
  kable(align = "c", "pipe")
```

In der @tbl-effect-trespe-abbott haben wir dann einmal den Wirkungsgrad nach Abbott `wg_abbott` aus dem mittleren Befall über die vier Blöcke berechnet. Wir lassen hier mal weg, dass wir in einigen Parzellen mehr Befall hatten und deshab negative Wirkunsggrade. Hier hilft dann der Mittelwert.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-effect-trespe-abbott
#| tbl-cap: "Trespebefall in einem randomisierten Blockdesign mit vier Behandlungsvarianten ergänzt um den Wirkungsgrad nach Abbott `wg_abbott` berechnet aus dem mittleren Befall über die Blöcke."

trespe_tbl |> 
  rowwise() |> 
  mutate(mean = mean(c_across(block_1:block_4))) |> 
  ungroup() |> 
  mutate(wg_abbott = round(1 - mean/mean[1], 2)) |> 
  kable(align = "c", "pipe")
```

Dann das Ganze auch nochmal händisch aufbereitet und berechnet. Wir erhalten die gleichen Werte wie auch oben in der Tabelle.

$$
\begin{align}
WG_2 &= \left(1 - \cfrac{2.8}{12}\right) = 0.77 \\
WG_3 &= \left(1 - \cfrac{3.2}{12}\right) = 0.73 \\
WG_4 &= \left(1 - \cfrac{7.0}{12}\right) = 0.42
\end{align}
$$

Jetzt wollen wir etwas tiefer in die Thematik einsteigen und müssen daher einmal die Daten umwandeln, da unsere Daten nicht im Long-Format vorliegen. Wir müssen die Daten erst noch anpassen und dann die Spalte `block` in einen Faktor umwandeln. Laden wir einmal den Datensatz in R und verwandeln das Wide-Format in das Long-Format. Dann natürlich wie immer alle Faktoren als Faktoren mutieren.

```{r}
trespe_tbl <- read_excel("data/raubmilben_data.xlsx", sheet = "trespe") %>% 
  pivot_longer(block_1:block_4,
               names_to = "block",
               values_to = "count") %>% 
  mutate(block = factor(block, labels = 1:4),
         variante = as_factor(variante))
```

In der @fig-effect-1 sehen wir nochmal die orginalen, untransformierten Daten sowie die $log$-transformierten Daten. Wir nutzen die Funktion `log1p()` um die Anzahl künstlich um 1 zu erhöhen, damit wir Nullen in den Zähldaten zum Logarithmieren vermeiden.

```{r}
#| echo: true
#| message: false
#| warning: false
#| fig-align: center
#| fig-height: 3
#| fig-width: 5
#| label: fig-effect-1
#| fig-cap: "Dotplot der Anzahl an Trespen je Sorte und Block."
#| fig-subcap: 
#|   - "Verteilung der Werte auf der originalen Skala."
#|   - "Verteilung der Werte auf der logarithmischen Skala. Beobachtungen mit einer 0 Zählung wurden auf 1 gesetzt."
#| layout-nrow: 1
#| column: page

ggplot(trespe_tbl, aes(variante, count, fill = block)) +
  theme_minimal() +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  labs(x = "Behandlungsvariante", y = "Anzahl", fill = "Block") +
  scale_fill_okabeito()
  
ggplot(trespe_tbl, aes(variante, log1p(count), fill = block)) +
  theme_minimal() +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  labs(x = "Behandlungsvariante", y = "log(Anzahl)", fill = "Block") +
  scale_fill_okabeito()
  
```

Nun geht es eigentlich ganz fix den Wirkungsgrad nach Abbott mit einer linearen Poisson Regression zu berechnen. Erstmal schätzen wir eine Poissonregression mit der Anzahl der Trespen als Outcome. Es ist wichtig in der Funktion `glm()` die Option `family = poisson` zu setzen. Sonst rechnen wir auch keine Poissonregression. Wir haben hier im Prinzip die effizientere Variante vorliegen, da wir durch die Modellierung auch noch andere Einflussfaktoren mit ins Modell nehemn können als nur den Block. Da wir hier aber nichts Weiteres gemessen haben, bleibt es bei dem einfachen Modell.

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

fit <- glm(count ~ variante + block, data = trespe_tbl, family = poisson)
```

Wir nutzen die Schätzer aus dem Modell um mit der Funktion `emmeans()` die Raten in jeder Variante gemittelt über alle Blöcke zu berechnen. Dann müssen wir nur noch die Formel nach Abbott nutzen um jede Rate in das Verhältnis zur Rate der Kontrolle `rate[1]` zu setzen. Wir erhalten dann den Wirkungsgrad nach Abbott für unsere drei Varianten. Wie wir sehen, sind auch hier die Werte gleich.

```{r}
#| message: false
#| warning: false
res_trespe_tbl <- fit %>% 
  emmeans(~variante, type = "response") %>% 
  tidy() %>% 
  mutate(WG_abbott = (rate[1] - rate)/rate[1],
         WG_abbott_per = percent(WG_abbott)) %>% 
  select(variante, rate, WG_abbott, WG_abbott_per)
res_trespe_tbl
```

Damit haben wir den Wirkungsgard $WG_{abbott}$ für unser Trepsenbeispiel einmal berechnet. Die Interpretation ist dann eigentlich sehr intuitiv. Wir haben zum Beispiel bei Variante 2 einen Wirkungsgrad von 76.7% der Kontrolle und somit auch eine Reduzierung von 76.7% der Trepsen auf unseren Parzellen im Vergleich zur Kontrolle.

Wir können auch einfach testen, ob sich die Kontrolle von den anderen Varianten sich unterscheidet. Wir nutzen dafür wieder die Funktion `emmeans()` und vergleichen alle Varianten zu der Kontrolle. Da die Kontrolle in der ersten Zeile steht bzw. das erste Level des Faktors `variante` ist, müssen wir noch die Referenz mit `ref = 1` setzen.

```{r}
fit %>% 
  emmeans(trt.vs.ctrlk ~ variante, type = "response", ref = 1) %>% 
  pluck("contrasts")
```

Wir wir sehen, sind alle Varianten, bis auf die Variante 4, signifikant unterschiedlich zu der Kontrolle. Nun könnten wir uns auch Fragen in wie weit sich die Wirkungsgrade untereinander unterscheiden. Dafür müssen wir dann die Anteile mit der Funktion `pairwise.prop.test()` paarweise Vergleichen. Da wir natürlich auf die Kontrolle standardisieren hat der Vergleich mit der Kontrolle nur so halbwegs Sinn, mag aber von Interesse sein.

```{r}
pairwise.prop.test(x = res_trespe_tbl$WG_abbott * 100, 
                   n = c(100, 100, 100, 100),
                   p.adjust.method = "none") 
```

Wir sehen, dass sich alle Varianten in ihrem Wirkungsgrad von der Kontrolle unterscheiden, die Kontrolle hat hier die ID 1. Die Variante 4 unterscheidet sich von allen anderen Varianten in ihrem Wirkungsgrad.

Im zweiten Beispiel wollen wir uns mit dem geometrischen Mittel $WG_{geometric}$ als Schätzer für den Wirkungsgrad beschäftigen. Hier kochen wir dann einmal die Veröffentlichung von @finner1989berechnung nach. Dafür brauchen wir einmal die Daten zu den Raubmilben, die ich schon als Exceldatei aufbereitet habe. Wie immer sind die Rohdaten im Wide-Format, wir müssen aber im Long-Format rechnen. Da bauen wir uns also einmal schnell die Daten um. Dann wollen wir noch die Anzahlen der Raubmilben logarithmieren, so dass wir jede Anzahl um 1 erhöhen um logarithmierte Nullen zu vermeiden. Das ganze machen wir dann in einem Rutsch mit der Funktion `log1p()`.

```{r}
mite_tbl <- read_excel("data/raubmilben_data.xlsx", sheet = "mite") %>% 
  pivot_longer(block_1:block_5,
               names_to = "block",
               values_to = "count") %>% 
  mutate(block = factor(block, labels = 1:5),
         sorte = as_factor(sorte),
         log_count = log1p(count))
```

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-effect-mite
#| tbl-cap: "Raubmilbenbefall auf acht Sorten und einer Kontrolle."

read_excel("data/raubmilben_data.xlsx", sheet = "mite") %>% 
  kable(align = "c", "pipe")
```

Schauen wir uns einmal die Daten in der @fig-effect-2 an. Wichtig ist, dass wir usn merken, dass in unseren Daten jetzt die Kontrolle an der neuten Position ist. Das ist eben in diesen Daen so, ich hätte die Sortierung dann anders gemacht. Aber wir arbeiten hier mit dem was wir vorliegen haben.

```{r}
#| echo: true
#| message: false
#| warning: false
#| fig-align: center
#| fig-height: 3
#| fig-width: 5
#| label: fig-effect-2
#| fig-cap: "Dotplot der Anzahl an Raubmilden je Sorte und Block."
#| fig-subcap: 
#|   - "Verteilung der Werte auf der originalen Skala."
#|   - "Verteilung der Werte auf der logarithmischen Skala. Beobachtungen mit einer 0 Zählung wurden auf 1 gesetzt."
#| layout-nrow: 1
#| column: page

ggplot(mite_tbl, aes(sorte, count, fill = block)) +
  theme_minimal() +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  labs(x = "Sorten", y = "Anzahl", fill = "Block") +
  scale_fill_okabeito()
  
ggplot(mite_tbl, aes(sorte, log1p(count), fill = block)) +
  theme_minimal() +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  labs(x = "Sorten", y = "log(Anzahl)", fill = "Block") +
  scale_fill_okabeito()
  
```

Wir brauchen jetzt eine Helferfunktion, die uns aus $Pr$ die Gegenwahrscheinlichkeit $1 - Pr$ berechnet. Auch wollen wir dann die Prozentangabe der Gegenwahrscheinlichkeit, also die Gegenwahrscheinlichkeit $1 - Pr$ multipliziert mit Einhundert. Dann brauchen wir als Variable noch die Gruppengröße $n_g$, die bei uns ja bei 5 liegt. Wir haben pro Sorte fünf Beobachtungen je Block.

```{r}
get_q <- function(x){100 * (1 - x)}
n_group <- 5
```

Wir nutezn jetzt das [geometrisches Mittel](#sec-desc-stat-geometric) um den Effekt der Behandlung bzw. Sorte im Verhältnis zur Kontrolle zu berechnen. Hierbei ist es wichtig sich zu erinnern, dass wir nicht alle paarweisen Vergleiche rechnen, sondern nur jede Sorte $j$ zu der Kontrolle $ctrl$ vergleichen. Dabei nutzen wir dann das Verhältnis der geometrisches Mittel um zu Beschreiben um wie viel weniger Befall mit Raubmilden wir in den Sorten $y_j$ im Verhältnis zur Kontrolle $y_{ctrl}$ vorliegen haben.

$$
\Delta_{geometric} = \left(\cfrac{\prod_{i=1}^n y_j}{\prod_{i=1}^n y_{ctrl}}\right)^{1/n_j} \mbox{ für Sorte } j 
$$

Berechnen wir also als erstes einmal das Produkt aller gezählten Raubmilden pro Sorte und speichern das Ergebnis in der Spalte `prod`.

```{r}
mite_wg_gemetric_tbl <- mite_tbl %>% 
  mutate(count = count + 1) %>% 
  group_by(sorte) %>% 
  summarise(prod = prod(count)) 
mite_wg_gemetric_tbl
```

Jetzt können wir im nächsten Schritt einmal das $\Delta_{geometric}$ berechnen und dann den Wirkungsgrad über die Gegenwahrscheinlichkeit. Wichtig ist hier, dass die Kontrolle in der neunten Zeile ist. Daher teilen wir immer durch das Produkt an neunter Position mit `prod[9]`. Wir sehen ganz klar, das wir ein Delta von 1 für die Kontrolle erhalten, da wir ja die Kontrolle ins Verhältnis zur Kontrolle setzen. Damit sollte der Rest auch geklappt haben. Den Wirkungsgrad der Sorten in Prozent gegen Raubmilbenbefall im Verhältnis zur Kontrolle können wir dann direkt ablesen. Die Funktion `percent()` berechnet uns aus den Gegenwahrscheinlichkeiten dann gleich die Prozent. Wir brauchen daher hier noch nicht unsere Funktion `get_q()`.

```{r}
mite_wg_gemetric_tbl %>% 
  mutate(delta_geometric = (prod/prod[9])^(1/n_group),
         WG_geometric = percent(1 - delta_geometric))
```

Soweit haben wir erstmal nur eine andere Variante des Wirkungsgrades berechnet. Im Gegensatz zu der Berechnung nach @abbott1925method können wir aber bei den geometrischen Wirkungsgrad auch ein einseitiges 95% Konfidenzintervall $[-\infty; upper]$ angeben. Dafür müssen wir erstmal den Exponenten $a$ berechnen und mit diesem dann die obere Konfidenzschranke $upper$. Dafür brauchen wir dann doch ein paar statistische Maßzahlen.

$$
\begin{align}
a &= \sqrt{2/n} \cdot s \cdot t_{\alpha=5\%} + \ln(1 - \Delta_{geometric})\\
upper &= 1 - e^a
\end{align}
$$

Wir brauchen zum einen die Freiheitsgrade der Residuen `df.residual` sowie den Standardfehler der Residuen `sigma` oder $s$. Beides erhalten wir aus einem linearen Modell auf den logarithmierten Anzahlen der Raubmilben.

```{r}
residual_tbl <- lm(log_count ~ sorte + block, data = mite_tbl) %>% 
  glance() %>% 
  select(df.residual, sigma)
residual_tbl
```

Mit den Freiheitsgraden der Residuen können wir jetzt den kritischen Wert $t_{\alpha = 5\%}$ aus der $t$-Verteilung berechnen.

```{r}
t_quantile <- qt(p = 0.05, df = residual_tbl$df.residual, lower.tail = FALSE)
t_quantile
```

Wir können jetzt unseren Datensatz `mite_wg_gemetric_tbl` um die Spalte mit den Werten des Exponenten $a$ ergänzen und aus diesem dann die obere Schranke des einseitigen 95% Konfidenzintervall berechnen. Die Werte für die Kontrolle ergeben keinen biologischen Sinn und sind ein mathematisches Artefakt. Wir kriegen halt immer irgendwelche Zahlen raus.

```{r}
#| message: false
#| warning: false
mite_res_tbl <- mite_wg_gemetric_tbl %>% 
  mutate(delta_geometric = (prod/prod[9])^(1/n_group),
         WG_geometric = get_q(delta_geometric),
         a = sqrt(2/n_group) * residual_tbl$sigma * t_quantile + log(delta_geometric),
         upper = get_q(exp(a)))
mite_res_tbl
```

Wir räumen nochmal die Ausgabe auf und konzentrieren uns auf die Spalten und runden einmal die Ergebnisse.

```{r}
#| message: false
#| warning: false
mite_res_tbl %>% 
  select(sorte, WG_geometric, upper) %>% 
  mutate(across(where(is.numeric), ~round(.x, 2))) 
```

Wir können wir jetzt das einseitige 95% Konfidenzintervall interpretieren? In der Sorte 1 erhalten wir einen Wirkungsgrad von 93.6% und mit 95% Sicherheit mindestens einen Wirkungsgrad von 82.1%. Damit haben wir auch eine untere Schranke für unseren Wirkungsgradschätzer. Wir berichten für die Sorte 1 einen Wirkungsgrad von 93.6% \[$-\infty$; 82.1%\].

Genauso können wir aber auch den Wirkungsgrad nach @abbott1925method nochmal auf den Daten berechnen. Dafür müssen wir nur eine Poissonregression auf der Anzahl der Raubmilben rechnen. Die angepassten Werte können wir dann verwenden um den $WG_{abbott}$ zu schätzen.

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

fit <- glm(count ~ sorte + block, data = mite_tbl, family = poisson)
```

Wir nutzen hier dann die Funktion `emmeans()` um die mittlere Anzahl an Raubmilben über alle Blöcke zu ermitteln. Dann können wir auch schon den $WG_{abbott}$ berechnen. Wichtig ist hier, dass die Referenzkategorie mit der Kontrolle an der neunten Stelle bzw. Zeile steht. Deshalb müssen wir auch hier durch `prod[9]` teilen. Achtung, die $p$-Werte haben hier keine tiefere Bedeutung im Bezug auf den Wirkungsgrad und deshalb schauen wir uns diese Werte auch gar nicht tiefer an.

```{r}
#| message: false
#| warning: false
res_mite_tbl <- fit %>% 
  emmeans(~sorte, type = "response") %>% 
  tidy() %>% 
  mutate(WG_abbott = (rate[9] - rate)/rate[9],
         WG_abbott_per = percent(WG_abbott)) %>% 
  select(sorte, WG_abbott, WG_abbott_per, rate)
res_mite_tbl
```

Wenn du verwirrt bist über die negativen Wirkungsgrade, dann musst du dir klar werden, dass wir immer den Wirkungsgrad im Verhältnis zur Kontrolle berechnen. Wenn du also negative Wirkungsgrade siehst, dann sind die Anzahlen in den Sorten *höher* als in der Kontrolle. Je nach Fragestellung macht dieses Ergebnis mehr oder weniger Sinn.

Auch hier können wir noch schnell einen statistischen Test rechnen und die Sorten zu der Kontrolle `ref = 9` vergleichen.

```{r}
fit %>% 
  emmeans(trt.vs.ctrlk ~ sorte, type = "response", ref = 9) %>% 
  pluck("contrasts")
```

Wir wir sehen, sind alle Varianten signifikant unterschiedlich zu der Kontrolle bis auf die Sorte 3.

Auch hier können wir die Wirkungsgrade miteinander vergleichen. Aber Achtung, wir können nur Wirkungsgrade zwischen 0 und 1 miteinander sinnvoll vergleichen. Deshalb fliegen bei uns einige an Sorten raus. Darüber hinaus musst du jetzt auch einmal schauen, welche Sorte sich hinter der numerischen ID aus dem `pairwise.prop.test()` verbirgt.

```{r}
res_mite_tbl %>% 
  filter(WG_abbott < 1 & WG_abbott >= 0) %$% 
  pairwise.prop.test(x = WG_abbott * 100, 
                   n = rep(100, 5),
                   p.adjust.method = "none") 
```

Das letzte Beispiel zeigt nochmal schön, dass wir uns immer überlegen müssen, was wir mit den Wirkungsgraden eigentlich zeigen wollen und welche Wirkungsgrade noch eine biologisch sinnvolle Bedeutung haben.

### Wirkungsgrad nach Sun-Shepard {#sec-effect-ss}

Der Wirkungsgrad nach Sun-Shepard bezieht sich auf die Mortalität \[%\] in unserem Versuch. Wir können müssen hier auch vor der Applikation und nach der Applikation in unserer Kontrolle die Mortalität bestimmen, sonst können wir nicht das $C_{\Delta}$ in der Kontrolle berechnen. Auf der [Webseite von Ehab Bakr](https://www.ehabsoft.com/ldpline/onlinecontrol.htm), auf die sich irgendwie in einem Zirkelschluss alle meine gefundenen Quellen beziehen, steht dann zwar $\pm$ aber ich denke, dass hier ausgesagt werden soll, dass $C_{\Delta}$ negativ oder positiv sein kann. Dann können wir auch auch besser nur Plus schreiben und damit ist dann die Formel auch in sich konsistent. Das passt dann auch zu der Durchführung und dem Beispiel auf der Webseite.

$$
WG_{SS} = \left(\cfrac{T + C_{\Delta}}{1 + C_{\Delta}}\right) 
$$

und

$$
C_{\Delta} = \left(\cfrac{C_a - C_b}{C_b}\right) 
$$

mit

-   $T$ Sterberate in der Behandlung
-   $C_b$ Anzahl lebend in der Kontrolle *nach* Applikation (eng. *after*)
-   $C_a$ Anzahl lebend in der Kontrolle *vor* Applikation (eng. *before*)

### Wirkungsgrad nach Schneider-Orelli {#sec-effect-so}

Der Wirkungsgrad nach Schneider-Orelli betrachtet nur die Mortalität nach dem Versuch und ist somit mit dem Wirkungsgrad nach Abbott zu vergleichen. Die Formel ist dabei sehr simpel

$$
WG_{SO} = \left(\cfrac{T_a - C_a}{1 - C_a}\right) 
$$

mit

-   $T_a$ Sterberate in der Behandlung *nach* Applikation (eng. *after*)
-   $C_a$ Sterberate in der Kontrolle *nach* Applikation (eng. *after*)

## Interrater Reliabilität {#sec-effect-interrater}

Worum geht es jetzt in diesem Abschnitt? Wir wollen uns anschauen, ob verschiedene Bewerter auf die gleichen Noten oder Einschätzungen kommen. Das ist eher selten im Kontext einer Abschlussarbeit, aber kommt häufiger bei der Bonitur vor. Wir können uns dabei zwei Fälle vorstellen. Zum einen wollen wir wissen ob zwei oder mehr Bewerter auf die gleichen Noten kommen, wenn sie das gleiche Subjekt bewerten. Wir haben also eine Rasenfläche oder ein Schwein vor uns und drei Forschende sollen bestimmen wie krank der Rasen oder das Schwein auf einer Skala von 1 bis 9 ist. Vorher wurden natürlich alle drei auf der Notensakal trainiert. Das würden wir dann die Interrater Reliabilität nennen. Benoten zwei oder mehr Bewerter ein Subjekt nach einem gegebenen Notenschlüssen und Kriterien gleich. Wichtig ist hier, dass es ein Kriterienkatalog gibt. Wir machen machen hier also nicht auf Geratewohl los.

Die andere Möglichkeit ist, dass wir wissen wollen, ob ein Bewerter selber konstante Benotungen an dem gleichen Subjekt liefert. Das können wir durch Fotos und zufällige Zuordnung machen oder eben eine Stunde oder Zeit dazwischen lassen. Wir sprechen dann von Intrarater Reliabilität. Diese Art der Reliabilität ist etwas seltener und kann auch gleich gelöst werden wie der Vergleich von unterschiedlichen Bewertern. Meistens reicht es dann auch aus und wir müssen die Korrelation durch den gleichen Bewerter eigentlich nicht weiter beachten.

Interrater Reliabilität

:   Unterscheiden sich die Bewertungen oder Benotungen an einer Beobachtung zwischen zwei oder mehr Bewertern?

Intrarater Reliabilität

:   Unterscheiden sich die Bewertungen oder Benotungen eines Bewerters an einer Beobachtung über die Zeit?

Reliabilität ist der Grad an Übereinstimmung einer Messung, wenn diese unter identischen Bedingungen jedoch von unterschiedlichen Bewertern wiederholt wird. Anhand der Interrater Reliabilität wird es somit möglich, den Grad zu ermitteln, in dem sich die erhaltenen Ergebnisse wiederholen lassen. Eine wunderbare Übersicht zur Interrater Reliabilität gibt auf Data Novia mit [Introduction to R for Inter-Rater Reliability Analyses](https://www.datanovia.com/en/lessons/introduction-to-r-for-inter-rater-reliability-analyses/). Ebenso liefert das Tutorium [Inter-rater-reliability and Intra-class-correlation](http://dwoll.de/rexrepos/posts/interRaterICC.html) einen guten Überblick über alle Methoden in R. Wir schauen uns jetzt einmal die häufigsten Maßzahlen für die Reliabilität an. Für das erste beschränken wir uns hier einmal auf Cohen's Kappan.

### ... für kategorielle Variablen

Fangen wir mit dem simpelsten Fall an. Wir haben zwei Bewerter $R_1$ und $R_2$ und wir wollen schauen, ob diese beiden Bewerter ein Subjekt richtig in $ja$ oder $nein$ einordnen. So kann die Frage sein, hat die Sonnenblume Mehltau? Oder aber die Frage, ist das Schwein erkrankt? Wir haben hier keine Wahrheit vorliegen, sondern wollen einfach nur wissen, ob unsere beiden Bewerter zum gleichen Schluss kommen. Dafür können wir dann Cohen's Kappa $\kappa$ nutzen. Dabei ist die Verwendung von Cohen's Kappa nicht unumstritten. @maclure1987misinterpretation diskutiert in seiner Arbeit [Misinterpretation and misuse of the kappa statistic](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=63e3fc4cf928bc0a12f6c62cace8d94aee34db62) diesen Sachverhalt. Fangen wir einmal mit einem Beispiel an. Die @tbl-2x2-kappa zeigt einmal das Ergebnis von siebzig bewerteten Subjekten. Wir haben hier zwei Bewerter vorliegen, die jeweils mit $ja$ oder $nein$ bewerten sollten. In R kannst du dir später so eine Tabelle einfach mit der Funktion `xtabs()` bauen.

|                |            |                   |                   |
|:--------------:|:----------:|:-----------------:|:-----------------:|
|                |            |  **Bewerter 1**   |                   |
|                |            |     *ja (1)*      |    *nein (0)*     |
| **Bewerter 2** |  *ja (1)*  | $25_{\;\Large a}$ | $10_{\;\Large b}$ |
|                | *nein (0)* | $15_{\;\Large c}$ | $20_{\;\Large d}$ |

: Eine 2x2 Tabelle als Beispiel für die Bewertung einer Erkrankung an siebzig Beobachtungen von zwei Bewertern. Für die folgenden Berechnungen sind die Zellen mit Buchstaben versehen. {#tbl-2x2-kappa}

Als Faustregel für Cohen's Kappa $\kappa$ dient dabei, dass die Diagonale von links oben nach rechts unten die Übereinstimmungen zeigt. Die Diagonale von links unten nach rechts oben zeigt dann die Nichtübereinstimmung. Das ist natürlich nur eine grobe Abschätzung im Folgenden einmal die Formel mit der du dann Cohen's Kappa $\kappa$ berechnen kannst.

$$
\kappa = \cfrac{p_0 - p_e}{1 - p_e}
$$

mit

-   $p_0$ Anteil der beobachteten Übereinstimmung
-   $p_e$ Anteil der veränderten Übereinstimmung

Das ist jetzt natürlich wieder sehr theoretisch. Deshalb wollen wir einmal die beiden Teile $p_0$ und $p_e$ aus den obigen Daten berechnen. Wenn du jetzt denkst, dass dir das doch irgendwie bekannt vorkommt, ja der $\mathcal{X}^2$-Test sieht sehr ähnlich aus. Wir können $p_0$ wie folgt aus der 2x2 Kreuztabelle berechnen. Wir schauen hier auf die Übereinstimmungen auf der Diagonalen von obeb links nach unten rechts.

$$
p_0 = \cfrac{a + d}{N} = \cfrac{25 + 20}{70} = 0.643 
$$

Für den Wert von $p_e$ brauchen wir zwei Werte, die wir dann aufaddieren. Einmal berechnen wir den Wert für $p_e$ von $ja$ und einmal den Wert von $p_e$ für $nein$. Das geht dann mit der Formel auch recht schnell. Wir schauen hier im Prinzip auf die Diagonale von unten links nach oben rechts.

$$
p_e^{ja} = \cfrac{a + b}{N} \cdot \cfrac{a + c}{N} = \cfrac{25+10}{70} \cdot \cfrac{25+15}{70} = 0.5 \cdot 0.57 = 0.285 
$$

$$
p_e^{nein} = \cfrac{c + d}{N} \cdot \cfrac{b + d}{N} = \cfrac{15+20}{70} \cdot \cfrac{10+20}{70} = 0.5 \cdot 0.428 = 0.214 
$$

Dann addieen wir die beiden Werte für $p_e$ einmal auf.

$$
p_e = p_e^{ja} + p_e^{nein} = 0.285 + 0.214 = 0.499
$$

Dann können wir auch schon Cohen's Kappa $\kappa$ für eine 2x2 Kreuztabelle berechnen. Ehrlicherweise, machen wir das dann natürlich gleich in R.

$$
\kappa = \cfrac{0.643 - 0.499}{1 - 0.499} = 0.28
$$

Dabei interpretieren wir dann Cohen's Kappa $\kappa$ wie folgt. Beachte dabei immer, dass es sich natürlich erstmal nur um eine grobe Abschätzung handelt. Wie gut oder wie schlecht dann diese Zahl von Cohen's Kappa $\kappa$ in deinem Bereich ist, liegt dann auch daran wie genau du arbeiten möchtest.

Cohen's Kappa $\kappa$ interpretieren

:   Cohen's Kappa kann von -1 (keine Übereinstimmung) bis +1 (perfekte Übereinstimmung) reichen. Wenn $\kappa = 0$ ist, ist die Übereinstimmung nicht besser als das, was durch Zufall erreicht werden würde. Wenn $\kappa$ negativ ist, ist die Übereinstimmung geringer als die zufällig erwartete Übereinstimmung. Wenn $\kappa$ positiv ist, übersteigt die Übereinstimmung der Bewerter die zufällige Übereinstimmung.

Was sind jetzt die exakten grenzen für die Cohen's Kappa Werte? Das kannst du dir fast wieder aussuchen wie du möchtest. Es gibt eine Reihe von Publikationen, die alle verschiedene Grenzwerte liefern. Gerne kannst du dir die Tabellen in dem Tutorium [Cohen's Kappa in R: For Two Categorical Variables](https://www.datanovia.com/en/lessons/cohens-kappa-in-r-for-two-categorical-variables/) einmal näher anschauen. Ich würde dann mal als Quintessenz folgende Regel von @mchugh2012interrater nehmen.

| Cohen's Kappa | Level of agreement | \% of data that are reliable |
|:-------------:|:------------------:|:----------------------------:|
|   0 - 0.20    |        None        |            0 - 4‰            |
|  0.21 - 0.39  |      Minimal       |           4 - 15%            |
|  0.40 - 0.59  |        Weak        |           15 - 35%           |
|  0.60 - 0.79  |      Moderate      |           35 - 63%           |
|  0.80 - 0.90  |       Strong       |           64 - 81%           |
|  Above 0.90   |   Almost Perfect   |          82 - 100%           |

: Interpretation des Cohen's Kappa Wert nach @mchugh2012interrater {#tbl-interpret-kappa}

Das [R Paket `{vcd}`](https://www.datavis.ca/courses/VCD/vcd-tutorial.pdf) liefert nun die Funktion `Kappa()` mit der wir auf einer 2x2 Kreuztabelle zum Beispiel aus der Funktion `xtabs()` dann Cohen's Kappa $\kappa$ berechnen können. Die Funktion `xtabs()` hilft dir aus zwei oder mehr Spalten eine Kreuztabelle zu bauen.

```{r}
xtab <- as.table(rbind(c(25, 10), c(15, 20)))
```

Dann können wir auch schon Cohen's Kappa $\kappa$ in R berechnen.

```{r}
kappa_res <- Kappa(xtab)
kappa_res
```

Die ungewichtete Version entspricht dem Cohen's Kappa, mit dem wir uns in diesem Kapitel beschäftigen. Das gewichtete Kappa sollte nur für ordinale Variablen betrachtet werden und wird weitgehend hier ignoriert. Dann noch die Konfidenzintervalle und alles ist soweit da, was wir brauchen.

```{r}
kappa_res |> 
  confint()
```

Das Ganze geht natürlich nicht nur für zwei Level (ja/nein) sondern dann in R auch für mehrere Level. Also wenn wir zum Beispiel dann Boniturnoten einmal bewerten und dann für zwei Bewerter vergleichen wollen.

Cohen's Kappa für mehr als zwei Bewerter liefert das [R Paket `{irr}`](https://cran.r-project.org/web/packages/irr/irr.pdf) mit der Funktion `kappam.light()` oder `kappam.fleiss()`, wenn es dann auch noch mehr Kategorien als zwei in den Spalten gibt. Das R Paket ist schon ziemlich alt und hat auch sehr schlechte bis gar keine Hilfeseiten. Wenn man weiß, was man braucht und die Funktion bedeuten, dann kannst du das Paket gut nutzen.

Hier helfen dann auch die passenden Tutorien weiter. Ich habe hier nochmal die beiden anderen Fälle aufgelistet, wenn du eben Noten oder noch andere Kategorien in deinen Daten hast. Im Prinzip änderst du nur deine Funktion die du dann nutzt. Die Bewertung der einzelnen Ausgaben schwankt dann auch leicht. Aber eine Eins bedeutet immer eine perfekte Übereinstimmung.

-   [Weighted Kappa in R: For Two Ordinal Variables](https://www.datanovia.com/en/lessons/weighted-kappa-in-r-for-two-ordinal-variables/) mit der Funktion `Kappa()` aus dem R Paket `{vcd}`. Wir nutzen dann die `Weighted` Zeile in der Ausgabe der Funktion.
-   [Fleiss’ Kappa in R: For Multiple Categorical Variables](https://www.datanovia.com/en/lessons/fleiss-kappa-in-r-for-multiple-categorical-variables/) und der Funktion `kappam.fleiss()` aus dem R Paket `{irr}`.

Damit wären wir hier soweit erstmal fertig. Im nächsten Abschnitt schauen wir uns dann einmal den Fall an, dass wir kontinuierliche Bewertungen haben. Wir haben nämlich gezählt, was ja auch eine Art der Bewertung ist. Nur eben nicht auf einer Notenskala sondern eben kontinuierlich.

### ... für kontinuierliche Variablen

jetzt betrachten wir einen Datensatz indem Entenküken auf verschiedenen Bildern von Studierenden gezählt werden sollten. Insgesamt haben wir vier Studierende als Bewerter, die die Küken dann in zwei voneinander getrennten Sessions bewertet haben. Wir laden also einmal den Datensatz und schauen uns die Zähldaten einmal an.

```{r}
duck_tbl <- read_excel("data/interrater_duck_count.xlsx")
```

In der @tbl-duck-data siehst du einmal einen Auszug aus den Entendaten. Dabei steht die `id` für ein Bild, die Spalte `session` beschreibt die beiden Termine, wo die Studierenden dann die Bilder ausgezählt haben.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-duck-data
#| tbl-cap: "Auszug aus Entenkükendatensatz."

rbind(head(duck_tbl, n = 5),
      rep("...", times = ncol(duck_tbl)),
      tail(duck_tbl, n = 5)) %>% 
  kable(align = "c", "pipe")
```

Da so eine Datentabelle dann imer schwer zu lesen ist, dann hier einmal in der @fig-effect-rel-duck die Visualisierung der Daten. Wir sehen, dass die zweite Session schonmals ehr viel weniger Varianz zeigt und sich die Bewerter sehr viel ähnlicher sind als bei der ersten Session. Darüber hinaus ist der dritte Bewerter in der ersten Session am schlechtesten, da er die größte Varianz und Abweichung zu den anderen Bewertern zeigt.

```{r}
#| echo: true
#| warning: false
#| message: false
#| label: fig-effect-rel-duck
#| fig-align: center
#| fig-height: 6
#| fig-width: 7
#| fig-cap: "Visualisierung der Zählungen an Entenküken über die beiden Sessions der vier Studierenden. Die Zahlen geben die jeweilige Bildnummer an."

duck_tbl |> 
  pivot_longer(cols = rater_1:rater_4,
               values_to = "count",
               names_to = "rater") |> 
  mutate(session = as_factor(session)) |> 
  ggplot(aes(x = rater, y = count, color = session, label = id)) +
  theme_minimal() +
  geom_text(position = position_dodge(0.5), show.legend = FALSE) +
  scale_color_okabeito() +
  labs(x = "Bewerter", y = "Gezählte Anzahl an Entenküken auf dem Bild")

```

Was erwarten wir nun? Zum einen sollte die zweite Session bessere Übereinstimmungen liefern als die erste Session. Immerhin haben ja die Studierenden geübt. Wir nutzen jetzt hier den [Intraclass Correlation Coefficient (ICC) in R](https://www.datanovia.com/en/lessons/intraclass-correlation-coefficient-in-r/) um auf die Reliabilität von Bewertern zurückzuschließen. Das Tutorium liefert nochmal eine Reihe an mehr Informationen. Ich halte hier die Sachlage relativ kurz. Wir wollen hier mal prinzipiell schauen, wie unsere Bewerter so abschneiden. Als zweite Überlegung wollen wir dann noch wissen, ob unsere Bewerter dann auch konsistent sind. Daher die erste Session genauso wie die zweite Session bewertet haben. Aber auch hier Achtung, wir wissen nicht die Wahrheit. Wir schauen nur, ob die beiden Sessions gleich waren. Ein Bewerter, der in der ersten Session versagt und in der zweiten alles richtig zählt, würde hier also schlechter wegkommen, als ein Bewerter, der nichts gelernt hat und in beiden Fällen schlecht zählt. Wir würden dann nur sehen, dass die Bewerter untereinander nicht passen.

Betrachten wir also einmal in den beiden folgenden Tabs die beiden Sessions getrennt voneinander über alle Bewerter hinweg und berechnen den *Intraclass Correlation Coefficient (ICC)*. Der ICC-Wert läuft von 0 bis 1, wobei 0 dann gar keine Übereinstimmung und 1 eine perfekte Übereinstimmung bedeutet. Wir würden sagen, dass ab 0.9 eine perfekte Übereinstimmung und unter 0.5 eine schlechte Übereinstimmung vorliegt. Wichtig ist hier, dass die Funktion `icc()` nur die Spalten mit den Bewertungen will und sonst rein gar nichts.

::: panel-tabset
## Session 1

```{r}
duck_tbl |> 
  filter(session == 1) |> 
  select(matches("rater")) |>
  icc(model = "twoway", type = "agreement", unit = "single")
```

Wir sehen, dass wir eine relativ schlechte Übereinstimmung haben. Der ICC-Wert ist mit 0.304 sehr nahe an der Null. Auch wenn der Test sagt, dasss wir einen signifkanten Unterschied von der Null haben, geht mir ein ICC-Wert von unter 0.8 nicht als Gleich durch. Das deckt sich auch mit unserer obigen Abbildung.

## Session 2

```{r}
duck_tbl |> 
  filter(session == 2) |> 
  select(matches("rater")) |>
  icc(model = "twoway", type = "agreement", unit = "single")
```

Wir sehen, dass wir eine relativ gute Übereinstimmung haben! Der ICC-Wert ist mit 0.806 nahe an der Eins. Das deckt sich auch mit unserer obigen Abbildung, wo wir sehen, dass die Bewerter in der zweiten Session viel mehr auf einer Linie liegen. Zwar ist 0.8 noch nicht perfekt, aber schonmal eins sehr viel besserer Wert.
:::

Jetzt möchte ich mir nochmal die Bewerter einzeln für die jeweiligen Sessions anschauen. Wie unterscheiden sich den die Bewerter, wenn sie ein Bild erneut sehen? Auch hier nochmal, wir haben keinen Goldstandard, wir schauen nur, ob die Bewerter in der zweiten Session genauso zählen wie in der ersten Session. Wir machen das jetzt im Folgenden auf zwei Arten. Einmal kannst du das händsich über `select()` machen. Du wählst einfach immer den Bewerter aus, den du möchtest und führst den Code aus. Du ersetzt also `rater_1` dann durch `rater_2` und dann so weiter. Die zweite Variante ist alles in einem Rutsch mit der Funktion `map()`.

::: panel-tabset
## Für einen Rater -- händisch

```{r}
duck_tbl |> 
  select(id, session, rater_1) |>
  pivot_wider(names_from = session, values_from = rater_1, names_prefix = "session_") |> 
  select(-id) |> 
  icc(model = "twoway", type = "agreement", unit = "single")
```

Wir sehen, dass der erste Studierende extrem konsistent die zweite Session im Vergleich zu der ersten bewertet hat. Der ICC-Wert ist fast bei Eins und somit fast ein perfektes Ergebnis. Auf einer numerischen Skala ist auch kaum noch mehr Spielraum nach oben.

## Für alle Rater -- `map()`

Jetzt spielen wir ein wenig. Zuerst bauen wir uns dann für jeden Bewerter einen Datensatz der nur aus den beiden Sessions besteht. Dafür nutzen wir die Funktion `pivot_wider()`. Dann entfernen wir noch die ID, da die Funktion `icc()` gleich die ID nicht mag. Wir nutzen die Funktion `map()` um für jeden Bewerter den Datensatz zu bauen und am Ende in einer Liste zu speichern.

```{r}
#| message: false
#| warning: false
rater_lst <- map(c("rater_1", "rater_2", "rater_3", "rater_4"), \(x) {
  duck_tbl |> 
  select(id, session, all_of(x)) |>
  pivot_wider(names_from = session, 
              values_from = x, 
              names_prefix = "session_") |> 
  select(-id)
})
```

Da wir alles einer Liste gespeichert haben, können wir jetzt die Liste mit der Funktion `map()` iterativ bearbeiten. Wir rufen dann für jeden der Listeneinträge dann einmal die Funktion `icc()` auf und extrahieren uns durch die Funktion `pluck()` den berechneten ICC Wert. Wenn du die letzte Zeile weglässt, dann kriegst du eben vier mal die Ausgabe der Funktion `icc()`, das war mir dann hier zu viel Platz.

```{r}
rater_lst |> 
  map(~icc(ratings = .x, model = "twoway", type = "agreement", unit = "single")) |> 
  map_dbl(pluck("value"))
```

Wir sehen hier einmal als Werte, was wir auch schon in der @fig-effect-rel-duck beobachten konnten. Je höher die Zahl, desto besser die Übereinstimmung zwischen der ersten und der zweiten Session. Die Reihenfolge ist die Reihenfolge der Bewerter.
:::

Der erste Studierende ist der konsitenteste Bewerter und der vierte Studierende hat die geringste Übereinstimmung. Hier muss man aber vorsichtig sein, da wir ja sehen, dass die erste Session teilweise sehr mies war. Dann kann auch die Übereinstimmung zwischen der ersten und zweiten Session nur schelcht sein, wenn sich der Studierende sehr verbessert hat. Denn eins musst du immer wissen, die Lernleistung oder Verbesserung kannst du hier nicht abbilden. Ein sehr schelchter Bewerter, der aber zügig lernt, wird hier immer mit einer sehr schelchten Übereinstimmung rauskommen. Dagegen müsstest du dann sehr viele Sessions laufen lassen, denn je mehr Sessions er macht, desto weniger unterscheiden sich diese.

## Referenzen {.unnumbered}