Übungszettel 9

1 Aufgaben

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()

npk <- as_tibble(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).

1.1 Übungsaufgaben

1. Werten Sie den Datzensatz npk mit einer passenden ANOVA aus.

   1. Verschaffen Sie sich eine (visuelle/statistische) Übersicht über die Daten.

   2. Prüfen Sie die Voraussetzungen der einfachen ANOVA.

   3. Wenden Sie das passende ANOVA-Modell an.

   4. Interpretieren Sie die Ergebnisse.

1.2 Lösungen

1.2.1 Werten Sie den Datzensatz npk mit einer passenden ANOVA aus.

1.2.1.1 Verschaffen Sie sich eine (visuelle/statistische) Übersicht über die Daten.

# 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()

• Anhand der unterschiedlichen Gruppierung durch Aufteilung (x-Achse) sowie der Färbung (Füll-und Randfarbe).
• Beispielsweise führt das Auslassen der Anwendung von Stickstoff sowie der Einsatz von Phosphat und Kalium zur kleinste Box (geringster Ertrag).
• Der größte Ertrag scheint nur bei dem Einsatz von Stickstoff erzielt worden zu sein.
• Die Boxen scheinen etwas zu variieren. Die Voraussetzungen testen wir in dem nächsten Schritt, um diesem auf den Grund zu gehen.

1.2.1.2 Prüfen Sie die Voraussetzungen der einfachen ANOVA.

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()

• Auch hier sehen wir (p > 0,05), dass eine Normalverteilung vorliegt

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.

1.2.1.3 Wenden Sie das passende ANOVA-Modell an.

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:

aov_ergebnis <- npk %>%
  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)

1.2.1.4 Interpretieren Sie die Ergebnisse.

# Für die Grafik benötigen wir dennoch die Tabelle der paarweisen Vergleichen:
tukey_ergebnis <- npk %>%
  tukey_hsd(yield ~ block + N * P * K) %>%
    filter(term == "N")
ggplot(npk) +
  geom_boxplot(aes(N, yield, fill = N)) +
  ggpubr::stat_pvalue_manual(tukey_ergebnis, label = "p.adj.signif", y.position = 72) +
  theme(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.

tukey_ergebnis <- npk %>%
  tukey_hsd(yield ~ block + N * P * K) %>%
    filter(term == "K")
ggplot(npk) +
  geom_boxplot(aes(K, yield, fill = K)) +
  ggpubr::stat_pvalue_manual(tukey_ergebnis, label = "p.adj.signif", y.position = 72) +
  theme(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