In diesem Teil wollen wir den Datensatz npk
auswerten:
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.3 ✔ forcats 0.5.2 ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
<- as_tibble(npk)
npk npk
Wie bereits erläutert handelt es sich bei dem Datensatz um ein
typischen vollständigen Versuchsplan aus dem Bereich der
Agrarwissenschaften. Der Datzensatz beschreibt den Einfluss des
Einsatzes von Stickstoff (N
), Phosphat (P
) und
Kalium (K
) auf den Ertrag von Erbsen
(yield
).
Werten Sie den Datzensatz npk
mit einer passenden
ANOVA aus.
Verschaffen Sie sich eine (visuelle/statistische) Übersicht über die Daten.
Prüfen Sie die Voraussetzungen der einfachen ANOVA.
Wenden Sie das passende ANOVA-Modell an.
Interpretieren Sie die Ergebnisse.
npk
mit einer passenden ANOVA aus.# Dimension, Datentypen:
npk
str(npk)
tibble [24 × 5] (S3: tbl_df/tbl/data.frame)
$ block: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
$ N : Factor w/ 2 levels "0","1": 1 2 1 2 2 2 1 1 1 2 ...
$ P : Factor w/ 2 levels "0","1": 2 2 1 1 1 2 1 2 2 2 ...
$ K : Factor w/ 2 levels "0","1": 2 1 1 2 1 2 2 1 1 2 ...
$ yield: num [1:24] 49.5 62.8 46.8 57 59.8 58.5 55.5 56 62.8 55.8 ...
dim(npk)
[1] 24 5
# Statistische Maßzahlen
summary(npk)
block N P K yield
1:4 0:12 0:12 0:12 Min. :44.20
2:4 1:12 1:12 1:12 1st Qu.:49.73
3:4 Median :55.65
4:4 Mean :54.88
5:4 3rd Qu.:58.62
6:4 Max. :69.50
Zuvor hatten wir bereits gesehen, dass der Versuch balanciert ist, d. h. jede Kombination von dem Einsatz von Stickstoff, Phosphat und Kalium liegt genau in dreifacher Wiederholung vor.
An dieser Stelle könnten wir den Ertrag noch zusätzlich visualisieren:
ggplot(npk) +
geom_boxplot(aes(N, yield, color = P, fill = K)) +
labs(y = "Ertrag") +
scale_color_viridis_d()
Ausreißer: In dem vorherigen Boxplots schienen keine Ausreißer vorhanden zu sein.
library(rstatix)
Attache Paket: ‘rstatix’
Das folgende Objekt ist maskiert ‘package:stats’:
filter
%>%
npk group_by(N, P, K) %>%
identify_outliers(yield)
Normalverteilung: Hierzu betrachten wir Histogramme sowie den Shapiro-Wilk-Test:
ggplot(npk, aes(yield, ..density..)) +
geom_histogram(fill = "ForestGreen", alpha = 0.9, bins = 8) +
geom_density(color = "tomato", fill = "tomato3", alpha = 0.2)
Die gesamte Tendenz der Variable yield
könnte
normalverteilt sein.
Einzelnen Gruppen können wir hier nur schwierig betrachten, da der Umfang zu gering ist:
ggplot(npk, aes(yield, ..density.., fill = N, color = N)) +
geom_histogram( alpha = 0.7, bins = 6) +
geom_density(alpha = 0.3) +
facet_grid(P ~ K)
Shapiro-Wilk-Test für alle Gruppen:
shapiro.test(npk$yield)
%>%
npk group_by(N, P, K) %>%
shapiro_test(yield)
Der Shapiro-Wilk-Test für alle Gruppen gemeinsam gibt mit p =
0,87 (> ⍺) an, dass yield
normalverteilt ist.
Wenn wir für jede einzelne Gruppe separat testen, scheint nur der Einsatz von Stickstoff und Kalium einen p-Wert unterhalb von 0,05 aufzuweisen. Somit sind die Daten tendenziell normalverteilt (wenn auch nicht perfekt, da wenig Daten pro Gruppe).
Falls mehrfaktorielle Daten vorliegen, die überhaupt nicht normalverteilt sind, kann man überlegen, diese zu transfomieren (log, Box-Cox, …).
Um genauer zu sein, sollten die Residuen des zugehörigen linearen Modells einer Normalverteilung folgen. Das ließe sich wie folgt bestimmen:
lm(yield ~ N * P * K, data = npk) %>%
residuals() %>%
shapiro.test()
Varianzhomogenität: Hierzu könnten wir die Größen der Boxplots betrachten, welche einigermaßen gleich groß sein sollten. Einfacher ist es, den Levene-Test zu verwenden:
%>%
npk levene_test(yield ~ N * P * K)
Bei dem Levene-Test wird immer das *
zwischen den
Faktoren verwendet, da auch die gemeinsamen Varianzen eine Rolle
spielen
Da der p-Wert (~0,93) deutlich über dem üblichen Signifikanzniveau liegt, können wir hier auch von Varianzhomogenität ausgehen.
Da die Daten grob normalverteilt sind, es keine Ausreißer gibt und die Varianzhomogenität ebenfalls gegeben ist, können wir die mehrfaktorielle ANOVA anwenden:
%>%
npk anova_test(yield ~ N * P * K)
Achten Sie darauf, dass zwischen jedem Faktor das *
steht, sodass Interaktionen der Faktoren betrachtet werden.
Die p-Werte zeigen auf, dass nur der Haupteffekt von dem Stickstoffeinsatz einen signifikanten Einfluss aufweist (Unterschied zwischen Einsatz und keinem Stickstoffeinsatz auf den Ertrag).
Man könnte nun einen Post-Hoc-Test durchführen, jedoch wird dieser eine Menge an nicht-signifikanten Paaren anzeigen:
%>%
npk tukey_hsd(yield ~ N * P * K) %>%
filter(p.adj < 0.05)
Wir sehen erneut, dass sich die Gruppe mit dem Einsatz von
Stickstoff und die Gruppe ohne dessen Einsatz (als Ganze) signifikant (p
< 0,05) unterscheiden. Dem ergänzend scheint sich der Vergleich von
dem Einsatz von Stickstoff und ohne Stickstoff bezüglich des Einsatzes
von Kalium zu unterscheiden (N:K
).
Ersteres scheint hierbei das hauptsächliche Ergebnis zu sein.
Übersichtlicher wäre unter Umständen hierbei die Fixierung der Kalium-und Phosphatgruppen:
%>%
npk group_by(P , K) %>%
anova_test(yield ~ N)
%>%
npk group_by(P, K) %>%
tukey_hsd(yield ~ N)
An dieser Stelle könnten wir bereits fertig sein, denn keine Interaktion ist signifikant und es gibt nur zwei Gruppen für jeden Faktor (Anwendung oder keine Anwendung)
Betrachten wir noch eine kleine Verfeinerung des Modells. Da der
Block, in denen die Erbsen angebaut wurden ebenfalls einen Einfluss auf
das Wachstum der Pflanze haben kann (z. B. durch leicht unterschiedliche
Mineralstoffe, Wässerung, usw.) können wir die Variable
block
als “Störfaktor” in das Modell aufnehmen. Hierzu wird
dieser einfach durch das +
vor die anderen
Faktoren gesetzt:
<- npk %>%
aov_ergebnis anova_test(yield ~ block + N * P * K)
aov_ergebnis
Interessant ist hierbei, dass der Block einen signifikanten Einfluss hat (das interessiert uns im Gesamtbild des Versuchs eigentlich nicht) sowie zusätzlich Kalium.
Allerdings scheinen die gesamte Interaktion N:K:P
sowie die Interaktionen N:P
, N:K
und
P:K
nicht signifikant zu sein.
Da keine Interaktionen signifikant sind, müssen wir uns keine paarweisen Vergleiche anschauen, da alle Faktoren nur zwei Ausprägungen haben:
# ANOVA und Tukey-Test erfüllen denselben Zweck (da nur zwei Ausprägungen vorhanden sind)
%>%
npk group_by(P, K) %>%
anova_test(yield ~ N)
%>%
npk group_by(P, K) %>%
tukey_hsd(yield ~ N)
# Für die Grafik benötigen wir dennoch die Tabelle der paarweisen Vergleichen:
<- npk %>%
tukey_ergebnis tukey_hsd(yield ~ block + N * P * K) %>%
filter(term == "N")
ggplot(npk) +
geom_boxplot(aes(N, yield, fill = N)) +
::stat_pvalue_manual(tukey_ergebnis, label = "p.adj.signif", y.position = 72) +
ggpubrtheme(legend.position = NULL) +
# Zusätzliche Information des Tests:
labs(title = "Stickstoffbehandlung weist einen positiven Einfluss auf das Erbsenwachstum auf",
subtitle = get_test_label(aov_ergebnis, detailed = TRUE, row = 2),
y = "Ertrag")
Der Einsatz von Stickstoff hat einen signifikanten (positiven) Einfluss auf das Wachstum von Erbsen
Wenn wir zusätzlich den Block in dem Modell berücksichtigen (Störfaktor), dann hat ebenfalls Kalium einen signifikanten Einfluss auf das Wachstum von Erbsen.
<- npk %>%
tukey_ergebnis tukey_hsd(yield ~ block + N * P * K) %>%
filter(term == "K")
ggplot(npk) +
geom_boxplot(aes(K, yield, fill = K)) +
::stat_pvalue_manual(tukey_ergebnis, label = "p.adj.signif", y.position = 72) +
ggpubrtheme(legend.position = NULL) +
# Zusätzliche Information des Tests:
labs(title = "Kalium weist einen negativen Einfluss auf das Erbsenwachstum auf",
subtitle = get_test_label(aov_ergebnis, detailed = TRUE, row = 4),
y = "Ertrag")
Bei Kalium scheint dieser Trend jedoch negativ zu sein.
Theoretisch könnten wir uns genauer anschauen, welche Gruppen diese Unterschiede aufweisen, aber das ist nicht ganz so wichtig:
%>%
npk tukey_hsd(yield ~ block + N * P * K)
Insgesamt sehen wir signifikante Haupteffekte von Stickstoff (positiv) und Kalium (negativ), d. h. unabhängig der anderen Faktoren sorgen Stickstoff und Kalium durchschnittlich für einen signifikanten Einfluss auf das Erbsenwachstum.
Es gab (abhängig von der Betrachtungweise) keine signifikante Interaktion, d. h. eine Kombination der Anwendung der drei Faktoren scheint keinen signifikanten Einfluss auf das Erbsenwachstum zu haben.
Wir sehen zusätzlich, dass die unterschiedlichen Blöcke einen hohen Einfluss auf die anderen Faktoren haben und das Experiment “stören”.
Abschließend kann man überlegen, warum Kalium einen negativen Einfluss hat und die Interaktionen keine Signfikanz aufzeigen, da üblicherweise Kaliummangel schlecht für Pflanzen ist und das Verhältnis von NPK entscheidend ist für Düngemittel.
Es mag sein, dass der verwendete Dünger zu viel Kalium enthalten hat, der Boden bereits mit Kalium gesättigt war oder Erbsen keinen Zusatz von Kalium benötigen, sodass der zusätzliche Einsatz von Kalium einen negativen Effekt hatte.
Stickstoff hatte wie erwartet einen positiven Effekt, der unabhängig von Kalium und Phosphat ist.
Da Phosphat keinen Effekt aufzeigt, auch nicht in Kombination mit Stickstoff und Kalium, ist es möglich, dass Phosphat nicht im passenden Maße eingesetzt wurde oder tatsächlich keine Wirkung auf Erbsenwachstum hat.
Bei wiederholtem Experiment könnte auch überlegt werden, die Bedingungen der veschiedenen Blöcke gleichmäßiger zu halten (Bewässerung, Sonne, Erde, usw.).
Weitere Interpretationen obliegen dem Expertenwissen im Bereich der Pflanzenzucht rund um die Erbse. Eine Optimierung der Nährstoffanwendung