-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path1f_rcbd.Rmd
121 lines (101 loc) · 6.32 KB
/
1f_rcbd.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
---
title: "1 Beh.faktor - RCBD"
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
rcbd <- data.table(john.alpha)
rcbd$row <- rep(rep(4:1),18)
rcbd$col <- sort(rep(seq(18:1),4))
levels(rcbd$rep) <- c("Rep1", "Rep2", "Rep3")
```
```{r, warning=FALSE, message=FALSE, error=FALSE}
library(data.table) # bessere Datenmanipulation
library(ggplot2); library(ggfortify) # bessere Plots
library(emmeans) # adjustierte Mittelwerte
```
Dieses Beispiel ist dem Beispiel "1f_crd" sehr ähnlich und baut darauf auf. In diesem Feldversuch wurde allerdings der Ertrag von 24 Sorten gemessen, welche auch in 3 Wiederholungen, aber als **randomisierte vollständige Blockanlage** (=**R**andomized **C**omplete **B**lock **D**esign) angelegt wurden. Demnach gibt es nun eine zusätzliche Spalte im Datensatz, die Informationen zu den vollständigen Blöcken enthält. Diese Spalte heißt hier "rep" (*replicate*), da vollständige Blöcke auch oft *Wiederholungen* genannt werden. [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"}
colors <- c(brewer.pal(8,"Blues"),
brewer.pal(8,"Oranges"),
brewer.pal(8,"YlGn"))
desplot(data=rcbd, form= gen ~ col+row, col.regions=colors,
text=gen, cex=0.4, shorten="no",
main="Feldplan",
xlab="", ylab="", show.key=F, out1=rep)
desplot(data=rcbd, form= rep ~ col+row|rep, col.regions=brewer.pal(3,"Reds"),
text=gen, cex=0.4, shorten="no",
main="Vollständige Blöcke",
xlab="", ylab="", show.key=F)
rcbd <- rcbd[,c("gen", "rep", "yield")]
```
</div> <div class = "col-md-6">
```{r}
print(rcbd, 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 dem Beispiel "1f_crd" erstellen wir diesmal auch einen Boxplot für die Blöcke.
<div class = "row"> <div class = "col-md-6">
```{r, eval=FALSE}
str(rcbd)
```
```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
str(rcbd, width=40, strict.width="cut")
```
```{r}
plot(y=rcbd$yield, x=rcbd$gen, las=2) # las=2 dreht Achsenbeschriftung
```
</div> <div class = "col-md-6">
```{r}
summary(rcbd)
plot(y=rcbd$yield, x=rcbd$rep)
```
</div> </div>
# Schließende Statistik
## Lineares Modell
Wir können uns nun entschließen die Daten mittels eines linearen Modells zu analysieren. Der Ertrag ist unsere metrische Zielvariable. 'Sorte' ist ein qualitativer Faktor. Außerdem haben wir den nun noch den qualitativen Faktor 'rep' im Modell.
```{r}
mod <- lm(yield ~ gen + rep, data=rcbd)
```
Zunächst sollten nun die Residuenplots (z.b. mit `autoplot(mod)`) evaluiert werden, was hier aber übersprungen wird. Erst dann ist eine Varianzanalyse zulässig.
## Varianzanalyse
```{r}
anova(mod)
```
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.
```{r}
means <- emmeans(mod, pairwise ~ gen, adjust="tukey") # adjust="none" für t-test
as.data.table(means$emmeans)[1:6,] # 6 der 24 Mittelwerte
```
```{r}
as.data.table(means$contrasts)[1:6,] # 6 der 276 Differenzen
```
Im Vergleich zum Beispiel "1f_crd" gibt es hier deutlich mehr Differenzen/paarweise Vergleiche. Das liegt daran, dass wir dort nur 3 Sortenmittelwerte miteinander verglichen haben, während es hier 24 sind. Mit steigender Sortenanzahl *n* steigt die Anzahl aller möglichen Sortenvergleiche *n(n-1)/2* sehr schnell an. Beim Betrachten der p-Werte fällt außerdem wieder auf, dass nicht alle Differenzen signifikant sind. Bei dieser größeren Anzahl Sorten wird deutlich wie hilfreich die Buchstabendarstellung ist, welche wir erneut mit dem `CLD()` statement erzeugen können.
```{r}
means <- CLD(means$emmeans, details=TRUE, Letters=letters)
as.data.table(means$emmeans)[1:6,]
```
# Ergebnisaufbereitung
Erneut wollen wir die Ergebnisse abschließend in einem Balkendiagramm darstellen. Diesmal nutzen wir aber mehr `ggplot`-Funktionen als im Beispiel "1f_crd". So kann man generell mit dem `theme()` statement [sehr viele](https://ggplot2.tidyverse.org/reference/theme.html) Dinge am Aussehen eines ggplots ändern. Um aber nicht jede Linie, Füllfarbe usw. einzeln bearbeiten zu müssen, kann man auch [vorgefertigte Themes](https://ggplot2.tidyverse.org/reference/ggtheme.html) nutzen, wie wir hier das `theme_bw()`. Desweiteren lassen wir die Fehlerbalken diesmal nicht ± Standardfehler abbilden, sondern die Spanne des 95% Konfidenzintervall, die auch schon im `emmeans()` statement berechnet wurden.
Außerdem ordnen wir die Sorten auf der x-Achse diesmal nicht nach Namen, sondern nach ihrem adjustierten Mittelwert. mit dem `reorder()` statement in der ersten Zeile. Schließlich lassen wir die Information zur Interpretation der Buchstabendarstellung diesmal nicht als `caption=` unter der Grafik erscheinen, sondern fügen mittles `annotate()`eine Textbox innerhalb des Diagramms ein. Der Text wird automatisch auf der angegebenen x/y-Koordinate zentriert. Im Satz fügen wir einen Absatz mit der Zeichenfolge *\\n* ein.
```{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.")
```