-
Notifications
You must be signed in to change notification settings - Fork 0
/
2f_splitplot.Rmd
122 lines (96 loc) · 8.07 KB
/
2f_splitplot.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 - Split-Plot Design"
output:
html_document:
includes:
after_body: footer.html
---
```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
library(data.table)
library(desplot) # plotte Feldplan
library(RColorBrewer) # Farben für Feldplan
load("D:/RKurse/Dokumentation/crashcouRse/datasets/rice split plot.rda")
# Mehr Infos
# Split Plot:
# Bailey (2008) 8.3 + 8.4
# Dean & Voss (1998) 2.44 + 19
# Quinn & Keough (2002) 11.1
```
# Datensatz
```{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
library(RColorBrewer) # bessere Farben
```
In diesem Experiment wurden die Auswirkungen von 6 Stickstoffdüngern (N) auf den Ertrag von 4 Reissorten (G) untersucht. Demnach gab es zwei Behandlungsfaktoren mit insgesamt 6x4=24 Behandlungsfaktorstufenkombinationen. Diese 24 Kombinationen wurden 3 mal in vollständigen Blöcken wiederholt. Innerhalb der vollständigen Blöcke wurden die Kombinationen allerdings nicht wie bei einem RCBD vollständig randomisiert. Grund hierfür war, dass es beim Anlegen des Feldversuchs und speziell beim Ausbringen der Dünger praktischer ist, wenn die Parzellen mit demselben Dünger in einer Fahrtrichtung liegen. Deshalb wurden die vollständigen Blöcke erst in 6 Spalten (=*main plots*) unterteilt, auf welche die 6 N-Stufen randomisert wurden. Erst dann, wurden die 4 Sorten den 4 Zeilen innerhalb jeder Spalte (=*sub plots*) zufällig zugeordnet. Solch eine Versuchsanlage nennt man **split-plot design**. [Hier mehr Infos zu Versuchsdesigns](appendix_designs.html)
```{r, echo=F, fig.height = 3, fig.width = 5, fig.align = "center"}
desplot(data=splitplot, form= rep ~ col+row|rep, flip=T,
text=G, cex=1, col=N,
col.regions=c("grey100", "grey90", "grey80"),
col.text=brewer.pal(6, "Dark2"),
main="Feldplan", key.cex = 0.6,
out1=row, out1.gpar=list(col="grey50", lwd=1, lty=1),
out2=col, out2.gpar=list(col="grey50", lwd=1, lty=1))
splitplot <- splitplot[,c("N", "G", "rep", "mainplot", "yield")]
```
<div class = "row"> <div class = "col-md-6">
```{r}
print(splitplot, nrows=10)
```
</div> <div class = "col-md-6">
```{r, eval=FALSE}
str(splitplot)
```
```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
str(splitplot, 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 Stufenkombinationen der zwei Behandlungsfaktoren. Mittels des `col=` statements der `boxplot()` Funktion können wir angeben welche Füllfarben die Boxen haben sollen. Wir wählen dieselben Farben wie im Feldplan - [mehr Infos dazu hier](datr_desplot.html).
```{r}
boxplot(yield ~ N*G, col=brewer.pal(6, "Dark2"), data=splitplot, las=2) # las=2 dreht Achsenbeschriftung
```
Es fällt auf, dass sich das Muster der Werte der verschiedenen N-Stufen für Sorte D deutlich von denen der anderen 3 Sorten unterscheiden. Dies deutet auf Wechselwirkungen zwischen den beiden Behandlungsfaktoren hin.
# Schließende Statistik
## Gemischtes Modell
Wir können uns nun entschließen die Daten mittels eines Modells zu analysieren. Der Ertrag ist unsere metrische Zielvariable. 'G' und 'N' sind beides qualitative Behandlungsfaktoren. Um eine möglich Wechselwirkung der beiden zu modellieren, ergibt sich der *treament part* des Modell als `G + N + G*N`.
Für die vollständigen Blöcke haben wir den festen Effekt 'rep' im Modell. Desweiteren müssen wir zufällige Effekte für die main plots ins Modell aufnehmen, da diese als zusätzliche Randomisationseinheit fungierten. Die andere, kleinste Randomisationseinheit des sub plots (=Parzelle) wird durch den üblichen Fehlerterm abgebildet. Der *design part* des Modells ist also `rep + (1|rep:mainplot)`.
Somit haben wir gleichzeitig feste und zufällige Effekte im Modell und demnach ein gemischtes Modell, welches wir wieder mit `lmer()` aus dem `lme4` package modellieren.
```{r}
mod <- lmer(yield ~ G + N + G:N + rep + (1|rep:mainplot), data=splitplot)
```
Zunächst sollten nun die Residuenplots evaluiert werden - beispielsweise mit `plot(mod)` und `qqnorm(resid(mod)); qqline(resid(mod))`. Erst danach ist eine Varianzanalyse zulässig.
> Die Spalte 'mainplot' ist in diesem Beispiel übrigens prinzipiell überflüssig. Das liegt daran, dass man die 3x6 main plots anstatt mit 'rep:mainplot' auch mit 'rep:N' hätte modellieren können, da ein main plot immer einer N-Stufe entspricht (siehe Feldplan oben). Um es aber leichter verständlich zu machen, nutzen wir hier eine extra 'mainplot' Spalte. Siehe außerdem die Anmerkung aus [dem vorigen Beispiel](1f_alpha.html).
## Varianzanalyse
```{r, warning=FALSE, message=FALSE, error=FALSE}
library(car)
Anova(mod, test.statistic="F", type="III")
```
Da dies Ergebnisse einer mehrfaktorielle Varianzanalyse sind, betrachten zuallererst den F-Test für die Wechselwirkung (mehr Infos [hier](appendix_interaktionen.html)). Der signifikante p-Wert von `G:N` deutet darauf hin, dass die beiden Behandlungsfaktoren nicht rein additiv wirken, sondern wechselwirken. Dies bestätigt die Vermutung aus dem Boxplot der Rohdaten oben. Das bedeutet außerdem, dass wir die adjustierten Mittelwerte der Wechselwirkungseffekte berechnen wollen und *nicht* die der beiden Haupteffekte.
## Multipler Mittelwertvergleich
Mit `emmeans()` können wir auch adjustierte Mittelwerte für Wechselwirkungseffekte berechnen. Das geht entweder mit `pairwise ~ N : G` oder `pairwise ~ N | G`. Im beiden Fällen werden wie erwartet Mittelwerte für alle Kombinationen errechnet. Der Unterschied zwischen den beiden Befehlen besteht darin welche Differenzen zwischen den Mittelwerten berechnet werden. Im ersten Fall (`N:G`) werden **alle** Mittelwerte miteinander verglichen. Im zweiten Fall (`N|G`), werden nur N-Stufen innerhalb einer Sorte miteinander verglichen. Natürlich können 'N' und 'G' hier auch ausgetauscht werden. Es kommt auf die Versuchsfrage an, welche Mittelwertvergleiche von Interesse sind. Wir wollen uns hier auf die Vergleiche der N-Mittelwerte getrennt pro Sorte fokussieren.
```{r}
means <- emmeans(mod, pairwise ~ N | G, adjust="tukey") # adjust="none" für t-test
means <- CLD(means$emmeans, details=TRUE, Letters=letters)
means$emmeans # N-Mittelwerte pro Sorte
means$comparisons # N-Mittelwertdifferenzen pro Sorte
```
# Ergebnisaufbereitung
Wir wollen die Ergebnisse abschließend in einem Balkendiagramm darstellen. Weil wir aber die N-Vergleiche isoliert pro Sorte konzentrieren wollen, können wir uns entscheiden ein Balkendiagramm pro Sorte zu erstellen. In solch einer Situation kann in der `ggplot()` Funktion das `facet_wrap()` statement genutzt werden.
Außerdem wollen wir die Balken wieder mit denselben Farben versehen, die die N-Stufen bereits im Feldplan und dem Boxplot hatten. Dazu fügen wir zuerst dem `geom_bar()` statement `fill=N` hinzu. Nun bekommt jede N-Stufe eine andere Farbe. Um nicht die default Farben zu nehmen, sondern selbst zu bestimmen welche Farben es sein sollen, schreiben wir zusätzlich noch `scale_fill_manual(values=...)` in die Funktion. Um die dann automatisch generierte Legende zu unterdrücken (weil sie in diesem Fall keine zusätzlichen Informationen bringt), fügen wir außerdem `guides(fill=FALSE)` hinzu.
```{r}
means.plot <- as.data.table(means$emmeans) # korrektes Format für ggplot
print(means.plot, nrows=2)
ggplot(data=means.plot, aes(x=N)) +
geom_bar(aes(y=emmean, fill=N), stat="identity", width=0.8) +
scale_fill_manual(values=brewer.pal(6, "Dark2")) + guides(fill=FALSE) +
geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=0.4) +
geom_text(aes(y=emmean+1500, label=.group)) +
facet_wrap(~G) + # ein plot pro Sorte
theme_bw() +
labs(main="Adj. N-Mittelwerte pro Sorte",
y="Adjustierter Ertragsmittelwert ± 95% Konfidenzintervall", x="N-Dünger",
caption="Getrennt für jede Sorte gilt: Mittelwerte, die mit einem gemeinsamen Buchstaben versehen sind,\n sind laut Tukey-test nicht signifikant voneinander verschieden.")
```