stat-modeling-poisson.qmd

```{r}
#| echo: false
#| message: false
#| warning: false
pacman::p_load(tidyverse, readxl, knitr, kableExtra, patchwork)
source("images/R/stat-modeling-preface-00.R")
source("images/R/stat-modeling-preface-02.R")

```

# Poisson Regression {#sec-poisson}

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

In diesem Kapitel wollen wir eine Poisson Regression rechnen. Wir müssen uns hier wieder überlegen, was ist eigentlich unser Outcome $y$ und was sind unsere Einflussvariablen $x$. Die Poisson Regression ist je nach Hintergrund des Anwenders eher selten. In der Ökologie, wo gerne mal gezaählt wird, wie oft etwas vorkommt, ist die Poisson Regression häufig vertreten. Sonst fristet die Poisson Regresson eher ein unbekanntes Dasein.

Ein häufig unterschätzter Vorteil der Poisson Regression ist, dass wir auch auch $0/1$ Daten eine Poisson Regression rechnen können. Moment, wirst du jetzt vielleicht denken, das machen wir doch mit der logistischen Regression. Ja, das stimmt, aber wir können auf Zahlen viel rechnen. Wenn wir auf ein $0/1$ Outcome eine Poisson Regression rechnen, dann kriegen wir nicht Odds Ratios $OR$ als Effektschätzer sondern Risk Ratios $RR$. Wir erhalten also keine Chancen sondern Wahrscheinlichkeiten. Unter der Annahme, dass das Modell auch konvergiert und wir sinnvolle Zahlen erhalten.

Ein weiteres Problem sind die zu vielen Nullen in dem Outcome $y$. Daher wir zählen über die Maßen viel *Nichts*. Wir nennen diesen Fall *zero inflation* und beschreiben damit die zu vielen Nullen in den Daten. Hier muss dann noch speziell modelliert werden. Eine Poisson Regression hat schon so seine speziellen Tücken.

In der folgenden @fig-scatter-modeling-poisson-01 siehst du nochmal den Zusammenhang zwischen einem poissonverteilten Outcome, wie Anzahl oder Auftreten, im Bezug zu einem kontinuierlichen oder kategoriellen $x$ als Einflussgröße. Wir sind dann eher im Bereich der Posthoc-Tests mit kategoriellen $x$ als in wirklich einer Interpretation einer Poisson Regression mit einem kontinuierlichen $x$. Dennoch müssen wir wissen, ob die Poisson Regression gut funktioniert hat. Die Überprüfung der Modellierung können wir dann in diesem Kapitel einmal durchgehen.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-scatter-modeling-poisson-01
#| fig-align: center
#| fig-height: 6
#| fig-width: 15
#| fig-cap: "Visueller Zusammenhang eines kontinuierlichen Outcomes ($y$) aus einer Poissonverteilung zu Zähldaten im Verhätnis zu verschiedenen Skalen der Einflussvariable ($x$). Ein Punkt stellt eine Beobachtung dar. Deutlich sind die Ebenen durch die absoluten Zählwerte zu erkennen. **(A)** $x$ ist kontinuierlich. **(B)** $x$ ist kategoriell mit zwei oder mehr Gruppen. **(C)** $x$ ist kategoriell mit zwei Gruppen. *[Zum Vergrößern anklicken]*"

p21 + p22 + p23 +
  plot_layout(ncol = 3) +
  plot_annotation(tag_levels = 'A')

```

Häufig wollen wir dann die Ergebnisse aus einer Poisson Regression dann im Falle eines kategoriellen $x$ dann als Barplots wie in der @fig-scatter-modeling-poisson-02 darstellen. Hier musst du dann einmal weiter unten in den Abschnitt zu den Gruppenvergleichen springen.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: fig-scatter-modeling-poisson-02
#| fig-align: center
#| fig-height: 5
#| fig-width: 5
#| fig-cap: "Visueller Zusammenhang eines gemittelten Outcomes ($y$) aus einer Poissonverteilung im Verhältnis zu der Einflussvariable ($x$) mit zwei oder mehr Kategorien anhand von Barplots. Hauptsächlich unterscheiden sich die Barplots durch die unterschiedlichen Einheiten auf der $y$-Achse. Die Fehlerbalken stellen den Standardfehler (*SE*) dar. *[Zum Vergrößern anklicken]*"

p02 +
  theme(axis.text.y = element_blank())

```

Im folgenden Kapitel zu der multiplen Poisson linearen Regression gehen wir davon aus, dass die Daten in der vorliegenden Form *ideal* sind. Das heißt wir haben weder fehlende Werte vorliegen, noch haben wir mögliche Ausreißer in den Daten. Auch wollen wir keine Variablen selektieren. Wir nehmen alles was wir haben mit ins Modell. Sollte eine oder mehre Bedingungen nicht zutreffen, dann schaue dir einfach die folgenden Kapitel an.

-   Wenn du fehlende Werte in deinen Daten vorliegen hast, dann schaue bitte nochmal in das @sec-missing zu Imputation von fehlenden Werten.
-   Wenn du denkst, dass du Ausreißer oder auffälige Werte in deinen Daten hast, dann schaue doch bitte nochmal in das @sec-outlier zu Ausreißer in den Daten.
-   Wenn du denkst, dass du zu viele Variablen in deinem Modell hast, dann hilft dir das @sec-variable-selection bei der Variablenselektion.

Daher sieht unser Modell wie folgt aus. Wir haben ein $y$ und $p$-mal $x$. Wobei $p$ für die Anzahl an Variablen auf der rechten Seite des Modells steht. Im Weiteren folgt unser $y$ einer Poissonverteilung. Das ist hier sehr wichtig, denn wir wollen ja eine multiple Poisson lineare Regression rechnen.

$$
y \sim x_1 + x_2 + ... + x_p 
$$

Wir gehen jetzt einmal die lineare Regression mit einem normalverteilten Outcome $y$ einmal durch. Wie schon oben erwähnt, können wir die Poisson Regression als vorangestelltes Modell für eine ANOVA oder dann eben einen Posthoc-Test in `{emmeans}` nutzen.

::: callout-tip
## Weitere Tutorien für die Poisson Regression

Wie immer gibt es auch für die Frage nach dem Tutorium für die Poisson Regression verschiedene Quellen. Ich kann noch folgende Informationen und Hilfen empfehlen.
:::

## Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

```{r echo = TRUE}
#| message: false
pacman::p_load(tidyverse, magrittr, broom,
               parameters, performance, MASS, pscl, see,
               modelsummary, scales, emmeans, multcomp,
               conflicted)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::filter)
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
                "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
```

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

## Daten

Im folgenden schauen wir uns ein Datenbeispiel mit Hechten an. Es handelt sich um langnasige Hechte in nordamerikanischen Flüssen. Wir haben uns insgesamt $n = 68$ Flüsse einmal angesehen und dort die Anzahl an Hechten gezählt. Im Weiteren haben wir dann noch andere Flussparameter erhoben und fragen uns nun, welche dieser Parameter einen Einfluss auf die Anzahl an Hechten in den Flussarmen haben. In @sec-example-longnose findest du nochmal mehr Informationen zu den Daten. Wir entfernen hier die Informationen zu den Flüssen, die brauchen wir in dieser Analyse nicht. Die Daten zu den langnasigen Hechten stammt von [Salvatore S. Mangiafico - An R Companion for the Handbook of Biological Statistics](https://rcompanion.org/rcompanion/e_05.html). Besuche gerne mal seine Webseite, dort findest du auch andere tolle Beispiele zu statistischen Analysen.

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

longnose_tbl <- read_csv2("data/longnose.csv") |> 
  select(-stream)
```

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-pois-longnose
#| tbl-cap: Auszug aus dem Daten zu den langnasigen Hechten.
#| column: page

rbind(head(longnose_tbl, n = 4),
      rep("...", times = ncol(longnose_tbl)),
      tail(longnose_tbl, n = 4)) |> 
  kable(align = "c", "pipe")
```

Im Folgenden werden wir die Daten nur für das Fitten eines Modells verwenden. In den anderen oben genannten Kapiteln nutzen wir die Daten dann anders. In @fig-pois-model-longnose sehen wir nochmal die Verteilung der Anzahl der Hechte in den Flüssen.

```{r}
#| echo: true
#| message: false
#| label: fig-pois-model-longnose
#| fig-align: center
#| fig-height: 5
#| fig-width: 5
#| fig-cap: "Histogramm der Verteilung der Hechte in den beobachteten Flüssen."

ggplot(longnose_tbl, aes(longnose)) +
  theme_minimal() +
  geom_histogram()

```

## Fit des Modells

In diesem Abschnitt wollen wir verschiedene Modelle für Zähldaten schätzen. Die Poissonverteilung hat keinen eignen Parameter für die Streung wie die Normalverteilung. Die Poissonverteilung ist mit $\mathcal{Pois}(\lambda)$ definiert und hat somit die Eigenschaft das die Varianz eins zu eins mit dem Mittelwert $\lambda$ der Poissonverteilung ansteigt. Es kann aber sein, dass wir in den Daten nicht diesen ein zu eins Zusammenhang von Mittelwert und Varianz vrliegen haben. Häufig ist die Varianz viel größer und steigt schneller an. Wenn die Varianz in Wirklichkeit sehr viel größer ist, dann würden wir die Varianz in unseren Modell unterschätzen.

-   Ein klassisches Poissonmodell `glm(..., familiy = poisson)` mit der Annahme keiner Overdisperison.
-   Ein Quasi-Poissonmodell `glm(..., family = quasipoisson)` mit der Möglichkeit der Berücksichtigung einer Overdispersion.
-   Ein negative Binomialmodell `glm.nb(...)` ebenfalls mit der Berücksichtigung einer Overdispersion.

Beginnen wollen wir aber mit einer klassischen Poissonregression ohne die Annahme von einer Overdispersion in den Daten. Wir nutzen dafür die Funktion `glm()` und spezifizieren die Verteilungsfamilie als `poisson`. Wir nehmen wieder alle Variablen in das Modell auf der rechten Seite des `~`. Auf der linken Seite des `~` kommt dann unser Outcome `longnose` was die Anzahl an Hechten erhält.

[Hier gibt es nur die Kurzfassung der *link*-Funktion. @dormann2013parametrische liefert hierzu in Kapitel 7.1.3 nochmal ein Einführung in das Thema.]{.aside}

Wir müssen für die Possionregression noch beachten, dass die Zähldaten von $0$ bis $+\infty$ laufen. Damit wir normalverteilte Residuen erhalten und einen lineren Zusammenhang, werden wir das Modell auf dem $\log$-scale fitten. Das heißt, wir werden den Zusammenhang von $y$ und $x$ logarithmieren. Wichtig ist hierbei der *Zusammenhang*. Wir transformieren nicht einfach $y$ und lassen den Rest unberührt. Das führt dazu, dass wir am Ende die Koeffizienten der Poissonregression exponieren müssen. Das können die gängigen Funktionen, wir müssen das Exponieren aber aktiv durchführen. Deshalb hier schon mal erwähnt.

```{r}
poisson_fit <- glm(longnose ~ area + do2 + maxdepth + no3 + so4 + temp,
                    longnose_tbl, family = poisson)
```

Wir schauen uns die Ausgabe des Modells einmal mit der `summary()` Funktion an, da wir hier einmal händisch schauen wollen, ob eine Overdispersion vorliegt. Sonst könnten wir auch die Funktion `model_parameters()` nehmen. Die nutzen wir später für die Interpretation des Modells, hier wollen wir erstmal sehen, ob alles geklappt hat.

```{r}
poisson_fit |> summary()
```

Wir schauen in die Summary-Ausgabe des Poissonmodells und sehen, dass dort steht, dass `Dispersion parameter for poisson family taken to be 1`. Wir modellieren also einen eins zu eins Zusammenhang von Mittelwert und Varianz. Wenn dieser Zusammenhang nicht in unseren Daten existiert, dann haben wir eine Overdispersion vorliegen.

Wir können die Overdispersion mit abschätzen indem wir die `Residual deviance` durch die Freiheitsgrade der `Residual deviance` teilen. Daher erhalten wir eine Overdispersion von $\cfrac{1590.04}{61} \approx 26.1$. Damit haben wir eine eindeutige Overdispersion vorliegen. Damit steigt die Varianz in einem Verhältnis von ca. 1 zu 26. Wir können auch die Funktion `check_overdispersion()` aus dem R Paket `{performance}` nutzen um die Overdispersion zu berechnen. Die Funktion kann das schneller und ist auch in der Abfolge einer Analyse besser geeignet.

```{r}
poisson_fit |> check_overdispersion()
```

Wenn wir Overdispersion vorliegen haben und damit die Varianz zu niedrig schätzen, dann erhalten wir viel mehr signifikante Ergebnisse als es in den Daten zu erwarten wäre. Schauen wir uns nochmal die Parameter der Poissonverteilung und die $p$-Werte einmal an.

```{r}
#| message: false

poisson_fit |> model_parameters()
```

In der Spalte `p` finden wir die $p$-Werte für alle Variablen. Wir sehen, dass fast alle Variablen signifikant sind und das wir eine sehr niedrige Varianz in der Spalte `SE` sehen. Das heißt unser geschätzer Fehler ist sehr gering. Das ahnten wir ja schon, immerhin haben wir eine Overdisperson vorliegen. Das Modell ist somit falsch. Wir müssen uns ein neues Modell suchen, was Overdispersion berückscihtigen und modellieren kann.

Die Quasi-Poisson Verteilung hat einen zusätzlichen, unabhänigen Parameter um die Varianz der Verteilung zu schätzen. Daher können wir die Overdispersion mit einer Quasi-Poisson Verteilung berückscihtigen. Wir können eine Quasi-Poisson Verteilung auch mit der Funktion `glm()` schätzen nur müssen wir als Verteilungsfamilie `quasipoisson` angeben.

```{r}
quasipoisson_fit <- glm(longnose ~ area + do2 + maxdepth + no3 + so4 + temp,
                        data = longnose_tbl, family = quasipoisson)
```

Nach dem Modellti können wir nochmal in der `summary()` Funktion schauen, ob wir die Overdispersion richtig berücksichtigt haben.

```{r}
quasipoisson_fit |> summary()
```

An der Zeile `Dispersion parameter for quasipoisson family taken to be 29.403319` in der Summary-Ausgabe sehen wir, dass das Modell der Quasi-Possion Verteilung die Overdispersion korrekt berücksichtigt hat. Wir können uns nun einmal die Modellparameter anschauen. Die Interpretation machen wir am Ende des Kapitels.

```{r}
#| message: false

quasipoisson_fit |> model_parameters()
```

Jetzt sieht unser Modell und die $p$-Werte zusammen mit dem Standardfehler `SE` schon sehr viel besser aus. Wir können also diesem Modell erstmal von der Seite der Overdispersion vertrauen.

Am Ende wollen wir nochmal das Modell mit der negativen Binomialverteilung rechnen. Die negativen Binomialverteilung erlaubt auch eine Unabhängigkeit von dem Mittelwert zu der Varianz. Wir können hier auch für die Overdispersion adjustieren. Wir rechnen die negativen Binomialregression mit der Funktion `glm.nb()` aus dem R Paket `{MASS}`. Wir müssen keine Verteilungsfamilie angeben, die Funktion `glm.nb()` kann nur die negative Binomialverteilung modellieren.

```{r}
negativebinomial_fit <- glm.nb(longnose ~ area + do2 + maxdepth + no3 + so4 + temp,
                               data = longnose_tbl)
```

Auch hier schauen wir mit der Funktion `summary()` einmal, ob die Overdisprsion richtig geschätzt wurde oder ob hier auch eine Unterschätzung des Zusammenhangs des Mittelwerts und der Varianz vorliegt.

```{r}
negativebinomial_fit |> summary()
```

Auch hier sehen wir, dass die Overdispersion mit dem Parameter $\theta$ berücksichtigt wird. Wir können die Zahl $1.67$ nicht direkt mit der Overdispersion aus einer Poissonregression verglechen, aber wir sehen dass das Verhältnis von `Residual deviance` zu den Freiheitsgraden mit $\cfrac{73.65}{61} \approx 1.20$ fast bei 1:1 liegt. Wir könnten also auch eine negative Binomialverteilung für das Modellieren nutzen.

```{r}
#| message: false

negativebinomial_fit |> model_parameters()
```

::: column-margin
Wie immer gibt es reichlich Tipps & Tricks welches Modell du nun nehmen solltest. [How to deal with overdispersion in Poisson regression: quasi-likelihood, negative binomial GLM, or subject-level random effect?](https://stats.stackexchange.com/questions/201903/how-to-deal-with-overdispersion-in-poisson-regression-quasi-likelihood-negativ/332250#332250) und das Tutorial [Modeling Count Data](https://online.stat.psu.edu/stat504/lesson/9/9.2-0). Auch ich mus immer wieder schauen, was am besten konkret in der Anwendung passen könnte und würde.
:::

Welches Modell nun das beste Modell ist, ist schwer zu sagen. Wenn du Overdisperion vorliegen hast, dann ist natürlich nur das Quasi-Poissonmodell oder das negative Binomialmodell möglich. Welche der beiden dann das bessere ist, hängt wieder von der Fragestellung ab. Allgemein gesprochen ist das Quasi-Poissonmodell *besser* wenn dich die Zusammenhänge von $y$ zu $x$ am meisten interessieren. Und das ist in unserem Fall hier die Sachlage. Daher gehen wir mit den Quasi-Poissonmdell dann weiter.

## Performance des Modells

In diesem kurzen Abschnitt wollen wir uns einmal anschauen, ob das Modell neben der Overdispersion auch sonst aus statistischer Sicht in Ordnung ist. Wir wollen ja mit dem Modell aus dem Fit `quasipoisson_fit` weitermachen. Also schauen wir uns einmal das pseudo-$R^2$ für die Poissonregression an. Da wir es mit einem GLM zu tun haben, ist das $R^2$ mit Vorsicht zu genießen. In einer Gaussianregression können wir das $R^2$ als Anteil der erklärten Varianz durch das Modell interpretieren. Im Falle von GLM's müssen wir hier vorsichtiger sein. In GLM's gibt es ja keine Varianz sondern eine Deviance.

```{r}
r2_efron(quasipoisson_fit)
```

Mit einem pseudo-$R^2$ von $0.33$ erklären wir ca. 33% der Varianz in der Anzahl der Hechte. Das ist zwar keine super gute Zahl, aber dafür, dass wir nur eine handvoll von Parametern erfasst haben, ist es dann auch wieder nicht so schlecht. Die Anzahl an Hechten wird sicherlich an ganz vielen Parametern hängen, wir konnten immerhin einige wichtige Stellschrauben vermutlich finden.

In @fig-pois-model-check schauen wir uns nochmal die Daten in den Modelgüteplots an. Wir sehen vorallem, dass wir vielelicht doch einen Ausreißer mit der Beobachtung 17 vorliegen haben. Auch ist der Fit nicht so super, wie wir an dem QQ-Plot sehen. Die Beobachtungen fallen in dem QQ-Plot nicht alle auf eine Linie. Auch sehen wir dieses Muster in dem Residualplot. Hiererwarten wir eine gerade blaue Linie und auch hier haben wir eventuell Ausreißer mit in den Daten.

```{r}
#| echo: true
#| message: false
#| label: fig-pois-model-check
#| fig-align: center
#| fig-height: 8
#| fig-width: 8
#| fig-cap: "Ausgabe ausgewählter Modelgüteplots der Funktion `check_model()`."

check_model(quasipoisson_fit, colors = cbbPalette[6:8], 
            check = c("qq", "outliers", "pp_check", "homogeneity")) 

```

## Interpretation des Modells

Um die Effektschätzer einer Poissonregression oder aber einer Quasipoisson-Regression interpretieren zu können müssen wir uns einmal einen Beispieldatensatz mit bekannten Effekten zwischen den Gruppen bauen. Im Folgenden bauen wir uns einen Datensatz mit zwei Gruppen. Einmal einer Kontrollgruppe mit einer mittleren Anzahl an $15$ und einer Behandlungsgruppe mit einer um $\beta_1 = 10$ höheren Anzahl. Wir haben also in der Kontrolle im Mittel eine Anzahl von $15$ und in der Behandlungsgruppe eine mittlere Anzahl von $25$.

```{r}
sample_size <- 100
longnose_small_tbl <- tibble(grp = rep(c(0, 1), each = sample_size),
                             count = 15 + 10 * grp + rnorm(2 * sample_size, 0, 1)) |>
  mutate(count = round(count),
         grp = factor(grp, labels = c("ctrl", "trt")))
```

In @tbl-pois-longnose-small sehen wir nochmal die Daten als Ausschnitt dargestellt.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-pois-longnose-small
#| tbl-cap: How much is the fish? Der Datensatz über $n = 1000$ Beobachtungen an dem wir überlegen wollen wie wir die Effektschätzer einer Poissonregression zu interpretieren haben. 

longnose_small_raw_tbl <- longnose_small_tbl |> 
  mutate(grp = as.character(grp))

rbind(head(longnose_small_raw_tbl, n = 4),
      rep("...", times = ncol(longnose_small_raw_tbl)),
      tail(longnose_small_raw_tbl, n = 4)) |> 
  kable(align = "c", "pipe")
```

Da sich die Tabelle schlecht liest hier nochmal der Boxplot in @fig-pois-model-small. Wir sehen den Grupenunterschied von $10$ sowie die unterschiedlichen mittleren Anzahlen für die Kontrolle und die Behandlung.

```{r}
#| echo: true
#| message: false
#| label: fig-pois-model-small
#| fig-align: center
#| fig-height: 5
#| fig-width: 5
#| fig-cap: How much is the fish? Der Boxplot über $n = 1000$ Beobachtungen an dem wir überlegen wollen wie wir die Effektschätzer einer Poissonregression zu interpretieren haben. 
#| fig-subcap: 
#|   - "Verteilung der Werte als Boxplot."
#|   - "Verteilung der Werte als Densityplot."
#| layout-nrow: 2

ggplot(longnose_small_tbl, aes(x = grp, y = count, fill = grp)) +
  theme_minimal() +
  geom_boxplot() +
  theme(legend.position = "none") +
  scale_fill_okabeito() 

ggplot(data = longnose_small_tbl, aes(x = count, fill = grp)) +
  theme_minimal() +
  geom_density(alpha = 0.75) +
  labs(x = "", y = "", fill = "Gruppe") +
  scale_fill_okabeito() +
  scale_x_continuous(breaks = seq(10, 30, by = 5), limits = c(10, 30)) 

```

Jetzt fitten wir einmal das simple Poissonmodell mit der Anzahl als Outcome und der Gruppe mit den zwei Leveln als $x$. Wir pipen dann das Ergebnis des Fittes gleich in die Funktion `model_parameters()` weiter um die Ergebnisse des Modellierens zu erhalten.

```{r}
#| message: false
#| warning: false
glm(count ~ grp, data = longnose_small_tbl, family = poisson) |>
  model_parameters(exponentiate = TRUE)
```

Als erstes fällt auf, dass wir die Ausgabe des Modells exponieren müssen. Um einen linearen Zusamenhang hinzukriegen bedient sich die Poissonregression den Trick, das der Zusammenhang zwischen dem $y$ und dem $x$ *transformiert* wird. Wir rechnen unsere Regression nicht auf den echten Daten sondern auf dem $\log$-scale. Daher müssen wir die Koeffizienten der Poissonregression wieder zurücktransfomieren, wenn wir die Koeffizienten interpretieren wollen. Das können wir mit der Option `exponentiate = TRUE` durchführen.

Gut soweit, aber was heißen den jetzt die Zahlen? Wir haben einen Intercept von $14.99$ das entspricht der mittleren Anzahl in der Kontrollgruppe. Und was sagt jetzt die $1.67$ vom Level `trt` des Faktors `grp`? Wenn wir $14.99 \cdot 1.67$ rechnen, dann erhalten wir als Ergebnis $25.03$, also die mittlere Anzahl in der Behandlungsgruppe. Was sagt uns das jetzt aus? Wir erhalten aus der Poissonregression eine Wahrscheinlichkeit oder aber ein Risk Ratio. Wir können sagen, dass die Anzahl in der Behandlungsgruppe $1.67$-mal so groß ist wie in der Kontrollgruppe.

Schauen wir uns nochmal das volle Modell an und interpretieren die Effekte der einzelnen Variablen.

```{r}
#| message: false

quasipoisson_fit |> 
  model_parameters(exponentiate = TRUE) 
```

So schön auch die Funktion `model_parameters()` ist, so haben wir aber hier das Problem, dass wir den Effekt von `area` nicht mehr richtig sehen. Wir kriegen hier eine zu starke Rundung auf zwei Nachkommastellen. Wir nutzen jetzt mal die Funktion `tidy()` um hier Abhilfe zu leisten. Ich muss hier noch die Spalte `estimate` mit `num(..., digits = 5)` anpassen, damit du in der Ausgabe auf der Webseite auch die Nachkommastellen siehst.

```{r}
#| message: false

quasipoisson_fit |> 
  tidy(exponentiate = TRUE, digits = 5) |> 
  select(term, estimate, p.value) |> 
  mutate(p.value = pvalue(p.value),
         estimate = num(estimate, digits = 5))
```

Schauen wir uns die Effekte der Poissonregression einmal an und versuchen die Ergebnisse zu interpretieren. Dabei ist wichtig sich zu erinnern, dass kein Effekt eine 1 bedeutet. Wir schauen hier auf einen Faktor. Wenn wir eine Anzahl mal Faktor 1 nehmen, dann ändert sich nichts an der Anzahl.

-   `(Intercept)` beschreibt den Intercept der Poissonregression. Wenn wir mehr als eine simple Regression vorliegen haben, wie in diesem Fall, dann ist der Intercept schwer zu interpretieren. Wir konzentrieren uns auf die Effekte der anderen Variablen.
-   `area`, beschreibt den Effekt der Fläche. Steigt die Fläche um ein Quadratmeter an, so erhöht sich die *Anzahl an Fischen* um den $1.00001$. Daher würde man hier eher sagen, erhöht sich die Fläche um jeweils 1000qm so erhöht sich die *Anzahl an Fischen* um den Faktor $1.1$. Dann haben wir auch einen besser zu interpretierenden Effektschätzer. Die Signifikanz bleibt hier davon unbetroffen.
-   `do2`, beschreibt den Partzialdruck des Sauerstoffs. Steigt dieser um eine Einheit an, so sehen wie eine Erhöhung der *Anzahl an Fischen* um den Faktor $1.25$. Der Effekt ist gerade *nicht* signifikant.
-   `maxdepth`, beschreibt die maximale Tiefe. Je tiefer ein Fluss, desto mehr Hechte werden wir beobachten. Der Effekt von $1.01$ pro Meter Tiefe ist signifikant.
-   `no3`, beschreibt den Anteil an Nitrat in den Flüssen. Je mehr Nitrat desto signifiant mehr Hechte werden wir beobachten. Hier steigt der Faktor auch um $1.20$.
-   `so4`, beschreibt den Schwefelgehalt und mit steigenden Schwefelgehalt nimmt die *Anzahl an Fischen* leicht ab. Der Effekt ist aber überhaupt nicht signifikant.
-   `temp`, beschreibt die Temperatur der Flüsse. Mit steigender Temperatur erwarten wir mehr Hechte zu beobachten. Der Effekt von $1.08$ Fischen pro Grad Erhöhung ist signifikant.

Was nehmen wir aus der Poissonregression zu den langnasigen Hechten mit? Zum einen haben die Fläche, die Tiefe und der Nitratgehalt einen signifikanten Einfluss auf die Anzahl an Hechten. Auch führt eine höhere Temperatur zu mehr gefundenen Hechten. Die erhöhte Temperatur steht etwas im Widerspuch zu dem Sauerstoffpartizaldruck. Denn je höher die Temperatur desto weniger Sauerstoff wird in dem Wasser gelöst sein. Auch scheint die Oberfläche mit der Tiefe korreliert. Allgemein scheinen Hechte große Flüsse zu mögen. Hier bietet sich also noch eine Variablenselektion oder eine Untersuchung auf Ausreißer an um solche Effekte nochmal gesondert zu betrachten.

## Zeroinflation

So eine Poissonregression hat schon einiges an Eigenheiten. Neben dem Problem der Overdispersion gibt es aber noch eine weitere Sache, die wir beachten müssen. Wir können bei einer Poissonregression auch eine Zeroinflation vorliegen haben. Das heißt, wir beobachten viel mehr Nullen in den Daen, als wir aus der Poissonverteilung erwarten würden. Es gibt also einen biologischen oder künstlichn Prozess, der uns Nullen produziert. Häufig wissen wir nicht, ob wir den Prozess, der uns die Nullen in den Daten produziert, auch abbilden. Das heißt, es kann sein, dass wir einfach nichts Zählen, weil dort nichts ist oder aber es gibt dafür einen Grund. Diesen Grund müssten wir dann irgendwie in unseren Daten erfasst haben, aber meistens haben wir das nicht.

Schauen wir usn dafür einmal ein Datenbeispiel von Eidechsen in der Lüneburgerheide an. Wir haben Eidechsen `lizard` in zwei verschiedenen Habitaten `grp` gezählt. Einmal, ob die Eidechsen eher im offenen Gelände oder eher im bedeckten Gelände zu finden waren. Im Weiteren haben wir geschaut, ob der Boden keinen Regen erhalten hatte, trocken war oder gar feucht. Mit trocken ist hier eine gewisse Restfeuchte gemeint. Am Ende haben wir noch bestimmt, ob wir eher nah an einer Siedlung waren oder eher weiter entfernt. Du kannst dir den Daten satz in der Datei `lizards.csv` nochmal anschauen. In @tbl-lizard-data sind die Daten nochmal dargestellt.

```{r}
#| message: false
#| warning: false
#| echo: false
#| tbl-cap: "Ausschnitt aus den Eidechsendaten für die zwei Habitate unter verschiedenen Feuchtigkeitsbedingungen und Nähe zur nächsten Siedlung."
#| label: tbl-lizard-data

set.seed(20210061)

lizard_zero_tbl <- expand_grid(grp = 1:2,
                               rain = 1:3,
                               pop = 1:2,
                               rep = 1:5) |> 
  mutate(lizard = 0 + 1.5 * grp - 2 * rain + 1.2 * pop + rnorm(n(), 0, 2),
         lizard = ifelse(round(lizard) < 0, 0, round(lizard)),
         grp = factor(grp, labels = c("open", "cover")),
         rain = factor(rain, labels = c("no", "dry", "wet")),
         pop = factor(pop, labels = c("near", "far"))) |> 
  select(-rep)

write_csv2(lizard_zero_tbl, "data/lizards.csv")

lizard_zero_tbl |> head(n = 7) |> kable(align = "c", "pipe")

```

In @fig-pois-lizard-data sehen wir die Zähldaten der Eidechsen nochmal als Histogramm dargestellt. Wenn wir an einem Punkt keine Eidechsen gefunden haben, dann haben wir keine fehlenden Werte eingetragen, sondern eben, dass wir keine Eidechsen gezählt haben. Wir sehen das wir sehr viele Nullen in unseren Daten haben. Ein Indiz für eine Inflation an Nullen oder eben einer Zeroinflation.

```{r}
#| echo: true
#| message: false
#| label: fig-pois-lizard-data
#| fig-align: center
#| fig-height: 5
#| fig-width: 5
#| fig-cap: "Histogramm der Verteilung der Hechte in den beobachteten Flüssen."

ggplot(lizard_zero_tbl, aes(lizard)) +
  theme_minimal() +
  geom_histogram() +
  labs(x = "Anzahl der gefundenen Eidechsen", y = "Anzahl") +
  scale_x_continuous(breaks = 0:7)

```

Um zu überprüfen, ob wir eine Zeroinflation in den Daten vorliegen haben, werden wir erstmal eine ganz normale Poissonregression auf den Daten rechnen. Wir ignorieren auch eine potenzielle Overdispersion. Das schauen wir uns dann in den Daten später nochmal an.

```{r}
#| message: false
#| warning: false
lizard_fit <- glm(lizard ~ grp + rain + pop, data = lizard_zero_tbl,
                  family = poisson)
```

Wie immer nutzen wir die Funktion `model_parameters()` um uns die exponierten Koeffizienten aus dem Modell wiedergeben zu lassen. Das Modell dient uns jetzt nur als Ausgangsmodell und wir werden das Poissonmodell jetzt nicht weiter tiefer verwenden.

```{r}
#| message: false
#| warning: false
lizard_fit |> model_parameters(exponentiate = TRUE)
```

Wir sehen, dass wir in der Variable `rain` eine starke Reduzierung der Anzahl an Eidechsen sehen. Vielleicht ist dies eine Variable, die zu viele Nullen produziert. Auch hat die Variable `pop`, die für die Nähe an einer Siedlung kodiert, einen starken positiven Effekt auf unsere Anzahl an Eidechsen. Hier wollen wir also einmal auf eine Zeroinflation überprüfen. Wir nutzen dazu die Funktion `check_zeroinflation()` aus dem R Paket `{performance}`. Die Funktion läuft nur auf einem Modellfit.

```{r}
#| message: true
#| warning: false
check_zeroinflation(lizard_fit)
```

Die Funktion gibt uns wieder, dass wir vermutlich eine Zeroinflation vorliegen haben. Das können wir aber Modellieren. Um eine Zeroinflation *ohne Overdispersion* zu modellieren nutzen wir die Funktion `zeroinfl()` aus dem R Paket `{pscl}`. Der erste Teil der Funktion ist leicht erkläret. Wir bauen uns wieder unswer Model zusammen, was wir fitten wollen. Dann kommt aber ein `|` und mit diesem Symbol `|` definieren wir, ob wir wissen, woher die Nullen kommen oder aber ob wir die Nullen mit einem zufälligen Prozess modellieren wollen.

Wenn wir das Modell in der Form `y ~ f1 + f2 | 1` schreiben, dann nehmen wir an, dass das Übermaß an Nullen in unseren Daten rein zufällig entstanden sind. Wir haben keine Spalte in de Daten, die uns eine Erklärung für die zusätzlichen Nullen liefern würde.

Wir können auch `y ~ f1 + f2 | x3` schreiben. Dann haben wir eine Variable `x3` in den Daten von der wir glauben ein Großteil der Nullen herrührt. Wir könnten also in unseren Daten annehmen, dass wir den Überschuss an Nullen durch den Regen erhalten haben und damit über die Spalte `rain` den Exzess an Nullen modellieren.

Man sollte immer mit dem einfachsten Modell anfangen, deshalb werden wir jetzt einmal ein Modell fitten, dass annimmt, dass die Nullen durch einen uns unbekannten Zufallsprozess entstanden sind.

```{r}
#| message: false
#| warning: false
lizard_zero_infl_intercept_fit <- zeroinfl(lizard ~ grp + pop + rain | 1, 
                                           data = lizard_zero_tbl) 
```

Wir schauen uns das Modell dann wieder einmal an und sehen eine Zweiteilung der Ausgabe. In dem oberen Teil der Ausgabe wird unsere Anzahl an Eidechsen modelliert. In dem unteren Teil wird der Anteil der Nullen in den Daten modelliert. Daher können wir über Variablen in dem `Zero-Inflation` Block keine Aussagen über die Anzahl an Eidechsen treffen. Variablen tauchen nämlich nur in einem der beiden Blöcke auf.

```{r}
#| message: false
#| warning: false
lizard_zero_infl_intercept_fit |> 
  model_parameters(exponentiate = TRUE)
```

Als erstes beobachten wir einen größeren Effekt der Variable `grp`. Das ist schon mal ein spannender Effekt. An der Signifikanz hat scih nicht viel geändert. Wir werden am Ende des Kapitels einmal alle Modell für die Modellierung der Zeroinflation vergleichen.

Nun könnte es auch sein, dass der Effekt der vielen Nullen in unserer Variable `rain` verborgen liegt. Wenn es also regnet, dann werden wir viel weniger Eidechsen beoabchten. Nehmen wir also `rain` als ursächliche Variable mit in das Modell für die Zeroinflation.

```{r}
#| message: false
#| warning: false
lizard_zero_infl_rain_fit <- zeroinfl(lizard ~ grp + pop | rain, 
                                      data = lizard_zero_tbl)
```

Wieder schauen wir uns einmal die Ausgabe des Modells einmal genauer an.

```{r}
#| message: false
#| warning: false
lizard_zero_infl_rain_fit |> model_parameters(exponentiate = TRUE)
```

Es ändert sich einiges. Zum einen erfahren wir, dass der Regen anscheined doch viele Nullen in den Daten produziert. Wir haben ein extrem hohes $OR$ für die Variable `rain`. Die Signifikanz ist jedoch eher gering. Wir haben nämlich auch eine sehr hohe Streuung mit den großen $OR$ vorliegen. Au der anderen Seite verlieren wir jetzt auch die Signifikanz von unseren Habitaten und dem Standort der Population. Nur so mäßig super dieses Modell.

Wir können jetzt natürlich auch noch den Standort der Population mit in den Prozess für die Entstehung der Nullen hineinnehmen. Wir schauen uns dieses Modell aber nicht mehr im Detail an, sondern dann nur im Vergleich zu den anderen Modellen.

```{r}
#| message: false
#| warning: false
lizard_zero_infl_rain_pop_fit <- zeroinfl(lizard ~ grp | rain + pop, 
                                          data = lizard_zero_tbl)
```

Die Gefahr besteht immer, das man sich an die Wand modelliert und vor lauter Modellen die Übersicht verliert. Neben der Zeroinflation müssen wir ja auch schauen, ob wir eventuell eine Overdispersion in den Daten vorliegen haben. Wenn das der Fall ist, dann müsen wir nochmal überlegen, was wir dann machen. Wir testen nun auf Ovrdisprsion in unserem *ursprünglichen* Poissonmodell mit der Funktion `check_overdispersion()`.

```{r}
#| message: true
#| warning: false
check_overdispersion(lizard_fit)
```

Tja, und so erfahren wir, dass wir auch noch Overdispersion in unseren Daten vorliegen haben. Wir müsen also beides Modellieren. Einmal modellieren wir die Zeroinflation und einmal die Overdispersion. Wir können beides in einem negativen binominalen Modell fitten. Auch hier hilft die Funktion `zeroinfl()` mit der Option `dist = negbin`. Mit der Option geben wir an, dass wir eine negative binominal Verteilungsfamilie wählen. Damit können wir dann auch die Ovrdispersion in unseren Daten modellieren.

```{r}
#| message: false
#| warning: false
lizard_zero_nb_intercept_fit <- zeroinfl(lizard ~ grp + rain + pop | 1, 
                                         dist = "negbin", data = lizard_zero_tbl)
```

Dann schauen wir usn einmal das Modell an. Zum einen sehen wir, dass der Effekt ähnlich groß ist, wie bei dem Intercept Modell der Funktion `zeroinfl`. Auch bleiben die Signifikanzen ähnlich.

```{r}
#| message: false
#| warning: false
lizard_zero_nb_intercept_fit |> model_parameters(exponentiate = TRUE)
```

Nun haben wir vier Modelle geschätzt und wolen jetzt wissen, was ist das beste Modell. Dafür hilft usn dann eine Gegenüberstellung der Modelle mit der Funktion `modelsummary()`. Wir könnten die Modelle auch gegeneinander statistsich Testen, aber hier behalten wir uns einmal den beschreibenden Vergleich vor. In @tbl-model-comp-zero sehen wir einmal die vier Modelle nebeneinander gestellt. Für eine bessere Übrsicht, habe ich aus allen Modellen den Intercept entfernt.

```{r}
#| message: false
#| echo: true
#| tbl-cap: "Modellvergleich mit den vier Modellen. Wir schauen in wie weit sich die Koeffizienten und Modelgüten für die einzelnen Modelle im direkten Vergleich zum vollen Modell verändert haben."
#| label: tbl-model-comp-zero

modelsummary(lst("ZeroInfl Intercept" = lizard_zero_infl_intercept_fit,
                 "ZeroInfl rain" = lizard_zero_infl_rain_fit,
                 "ZeroInfl rain+pop" = lizard_zero_infl_rain_pop_fit,
                 "NegBinom intercept" = lizard_zero_nb_intercept_fit),
             statistic = c("conf.int",
                           "s.e. = {std.error}", 
                           "t = {statistic}",
                           "p = {p.value}"),
             coef_omit = "Intercept", 
             exponentiate = TRUE)
```

Die beiden Intercept Modelle haben die kleinsten $AIC$-Werte der vier Modelle. Darüber hinaus haben dann beide Modelle auch die höchsten $R^2_{adj}$ Werte. Beide Modelle erklären also im Verhältnis viel Varianz mit 58.5%. Auch ist der $RMSE$ Wert als Fehler bei beiden Modellen am kleinsten. Damit haben wir die Qual der Wahl, welches Modell wir nehmen. Ich würde das negative binominal Modell nehmen. Wir haben ins unseren Daten vermutlich eine Zeroinflation sowie eine Overdispersion vorliegen. Daher bietest es sich an, beides in einer negativen binominalen Regression zu berücksichtigen. Zwar sind die beiden Intercept Modelle in diesem Beispielfall von den Koeffizienten fast numerisch gleich, aber das hat eher mit dem reduzierten Beispiel zu tun, als mit dem eigentlichen Modell. In unserem Fall ist die Overdispersion nicht so extrem.

Wie sehe den unser negative binominal Modell aus, wenn wir mit dem Modell einmal die zu erwartenden Eidechsen vorhersagen würden? Auch das kann helfen um abzuschätzen, ob das Modelle einigermaßen funktioniert hat. Wir haben ja hier den Vorteil, dass wir nur mit kategorialen Daten arbeiten. Wir haben keine kontiniuerlichen Variablen vorliegen und darüber hinaus auch nicht so viele Variablen insgesamt.

Daher bauen wir uns mit `expand_grid()` erstmal einen Datensatz, der nur aus den Faktorkombinationen besteht. Wir haben also nur eine Beobachtung je Faktorkombination. Danach nutzen wir die Daten einmal in der Funktion `predict()` um uns die vorhergesagten Eidechsen nach dem gefitten Modell wiedergeben zu lassen.

```{r}
newdata_tbl <- expand_grid(grp = factor(1:2, labels = c("open", "cover")),
                           rain = factor(1:3, labels = c("no", "dry", "wet")),
                           pop = factor(1:2, labels = c("near", "far")))

pred_lizards <- predict(lizard_zero_nb_intercept_fit, newdata = newdata_tbl) 
  
newdata_tbl <- newdata_tbl |> 
  mutate(lizard = pred_lizards)
```

Nachdem wir in dem Datensatz `newdata_tbl` nun die vorhergesagten Eidechsen haben, können wir uns jetzt in der @fig-pred-lizard die Zusammenhänge nochmal anschauen.

```{r}
#| echo: true
#| message: false
#| label: fig-pred-lizard
#| fig-align: center
#| fig-height: 5
#| fig-width: 10
#| fig-cap: "Scatterplot der vorhergesagten Eidechsen in den Habitaten (`grp`), der Feuchtigkeit des Bodens nach Regen und dem Abstand zur nächsten Ortschaft."


ggplot(newdata_tbl, aes(x = rain, y = lizard, colour = grp, group = grp)) +
  theme_minimal() +
  geom_point() +
  geom_line() +
  facet_wrap(~ pop) +
  labs(x = "Feuchtigkeit nach Regen", y = "Anzahl der gezählten Eidechsen",
       color = "Gruppe") +
  scale_color_okabeito()
```

Wir erkennen, dass mit der Erhöhung der Feuchtigkeit die Anzahl an aufgefundenen Eidechsen sinkt. Der Effekt ist nicht mehr so stark, wenn es schon einmal geregnet hat. Ebenso macht es einen Unterschied, ob wir nahe einer Siedlung sind oder nicht. Grundsätzlich finden wir immer mehr Eidechsen in geschützten Habitaten als in offenen Habitaten.

## Gruppenvergleich {#sec-mult-comp-pois-reg}

Häufig ist es ja so, dass wir das Modell für die Poisson Regression nur schätzen um dann einen Gruppenvergleich zu rechnen. Das heißt, dass es uns interessiert, ob es einen Unterschied zwischen den Leveln eines Faktors gegeben dem Outcome $y$ gibt. Da wir hier in unserem Beispiel zu den Flusshechten keine Gruppe drin haben, zeige ich dir das Prinzip einmal an einem ausgedachten Datensatz. Wir nehmen hier einen Datensatz zu Katzen-, Hunde- und Fuchsflöhen. Dabei erstellen wir uns einen Datensatz mit mittleren Anzahlen an Flöhen pro Tierart. Ich habe jetzt eine mittlere Flohanzahl von $10$ Flöhen bei den Katzen, eine mittlere Anzahl von $30$ Flöhen bei den Hunden und eine mittlere Anzahl von $10$ Flöhen bei den Füchsen gewählt. Wir generieren uns jeweils $20$ Beobachtungen je Tierart. Damit haben wir dann einen Datensatz zusammen, den wir nutzen können um einmal die Ergebnisse eines Gruppenvergleiches mit Zähldaten zu verstehen.

```{r}
set.seed(20231202)
n_rep <- 20
flea_count_tbl <- tibble(animal = gl(3, n_rep, labels = c("cat", "dog", "fox")),
                         count = round(c(rnorm(n_rep, 20, 1), 
                                         rnorm(n_rep, 30, 1), 
                                         rnorm(n_rep, 10, 1))))
```

Wenn du gerade hierher gesprungen bist, nochmal das simple Modell für unseren Gruppenvergleich unter einer Poisson Regression. Wir haben hier nur einen Faktor `animal` mit in dem Modell. Am Ende des Abschnitts findest du dann noch ein Beispiel mit zwei Faktoren zu Thripsen auf Apfelbäumen. Wir modellieren hier eigentlich nicht die Zähldaten an sich, sondern die Raten. Du wirst das aber gleich bei der Ausgabe von `emmeans()` sehen.

```{r}
pois_fit <- glm(count ~ animal, data = flea_count_tbl, family = "poisson") 
```

Eigentlich ist es recht einfach, wie wir anfangen. Wir rechnen jetzt als erstes die ANOVA. Hier müssen wir dann einmal den Test angeben, der gerechnet werden soll um die p-Werte zu erhalten. Dann nutze ich noch die Funktion `model_parameters()` um eine schönere Ausgabe zu erhalten.

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

pois_fit |> 
  anova(test = "Chisq") |> 
  model_parameters(drop = "NULL")
```

Im Folgenden nutzen wir das R Paket `{emmeans}` wie folgt. Wenn wir die mittleren Anzahlen benötigen, dann müssen wir die Option `type = "response"` verwenden. Sonst würdest du Werte auf der Link-Skala wiederbekommen, die dir hier nicht helfen. Ich nutze später die `emmeans` Ausgabe um ein Säulendiagramm mit den mittleren Anzahlen zu erstellen. Wir du gleich sehen wirst, werden wir aber nicht die Anzahlen vergleichen sondern die *Raten* in den jeweiligen Gruppen.

```{r}
emm_obj <- pois_fit |> 
  emmeans(~ animal, type = "response")
emm_obj
```

Jetzt steht hier zwar `rate` aber was ist das denn nun? Dafür berechnen wir mal die Mittelwerte der Anzahlen für die drei Tierarten über die Funktion `summarise()`. Dann vergleichen wir einmal die Ausgaben und schauen, was wir in `emmeans()` berechnet haben.

```{r}
flea_count_tbl |> 
  group_by(animal) |> 
  summarise(mean_count = mean(count))
```

Wunderbar, es sind die Mittelwerte der Anzahlen, die in der Spalte `rate` in der Ausgabe von `emmeans()` stehen. Damit können wir dann arbeiten und die Ausgabe nutzen um ein Säulendigramm zu erstellen. Bitte schaue für die Umsetzung in das Beispiel am Ende des Abschnitts. Da zeige ich dir dann wie du ein Säulendigramm zu den Anzahlen von Thripsen auf Apfelbäumen erstellst.

Wenn du die mittleren Anzahlen pro Tierart jetzt für jede Tierart testen willst, dann erhälst du aber keinen Mittelwertsunterschied der Anzahlen, sondern eine *Rate*. Rechnen wir aber erstmal den paarweisen Vergleich und schauen uns dann an, was wir erhalten haben. Wir kriegen am Ende wiederum einen Vergleich von Wahrscheinlichkeiten, aber eben hier von *Raten*. Wie immer, es kommt darauf an, was du willst. Hier einmal die paarweisen Vergleiche, darunter dann gleich das *compact letter display*.

```{r}
emm_obj |> 
  pairs(adjust = "bonferroni")
```

Okay, jetzt haben wir statt Raten also eine *Ratio*, also ein Verhältnis von Raten. Das schauen wir uns doch gleich mal in dem Beispiel an. Wir rechnen also einmal die mittlere Anzahl für jede Tierart durch die entsprechende andere mittlere Anzahl der anderen Tierart. Wir erhalten also folgende *Ratios* der mittleren Anzahlen.

$$
\cfrac{cat}{dog} = \cfrac{20.1}{30.1} = 0.67 
$$

$$
\cfrac{cat}{fox} = \cfrac{20.1}{10.2} = 1.97 
$$

$$
\cfrac{dog}{fox} = \cfrac{30.1}{10.2} = 2.95 
$$

Wie interpretieren wir jetzt die *Ratios*? Eigentlich ist das intuitiver als man auf den ersten Blick denkt. Wir haben zum Beispiel fast doppelt so viele Flöhe bei Katzen wie bei Füchsen. Die Anzahl von Hundeflöhen ist fast dreimal so hoch wie die Anzahl an Fuchsflöhen. Wir haben fast nur halb so viele Flöhe auf Katzen gefunden als bei den Hunden. Damit können wir dann schon arbeiten und eine Aussage treffen.

Und fast am Ende können wir uns auch das *compact letter display* erstellen. Auch hier nutzen wir wieder die Funktion `cld()` aus dem R Paket `{multcomp}`. Hier erhälst du dann die Information über die mittlere Anzahl der jeweiligen Tierarten und ob sich die mittlere Anzahl unterscheidet. Ich nutze dann die Ausgabe von `emmeans()` um mir dann direkt das Säulendiagramm mit den Fehlerbalken und dem *compact letter display* zu erstellen. Mehr dazu dann im Kasten weiter unten zu dem Beispiel zu Thripsenanzahl auf Apfelbäumen.

```{r}
#| message: false
#| warning: false
emm_obj |>
  cld(Letters = letters, adjust = "none")
```

Auch hier sehen wir, dass sich alle drei Gruppen signifikant unterschieden, keine der Tierarten teilt sich einen Buchstaben, so dass wir hier von einem Unterschied zwischen den mittleren Anzahlen an Flöhen auf den drei Tierarten ausgehen können.

::: callout-note
## Offset in einer Poisson Regression

Bei einem geplanten Experiment zählen wir meist in festgelegten Zeiteinheiten oder aber unser $n$ ist vorgeben. Jetzt kann es aber sein, dass wir nicht mit festen Zeitintervallen oder aber einer Fallzahl zählen. Beispielsweise können wir die Anzahl der Baumarten in einem Wald zählen: Ereignisse (eng. *event*) wären Baumbeobachtungen, die Exposition (eng. *exposure*) wäre eine Flächeneinheit und die Rate wäre die Anzahl der Arten pro Flächeneinheit. Wir können die Sterberaten in geografischen Gebieten als die Anzahl der Todesfälle geteilt durch die Personenjahre modellieren. Allgemeiner ausgedrückt können Ereignisraten als Ereignisse pro Zeiteinheit berechnet werden, wobei das Beobachtungsfenster für jede Einheit variieren kann. In diesen Beispielen ist die Exposition jeweils eine Einheit Fläche, Personenjahre und Zeiteinheit. So sollten beispielsweise sechs Fälle innerhalb eines Jahres nicht den gleichen Wert haben wie sechs Fälle innerhalb von 10 Jahren. In der Poisson-Regression wird dies als *Offset* behandelt. Mehr dazu dann gerne im Wikipediartikel zu ["Exposure" and offset](https://en.wikipedia.org/wiki/Poisson_regression#.22Exposure.22_and_offset) oder aber die Antwort [When to use an offset in a Poisson regression?](https://stats.stackexchange.com/questions/11182/when-to-use-an-offset-in-a-poisson-regression).

In `{emmeans}` können wir [Modelle mit Offset](https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html#offsets) relativ einfach modellieren. Dazu nutzen wir die Funktion `offset()` in unserem `glm()` Modell. Aber zuerst einmal ein Spieldatensatz. Wir haben wir unseren Faktor Tierart und dann aber noch das Alter als Faktor mit zwei Stufen. Wir haben jetzt aber nicht immer die gleiche Anzahl an Hunden und Katzen sowie Füchsen ausgezählt. Unsere Anzahl an Flöhen basiert also auf einer anderen $n$ Anzahl. Die Anzahl $n$ modellieren wir dann als `offset()`.

```{r}
toy_count_tbl <- tibble(n = c(500, 1200, 100, 400, 500, 300),
                        animal = factor(rep(1:3,2), labels = c("cat","dog","fox")),
                        age = gl(2, 3),
                        count = c(42, 37, 1, 101, 73, 14))
```

Wir nehmen die Anzahl an ausgezählten Tieren dann als Offset mit Logarithmus in das Modell. Deshalb steht dann da auch `log(n)` im Offset.

```{r}
toy_count_fit <- glm(count ~ animal + age + offset(log(n)), 
                     data = toy_count_tbl, family = "poisson")
```

Dann können wir uns einmal das Modell anschauen, wie `{emmeans}` den Vergleich rechnen würde.

```{r}
ref_grid(toy_count_fit)
```

Der Offset ist also grob $log(500) \approx 5.97$. Das stimmt numerisch nicht hundertprozentig, hat aber noch mit einer internen Korrektur zu tun, die uns hier aber nicht weiter interessiert. Das Modell können wir dann wie gewohnt in `emmeans()` stecken und dann weiter rechnen. In dem Modell wird dann berücksichtigt, dass sich die Anzahlen für das *Exposure*, also dem Nenner der Raten, unterscheiden.
:::

Damit sind wir einmal mit unserem Gruppenvergleich für die Poisson Regression auf Zähldaten durch. In dem Kapitel zu den [Multiple Vergleichen oder Post-hoc Tests](#sec-posthoc) findest du dann noch mehr Inspirationen für die Nutzung von `{emmeans}`. Hier war es dann die Anwendung auf Zähldaten zusammen mit einem Faktor. Wenn du dir das Ganze nochmal an einem Beispiel für zwei Faktoren anschauen möchtest, dann findest du im folgenden Kasten ein Beispiel für die Auswertung von Thripsen auf Apfelbäumen nach Gabe verschiedener Dosen eines Insektizids und Zeitpunkten.

::: callout-tip
## Anwendungsbeispiel: Zweifaktorieller Gruppenvergleich für Thripsenbefall

Im folgenden Beispiel schauen wir uns nochmal ein praktische Auswertung von einem agrarwissenschaftlichen Beispiel mit jungen Apfelbäumen an. Wir haben uns in diesem Experiment verschiedene Dosen `trt` von einem Insektizid aufgebracht sowie verschiedene Startanzahlen von Raubmilben als biologische Alternative untersucht. Dann haben wir noch fünf Zeitpunkte bestimmt, an denen wir die Anzahl an Thripsen auf den Blättern gezählt haben. Wir haben nicht die Blätter per se gezählt sondern Fallen waagerecht aufgestellt. Dann haben wir geschaut, wie viele Thripsen wir über `above` und unter `below` von den Fallen gefunden haben. In unserem Fall beschränken wir uns auf die obere Anzahl an Thripsen und schauen uns auch nur die Behandlung mit dem Insektizid an.

```{r}
insects_tbl <- read_excel("data/insects_count.xlsx") |> 
  mutate(timepoint = factor(timepoint, labels = c("1 Tag", "4 Tag", "7 Tag", "11 Tag", "14 Tag")),
         rep = as_factor(rep),
         trt = as_factor(trt)) |>
  select(timepoint, trt, thripse = thripse_above) |> 
  filter(trt %in% c("10ml", "30ml", "60ml"))
```

Dann können wir auch schon die Poisson Regression mit `glm()` rechnen. Auch hier wieder darauf achten, dass wir dann als Option `family = poisson` oder `family = quasipoisson` wählen. Es hängt jetzt davon ab, ob du in deinen Daten Overdispersion vorliegen hast oder nicht. In den beiden folgenden Tabs, rechne ich dann mal beide Modelle.

::: panel-tabset
## `family = poisson`

Als Erstes rechnen wir eine normale Poisson Regression und schauen einmal, ob wir Overdispersion vorliegen haben. Wenn wir Overdispersion vorliegen haben, dann können wir keine Poisson Regression rechnen, sondern müssen auf eine Quasipoisson Regression ausweichen. Das ist aber sehr einfach, wie du im anderen Tab sehen wirst.

```{r}
insects_poisson_fit <- glm(thripse ~ trt + timepoint + trt:timepoint, 
                           data = insects_tbl, 
                           family = poisson) 
```

Bevor wir uns das Modell mit `summary()` überhaupt anschauen, wollen wir erstmal überprüfen, ob wir überhaupt Overdispersion vorliegen haben. Wenn ja, dann können wir uns die `summary()` hier gleich sparen. Also einmal geguckt, was die Overdispersion macht.

```{r}
insects_poisson_fit |> check_overdispersion()
```

Wir haben sehr starke Overdispersion vorliegen und gehen daher rüber in den anderen Tab und rechnen eine Quasipoisson Regression. Nur wenn du keine Overdispersion vorliegen hast, dann kannst du eine eine Poisson Regression rechnen.

## `family = quasipoisson`

Entweder hast du in deinen Daten eine Overdispersion gefunden oder aber du meinst, es wäre besser gleich eine Quasipoisson zu rechnen. Beides ist vollkommen in Ordnung. Ich rechne meistens immer eine Quasipoisson und schaue dann nur, ob die Overdispersion sehr groß war. In den seltensten Fällen hast du eine Overdispersion vorliegen, die eher klein ist. Daher mache ich erst die Lösung und schaue, ob das Problem dann da war.

```{r}
insects_quasipoisson_fit <- glm(thripse ~ trt + timepoint + trt:timepoint, 
                                data = insects_tbl, 
                                family = quasipoisson) 
```

Du kannst in der `summary()` Ausgabe direkt sehen, ob du Overdispersion vorliegen hast. Du musst nur relativ weit unten schauen, was zu dem `Dispersion parameter` in den Klammern geschrieben ist. Wenn da eine Zahl größer als 1 drin steht, dann hast du Overdispersion.

```{r}
insects_quasipoisson_fit |> 
  summary()
```

Wir haben hier auf jeden Fall Overdispersion vorliegen. Daher nutze ich dann auch das Modell hier mit der Annahme an eine Quasipoissonverteilung. Dann stimmt es auch mit unseren Varianzen und wir produzieren nicht zufällig zu viele signifikante Ergebnisse, die es dann gar nicht gibt.
:::

Ich habe mich gerade in den obigen Tabs für eine Quasipoisson Regression entschieden, da wir Overdispersion vorliegen haben. Damit mache ich dann mit dem `insects_quasipoisson_fit` Modell weiter. In den beiden folgenden Tabs findest du dann einmal das Ergebnis für die ANOVA und einmal für den Gruppenvergleich mit dem R Paket `{emmeans}`. Bitte beachte, dass die ANOVA für ein `glm()`-Objekt nicht ganz gleich wie für ein `lm()`-Objekt ist. Du kannst aber die ANOVA erstmal ganz normal interpretieren, nur haben wir hier nicht die Möglichkeit ein $\eta^2$ zu bestimmen. Dann nutzen wir `{emmeans}` für den Gruppenvergleich. Nochmal, weil wir Overdispersion festgestellt haben, nutzen wir das Objekt `insects_quasipoisson_fit` mit der Berücksichtigung der Overdispersion.

::: panel-tabset
## ANOVA mit `anova()`

Wir rechnen hier einmal die ANOVA und nutzen den $\mathcal{X}^2$-Test für die Ermittelung der p-Werte. Wir müssen hier einen Test auswählen, da per Standardeinstellung kein Test gerechnet wird. Wir machen dann die Ausageb nochmal schöner und fertig sind wir.

```{r}
insects_quasipoisson_fit |> 
  anova(test = "Chisq") |> 
  model_parameters(drop = "NULL")
```

Wir sehen, dass der Effekt für die Behandlung signifikant ist, jedoch die Zeit und die Interaktion keinen signifikanten Einfluss haben. Wir haben aber also keine Interaktion vorliegen. Daher können wir dann die Analyse gemeinsam über alle Zeitpunkte rechnen.

## Gruppenvergleich mit `emmeans()`

Im Folgenden rechnen wir einmal über alle Faktorkombinationen von `trt` und `timepoint` einen Gruppenvergleich. Dafür nutzen wir die Opition `trt * timepoint`. Wenn du die Analyse *getrennt* für die Zeitpunkte durchführen willst, dann nutze die Option `trt | timepoint`. Wir wollen die Wahrscheinlichkeiten für das Auftreten einer Beschädigung von wiedergegeben bekommen, deshalb die Option `regrid = "response`. Dann adjustieren wir noch nach Bonferroni und sind fertig.

```{r}
emm_obj <- insects_quasipoisson_fit |> 
  emmeans(~ trt * timepoint, regrid = "response") |>
  cld(Letters = letters, adjust = "bonferroni")
emm_obj
```

Das `emm_obj` Objekt werden wir dann gleich einmal in `{ggplot}` visualisieren. Die `rate` stellt die mittlere Anzahl an Thripsen je Faktorkombination dar. Dann können wir auch das *compact letter display* anhand der Abbildung interpretieren.
:::

In der @fig-log-mod-pois-insect siehst du das Ergebnis der Auswertung in einem Säulendiagramm. Hier unbedingt `SE` als den Standardfehler für die Fehlerbalken nutzen, da wir sonst Fehlerbalken größer und kleiner als $0$ erhalten, wenn wir die Standardabweichung nutzen würden. Das ist in unserem Fall nicht so das Problem, aber wenn du eher kleine Anzahlen zählst, kann das schnell zu Werten kleiner Null führen. Wir sehen einen klaren Effekt der Behandlung `60ml`. Die Zeit hat keinen Effekt, was ja schon aus der ANOVA klar war, die Säulen sehen für jeden Zeitpunkt vollkommen gleich aus. Gut etwas Unterschied ist ja immer.

```{r}
#| echo: true
#| warning: false
#| message: false
#| label: fig-log-mod-pois-insect
#| fig-align: center
#| fig-height: 4
#| fig-width: 6
#| fig-cap: "Säulendigramm der mitleren Zahl der Thripsen aus einer Poisson Regression. Das `glm()`-Modell berechnet die mittlere Anzahl in jeder Faktorkombination. Das *compact letter display* wird dann in `{emmeans}` generiert."

emm_obj |> 
  as_tibble() |> 
  ggplot(aes(x = timepoint, y = rate, fill = trt)) +
  theme_minimal() + 
  labs(y = "Mittlere Anzahl an Thripsen", x = "Messzeitpunkte der Zählungen",
       fill = "Dosis") +
  geom_bar(stat = "identity", 
           position = position_dodge(width = 0.9, preserve = "single")) +
  geom_text(aes(label = .group, y = rate + SE + 0.01),  
            position = position_dodge(width = 0.9), vjust = -0.25) +
  geom_errorbar(aes(ymin = rate-SE, ymax = rate+SE),
                width = 0.2,  
                position = position_dodge(width = 0.9, preserve = "single")) +
  scale_fill_okabeito()
```
:::

## Referenzen {.unnumbered}