1f_latsq_poisson.Rmd

---
title: "1 Beh.faktor - Lat.Quadrat - Poissondaten"
output: 
  html_document:
    includes:
      after_body: footer.html 
---

# Datensatz

```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
library(data.table)
library(desplot) # plotte Feldplan
library(RColorBrewer) # Farben für Feldplan
trtfarben <- brewer.pal(5, "BuPu")
load("D:/RKurse/Dokumentation/crashcouRse/datasets/bugs poisson.rda")
```

```{r, warning=FALSE, message=FALSE, error=FALSE}
library(data.table) # bessere Datenmanipulation
library(ggplot2); library(ggfortify) # bessere Plots
library(emmeans) # adjustierte Mittelwerte
```

In diesem Beispiel wurden 5 verschiedene Begasungsbehandlungen für  Schnellkäfer (*Elateridae*) mit 5 Wiederholungen getestet. Das Versuchsdesign war ein lateinisches Quadrat - sowohl die Spalten als auch die Zeilen sind also vollständige Blöcke. In jeder Parzelle wurden im Jahr nach der Behandlung die Schnellkäferlarven gezählt. 

<div class = "row"> <div class = "col-md-6">
```{r, echo=F, fig.height = 2, fig.width = 4, fig.align = "center"}
desplot(data=bugs.poiss, form= trt ~ col+row,
        col.regions=trtfarben,
        text=trt, shorten="no", cex=0.8,
        out1=Col, out1.gpar=list(col="black", size=1.5),
        out2=Row, out2.gpar=list(col="black", size=1.5),
        main="Behandlung", show.key=F)
```
</div> <div class = "col-md-6">
```{r, echo=F, fig.height = 2, fig.width = 4, fig.align = "center"}
grays <- colorRampPalette(c("white","#252525"))
desplot(data=bugs.poiss, form= bugs ~ col+row,
        col.regions=grays(18), at=0:18-1,
        text=bugs, shorten="no", cex=0.8,
        main="Anzahl Käferlarven", show.key=F)
bugs <- bugs.poiss[, c("trt", "Row", "Col", "bugs")]
names(bugs) <- c("trt", "row", "col", "bugs")
bugs$row <- as.factor(bugs$row)
bugs$col <- as.factor(bugs$col)
```
</div> </div>

<div class = "row"> <div class = "col-md-6">
```{r}
print(bugs, nrows=10)
```
</div> <div class = "col-md-6">
```{r, eval=FALSE}
str(bugs)
```
```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
str(bugs, width=40, strict.width="cut")
```
</div> </div>

# Deskriptive Statistik

Erst wollen wir ein Gefühl für den Datensatz bekommen und betrachten einen Boxplot für die Anzahl gefundener Larven pro Behandlung. Die Füllfarben die Boxen sind dieselben Farben wie im Feldplan und wurden zuvor von uns im Vektor `trtfarben` gespeichert - [mehr Infos dazu hier](datr_desplot.html). 

```{r, fig.height = 4, fig.width = 4, fig.align = "center"}
boxplot(bugs ~ trt, col=trtfarben, data=bugs, las=2) # las=2 dreht die Achsenbeschriftung
```

# Schließende Statistik
## Generalisiertes Modell

Wie zu erkennen, ist das Auftreten der Larven ein relativ seltenes Ereignis, sodass die meisten Werte nahe 0 liegen. Wir haben also *Zählwerte ohne feste Beschränkung nach oben* und können uns entscheiden eine Poissonverteilung für diese Daten anzunehmen und ein entsprechendes generalisiertes lineares Modell anpassen ([mehr dazu hier](appendix_nichtnormalverteilt.html)). Bezüglich der Effekte im Modell verfahren wir ähnlich einer [einfaktoriellen ANOVA in einem RCBD](1f_rcbd.html), allerdings mit noch einem zweiten vollständigen Block Effekt, sodass wir Effekte für Zeilen und Spalten haben.  

```{r}
mod <- glm(bugs ~ trt + row + col, family=quasipoisson(link="log"), data=bugs)
```

Wie im [vorangegangenen Beispiel](1f_rcbd_binomial.html), legen wir die erwartete Verteileung mit dem `family=` statement fest. Ebenso entscheiden wir uns direkt für `quasipoisson` anstelle der einfacheren `poisson` Verteilung, da diese in gewisser Hinischt flexibler/robuster ist ([mehr dazu hier](intro_glm_carrot.html)). Auch hier ist der Zusatz `(link="log")` eigentlich überflüssig, weil dies die standardmäßige Linkfunktion für poissonverteilte Daten ist. Sie wird hier aber zur Verdeutlichung mit ins Modell geschrieben.

```{r, fig.height = 2, fig.width = 4, fig.align = "center"}
autoplot(mod)[1:2] # Residuenplots
```

## Test der festen Effekte

```{r}
anova(mod, test="Chisq")
car::Anova(mod, test.statistic="Wald")
```

Der p-Wert für den Behandlungseffekt weist in diesem Fall auf unterschiedliche Häufigkeiten der gefundenen Larven hin.

## Multipler Mittelwertvergleich

Mittels des `emmeans()` package lassen wir uns [wie im vorigen Beispiel](1f_rcbd_binomial.html) direkt die rücktransformierten Mittelwerte schätzen. Dazu muss lediglich `type="response"` hinzugefügt werden.

```{r, warning=FALSE, message=FALSE}
means <- emmeans(mod, pairwise ~ trt, type="response") # Mittelwertvergleiche
means <- CLD(means$emmeans, Letters = letters) # Buchstabenddarstellung
means$.group <- gsub(" ", "", means$.group, fixed = TRUE) # Entferne Leerzeichen
means
```

# Ergebnisaufbereitung

```{r}
ggplot() + theme_classic() +
   # Rohdaten (bugs)
   geom_boxplot(data=bugs, aes(x=trt, y=bugs), outlier.shape=NA, width=0.6, fill=trtfarben) +
   geom_jitter(data=bugs, aes(x=trt, y=bugs), width=0.1, shape=1, size=2) +
   # Ergebnisse (means)
   geom_point(data=means, aes(x=as.numeric(trt)+0.4, y=rate), col="red", shape=16, size=2) +
   geom_errorbar(data=means, aes(x=as.numeric(trt)+0.4, ymin=asymp.LCL, ymax=asymp.UCL), col="red", width=0.1) +
   geom_text(data=means, aes(x=trt, y=11, label =.group), col="red") +
  labs(y="Schwarz: Beobachtete Anzahl Larven\nRot: Rücktransformierter, adjustierter Mittelwert für die \nHäufigkeit ± 95% asympt. Konfidenzintervall", x="Behandlung",
       caption="Mittelwerte, die mit einem gemeinsamen Buchstaben versehen sind,\n sind laut Tukey-test auf der Log-Skala nicht signifikant voneinander verschieden.")
```