1f_alpha.Rmd

---
title: "1 Beh.faktor - Alpha-Design"
output: 
  html_document:
    includes:
      after_body: footer.html 
---

# Datensatz

```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
library(agridat) # agrarwissenschaftliche Beispieldatensätze
library(data.table)
library(desplot) # plotte Feldplan
library(RColorBrewer) # Farben für Feldplan
alpha <- data.table(john.alpha)
alpha$row <- rep(rep(4:1),18)
alpha$col <- sort(rep(seq(18:1),4))
names(alpha)[3] <- "inc.block"
levels(alpha$rep) <- c("Rep1", "Rep2", "Rep3")
```

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

Dieses Beispiel ist den Beispielen "1f_crd" und vor allem "1F_rcbd" sehr ähnlich und baut darauf auf. Tatsächlich handelt es sich um den gleichen Datensatz wie in "1F_rcbd". Es ist ein Beispieldatensatz mit dem Namen `john.alpha` aus dem `agridat` package.

Dieses Experiment wurde als ein $\alpha$-lattice angelegt. Dies wurde im letzten Beispiel verschwiegen und stattdessen behauptet, es wäre eine randomisierte vollständige Blockanlage. In der Tat besitzt das eigentliche Experiment dieselben vollständigen Blöcke, allerdings gibt es innerhalb dieser noch jeweils 6 unvollständige Blöcke. Die 'rep' Spalte steht also weiterhin für die vollständigen Blöcke, zusätzlich gibt es jedoch noch eine 'inc.block' (*incomplete block*) Spalte für die unvollständigen Blöcke. [Hier mehr Infos zu Versuchsdesigns](appendix_designs.html)

<div class = "row"> <div class = "col-md-6">
```{r, echo=F, fig.height = 2, fig.width = 4, fig.align = "center"}
colors24 <- c(brewer.pal(8,"Blues"),
              brewer.pal(8,"Oranges"),
              brewer.pal(8,"YlGn"))
colors18 <- c(brewer.pal(6,"Blues"),
              brewer.pal(6,"Oranges"),
              brewer.pal(6,"YlGn"))
alpha$real.inc.block <- paste0(alpha$rep, alpha$inc.block)
desplot(data=alpha, form= gen ~ col+row, col.regions=colors24,
        text=gen, cex=0.4, shorten="no",
        main="Feldplan",
        xlab="", ylab="", show.key=F, out1=rep,
        out2=inc.block, out2.gpar=list(col="black", lwd=1, lty=1))
desplot(data=alpha, form=real.inc.block~col+row|rep,
        col.regions=colors18,
        text=gen, cex=0.4, shorten="no",
        main="Vollst. und unvollst. Blöcke",
        xlab="", ylab="", show.key=F, 
        out1=real.inc.block, out1.gpar=list(col="black", lwd=1, lty=1))
alpha <- alpha[,c("gen", "rep", "inc.block", "yield")]
```
</div> <div class = "col-md-6">
```{r}
print(alpha, nrows=10)
```
</div> </div>

# Deskriptive Statistik
Erst wollen wir ein Gefühl für den Datensatz bekommen und betrachten einige Kennzahlen zu den Daten, sowie zwei Plots. Im Vergleich zu den vorangegangenen '1wayANOVA' Beispielen erstellen wir diesmal auch einen komplexeren Boxplot via `boxplot()` für die unvollständigem Blöcke.

<div class = "row"> <div class = "col-md-6">
```{r, eval=FALSE}
str(alpha)
plot(y=alpha$yield, x=alpha$gen)
```
```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
str(alpha, width=40, strict.width="cut")
plot(y=alpha$yield, x=alpha$gen)
```
</div> <div class = "col-md-6">
```{r}
summary(alpha)
plot(y=alpha$yield, x=alpha$rep)
```
</div> </div>

```{r}
boxplot(yield ~ rep + inc.block, data=alpha, las=2) #las=2 dreht Achsenbeschriftung
```

# Schließende Statistik {#Anmerkung1}
## Gemischtes Modell
Wir können uns nun entschließen die Daten mittels eines Modells zu analysieren. Der Ertrag ist unsere metrische Zielvariable. 'Sorte' ist ein qualitativer Faktor. Außerdem haben wir den qualitativen Faktor 'rep' im Modell. Im Gegensatz zu den vorangegangenen 1wayANOVA Beispielen müssen wir diesmal auch Effekte für die unvollständigen Blöcke modellieren. Diese sind ebenso qualitativ, sollten aber in diesem Fall als zufällige Effekte genommen werden. Somit haben wir gleichzeitig feste und zufällige Effekte im Modell und demnach ein gemischtes Modell. Die Funktion `lm()` ist nicht in der Lage gemischte Modelle anzupassen, sodass wir in diesem Beispiel zu `lmer()` aus dem `lme4` package wechseln, welches immer zusammen mit dem `lmerTest()` package geladen werden sollte. Die Syntax ist dem von `lm()` recht ähnlich, mit dem Unterschied, dass zufällige Effekte generell wie folgt codiert werden:

```{r}
mod <- lmer(yield ~ gen + rep + (1|rep:inc.block), data=alpha)
```

Zunächst sollten nun die Residuenplots evaluiert werden, was bei einem `lmer()` Objekt nicht mit `autoplot(mod)`, sondern beispielsweise mit `plot(mod)` und `qqnorm(resid(mod)); qqline(resid(mod))` funktioniert. Erst danach ist eine Varianzanalyse zulässig.

> Es fällt auf, dass wir den zufälligen Effekt für die unvollständigen Blöcke zusammen mit 'rep' und getrennt durch einen Doppelpunkt geschrieben haben. Das liegt daran wie die 'inc.block' Spalte in diesem Datensatz angelegt wurde. Sie enthält lediglich die 6 verschiedenen Einträge B1-B6. Es gab allerdings in jedem vollständigem Block 6, also insgesamt 18 unvollständige Blöcke. Ähnlich wie bei dem Boxplot oben muss also die Kombination aus 'rep' und 'inc.block' angegeben werden um im Modell die korrekten 18 zufälligen Effekte für die unvollständigen Blöcke anzupassen. Als Alternative hätte man auch direkt B1-B18 in die 'inc.block' Spalte und demnach `(1|inc.block)` ins Modell schreiben können. In den hier präsentierten Beispielen, werden wir solche Effekte allerdings immer so darstellen wie im Modell oben, da dies verdeutlicht, dass die unvollständigen Blöcke geschachtelt innerhalb der vollständigen Blöcke vorkommen. Hätte man außerdem fälschlicherweise `(1|inc.block)` ins Modell geschrieben obwohl 'inc.block' nur B1-B6 enthält, so wären nur 6 Effekte modelliert worden. B1 z.B. wäre dann ein Effekt über alle drei der jeweils ersten Spalte von links in den drei vollständigen Blöcken. `lmer` würde keine Fehlermeldung ausgeben, sondern mit dem falschen Modell falsche Ergebnisse ausgeben.  

## Varianzanalyse
```{r}
anova(mod)
```
Obwohl der Befehl `anova(mod)` exakt derselbe scheint wie in den vorangegangenen 1wayANOVA Beispielen, wird er nun durch das `lmerTest` package durchgeführt, sodass andere Berechnungen durchgeführt werden, die besser für Varianzanalysen von gemischten Modellen geeignet sind.
Der F-Test für den Faktor 'Sorte' zeigt einen p-Wert < 0.05 und somit signifikante Unterschiede zwischen den Sorten. Demnach wissen wir nun, **dass** es mindestens einen signifikanten Unterschied zwischen den Sorten gibt, aber nicht zwischen **welchen** Sorten. Um dies herauszufinden ist es üblich multiple Mittelwertvergleiche durchzuführen (z.B. t-test oder Tukey-test). 

## Multipler Mittelwertvergleich
Mit `emmeans()` erhalten wir in einem Zug sowowhl die adjustierten Mittelwerte für jede Sorte, als auch die Differenzen zwischen allen Sortenmittelwerten. Es ist sehr praktisch, dass das `emmeans` package nicht nur `lm()` Objekte, sondern auch `lmer()` Objekte, also gemischte Modelle, als Grundlage für die adjustierte Mittelwertschätzung verwenden kann. Auch die Buchstabendarstellung erhalten wir mit dem selben Code wie im letzten Beispiel.

```{r}
means <- emmeans(mod, pairwise ~ gen, adjust="tukey")  # adjust="none" für t-test
means <- CLD(means$emmeans, details=TRUE, Letters=letters)
as.data.table(means$emmeans)[1:6,] # 6 der 24 Mittelwerte 
as.data.table(means$comparisons)[1:6,] # 6 der 276 Differenzen
```


# Ergebnisaufbereitung

Erneut wollen wir die Ergebnisse abschließend in einem Balkendiagramm darstellen. 

```{r}
ggplot(data=means$emmeans, aes(x=reorder(gen, emmean))) +
  geom_bar(aes(y=emmean), stat="identity", width=0.8) +
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=0.4) +
  geom_text(aes(y=(upper.CL)*1.1, label=.group, angle=90)) +
  labs(y="Adjustierter Ertragsmittelwert mit 95% Konfidenzintervall", x="Sorte") +
  theme_bw() +
  annotate(geom="label", y=1, x=12, size=3, color="grey50", fill="white",
  label="Mittelwerte, die mit einem gemeinsamen Buchstaben versehen sind,\n sind laut Tukey-test nicht signifikant voneinander verschieden.")
```