2.1 Was ist R?

R ist eine Programmiersprache, die zur Datenanalyse entwickelt wurde und damit eine Alternative zu Programmen wie SPSS, SAS und den Statistik-Funktionen von Excel. Die zwei größten Vorteile von R gegenüber point-and-click Programmen sind:

1. Der Funktionsumfang ist grundsätzlich nicht beschränkt → durch Programmierung unendlich erweiterbar → Über 10.000 Software-Pakete, die verschiedenste Analyseverfahren ermöglichen
2. Die Analysen mit R sind reproduzierbar und transparent! Analyse kann ohne Aufwand wiederholt werden, z.B. nach Fehler in den Rohdaten. Jeder, der den Code sieht (und versteht), kann nachvollziehen, wie die Daten verwendet worden sind.

Weitere Vorteile sind:

• die Software ist open source also kostenfrei verfügbar
• sie läuft auf allen gängigen Betriebssystemen
• es gibt verschiedene Editoren, die die Nutzung erleichtern, der bekannteste ist R-Studio
• es gibt eine sehr aktive Community im Netz, die Hilfe bietet, oft innerhalb von Minuten bis Stunden

R ist auch eine Investition in Ihre (berufliche) Zukunft. Außer Daten auswerten kann man damit

• komplexe und schöne Graphiken erzeugen (siehe Abbildung unten)
• Berichte, Paper und Bücher schreiben (auch Bachelor- und Masterarbeiten)
• Web-Seiten erstellen (wie diese hier)
• Apps entwickeln
• Simulationsmodelle implementieren
• Eigene Methoden entwickeln und veröffentlichen

Beispiele für R-Graphiken

Ganz allgemein gilt R zunehmend als die Standardsprache für statistische Problemstellungen sowohl in der Wirtschaft als auch in der Wissenschaft (Amirtha et al., 2014).

2.2 R und R-Studio installieren

Im Rahmen dieses Seminars erwarten wir nicht zwingend, dass Sie R installieren und damit arbeiten. Tun Sie es am besten trotzdem - es wird sich lohnen!

1. Zuerst sollten Sie die eigentliche Programmiersprache R installieren. Auf dieser Seite finden Sie Downloads für verschiedene Betriebssysteme:

https://www.r-project.org

Sie gelangen hier auf eine Seite mit einem Link zum R-download. Folgen Sie diesem Link und suchen Sie dann einen CRAN-mirror aus. Im Prinzip ist egal, welchen Sie wählen (auf allen sind die gleichen Ressources verfügbar, deshalb ‘Spiegel’), es macht aber Sinn, einen Mirror in der Nähe zu wählen.

2. Dann (und erst dann! sonst kann es zu Konfigurationsproblemen kommen) installieren Sie auch noch einen Editor für R. Der Editor ist hauptsächlich dafür da, den Code, zum Beispiel für statistische Analysen, zu schreiben, zu kommentieren und zu speichern. R-Studio ist der bekannteste und wahrscheinlich beste Editor:

https://www.rstudio.com/products/rstudio/download/#download

Wenn Sie R-Studio öffnen, sehen Sie vier Fenster: den eigentlich Editor, in dem Sie R-code schreiben, die R Konsole, in der die aktive R-Session läuft, eine Übersicht über den Arbeitsspeicher und unten rechts ein Fenster in dem je nach Aktion Plots, Hilfedateien oder die Ordnerstruktur und Dateien des Verzeichnisses, in dem der R-Prozess gestartet wurde, zu sehen sind. Eventuell fehlt auch der Editor oben links noch, das ändert sich aber, wenn Sie ein neues Skript anlegen (siehe ‘Erste Schritte mit R’)

Ansicht von R-Studio mit 4 Fenstern

2.3 Erste Schritte mit R

• Ein neues Skript im Editor öffnen, um Code schreiben zu können (File → new file → R-Script)
• R als Taschenrechner benutzen: In den Editor

5 + 7

schreiben → auf den ‘Run’ Button klicken (über dem Editor) → der Code wird an den laufenden R Prozess in der Konsole geschickt und das Ergebnis in der Konsole ausgegeben

• Eine Variable anlegen: im Editor

Ergebnis1 <- 5 + 7

eingeben → Run-Button anklicken → nichts passiert außer dass die Zeile in der Konsole auftaucht? → das Ergebnis wird in der Variable Ergebnis1 gespeichert, diese Variable wird jetzt auch in der Übersicht des Arbeitsspeichers gelistet. Wenn man in der Konsole

Ergebnis1

eingibt, wird der Wert der Variable ausgegeben. Anstelle des Pfeils kann auch ein = verwendet werden.

• Kommentare schreiben: Damit andere oder unsere zukünftigen Ichs später noch wissen, was wir warum gemacht haben, kann und sollte man in das R-Skript Kommentare schreiben. Wenn solche Zeilen zum R-Prozess geschickt werden, werden sie einfach ignoriert. Ein Kommentar wird mit einem Raute-Zeichen (hash) eingeleitet.

# Hier benutzen wir R als Taschenrechner
Ergebnis1 <- 5 + 7

• das Skript speichern: File → Save Dabei werden nicht die Variablen, Ergebnisse und Daten gespeichert (das lernen wir später). Sie können aber mit einem Knopfdruck das gesamte Skript laufen lassen und so alle Arbeitsschritte wiederholen → alle Ergebnisse wieder herstellen.

2.4 Importieren von Daten aus Excel/LibreOffice nach R

1. Daten in Excel/LibreOffice eingeben
   • Häufige Fehlerquelle! Am besten zu zweit: einer liest vor, einer tippt. Zwei/drei durch zwo/drei disambiguieren. Wenn nur eine Person die Daten eingibt, am besten den Nummernblock der Tastatur verwenden und Zahlen blind eingeben, also nur auf Zettel und Bildschirm schauen. Eingegebene Daten danach nochmal prüfen.
   • Keine leeren Zellen im Datenblock → wenn die Werte fehlen: “NA” (= not available)
   • Erste Reihe wird per default als Spaltenüberschrift interpretiert, hier aussagekräftige aber nicht zu lange Überschriften wählen
   • Zeichenketten (wie Überschriften) sollten aus einem einzigen Wort bestehen. Wenn nötig, zum Beispiel Unterstriche verwenden um einzelne Wörter zu verbinden
   • Als Dezimaltrenner sollte ein Punkt genutzt werden → “7.5” nicht “7,5”
     
   • Pro Beobachtung sollte es eine Zeile geben
2. Die Daten in csv-Format speichern
3. Die Datei mit der Funktion read.csv() importieren.

kirschen <- read.csv(“Kirschen.csv”)

Die Datei ‘Kirschen.csv’ können Sie hier herunterladen:
https://ecampus.uni-bonn.de/goto_ecampus_file_2786374_download.html

Wenn die R-Session nicht im gleichen Ordner läuft, in dem auch die Datei gespeichert ist, wird R eine Fehlermeldung ausgeben, dass es die Datei nicht finden kann. Entweder, Sie geben dann den gesamten Pfad bis zur Datei ein (z.B. “~/Downloads/Kirschen.csv”) oder - am einfachsten - folgenden Code:

kirschen <- read.csv(file.choose())

Es öffnet sich dann ein Browser-Fenster (manchmal nicht im Vordergrund), in dem Sie die Datei suchen und auswählen können. Das ist nicht ganz ideal, weil nicht unbedingt reproduzierbar (eventuell gibt es zwei oder mehr verschiedene Dateien mit dem Namen “Kirschen.csv” auf Ihrem Computer…), ist aber immerhin leicht zu händeln. Sauberere Lösungen zeigen wir Ihnen im nächsten Kapitel.

Beispiel für richtig eingegebene Daten

2.5 Daten kontrollieren

Wenn Sie die Daten in R eingelesen haben, sollten Sie als erstes kontrollieren, ob der Import geklappt hat. Eine erste Übersicht bekommen Sie mit den Funktionen summary() und head()

summary(kirschen)
# gibt eine Zusammenfassung der Daten mit Mittelwerten etc. aus

head(kirschen)
# Zeigt die ersten 6 Zeilen des Datensatzes, wie er in R eingelesen wurde

Hier prüft man unter anderem, ob die Daten den richtigen Typ haben, zum Beispiel die Blockbezeichnungen als Faktoren und nicht als Zahlen interpretiert werden. Ist das nicht der Fall, kann es jetzt geändert werden (→ späteres Kapitel).

Auch einfache Plots sind ein gutes Mittel, um eine erste Vorstellung der Daten zu bekommen und Auffälligkeiten wie extreme Werte (verrutscher Punkt in den Zahlen) zu entdecken. Dafür eignen sich

• die einfach plot-Funktion

  plot(kirschen$Inokulat, kirschen$Frischgewicht, 
   main = "Inokulat", ylab = "Frischgewicht der Blätter [g]")

• oder eine schönere Version mit dem Paket ‘ggplot2’. Dazu muss das Paket ‘ggplot2’ erst installiert werden. In R-Studio geht das über ‘Tools -> Install Packages’ -> a search in the ‘Packages’ field.

  library(ggplot2)    
  ggplot(kirschen, aes(x=Inokulat, y=Frischgewicht)) +  
    geom_boxplot(notch=TRUE)

2.6 Hilfe bekommen

Hilfe zu finden ist zu Beginn des Arbeitens mit R die wichtigste Kompetenz überhaupt. Am einfachsten ist es natürlich, wenn Sie direkte Unterstützung durch jemanden bekommen, der sich mit R auskennt (und das ist bei der Bachelor-Arbeit - zumindest im Gartenbau - natürlich der Fall). Es gibt aber auch eine Reihe anderer Möglichkeiten:

1. in R selbst
   • ein “?” vor den Funktionsnamen setzen ?sd() → es öffnet sich eine Hilfedatei,
   • help.start() → öffnet Hilfeseite mit wichtigen Dokumenten
2. im Internet
   • Bei der Suche in der Suchmaschine ein R zu der Frage oder der Fehlermeldung hinzufügen (Fehlermeldungen gründlich lesen!)
   • Nach Tutorials und online Büchern suchen, zum Beispiel YaRrr - the pirates guide to R https://bookdown.org/ndphillips/YaRrr/
   • Cheat cheets
     …bieten einen Überblick der verfügbaren Funktionen und kurze Beschreibungen zu einem bestimmten Thema. Hier finden Sie eine Sammlung von cheat cheets: https://www.rstudio.com/resources/cheatsheets/ Hier ein cheat sheet mit den grundlegenden Funktionen: https://iqss.github.io/dss-workshops/R/Rintro/base-r-cheat-sheet.pdf

Und hier der Link zum ersten Quiz: https://ecampus.uni-bonn.de/goto_ecampus_tst_2786376.html

3.1 R-Projekte

R-Projekte sind gewissermaßen das zu Hause von allen R-Skripten und Daten, die zu einem Projekt (zum Beispiel einem Experiment) gehören. Dabei ist ein Projekt aber mehr als nur ein Ordner. Die Vorteile eines Projekts im Vergleich zu einer einzelnen neuen R Script-Datei sind, dass

• beim Öffnen des Projekts das Arbeitsverzeichnis automatisch auf den Ordner gesetzt wird, in dem sich die Projekt-Datei befindet. Man kann dann Dateien direkt aufrufen/einlesen ohne den vollständigen Pfad angeben zu müssen (höchstens Unterordner, wenn man denn welche angelegt hat). In Bezug auf das letzte Kapitel sollte jetzt also kirschen <- read.csv("Kirschen.csv") immer funktionieren, wenn die Datei ‘Kirschen.csv’ im gleicher Ordner, wie die .Rproj-Datei liegt.
• RStudio bei einem Neustart alle Dateien wieder öffnet, die beim letzten Schließen offen waren
• man gleichzeitig mehrere Projekte offen haben kann, ohne durcheinander zu kommen
• man das ganze Projekt komplett mit richtigen Verknüpfungen zwischen Daten/Dateien weiterreichen kann und die Skripte (zumindest bei gleicher R-Version) auf jedem Rechner reproduziert werden können.

Sie können ein RStudio-Projekt sowohl in einem neuen Verzeichnis als auch in einem vorhandenen Verzeichnis, in dem Sie bereits R-Code und Daten haben, erstellen. Um ein neues Projekt zu erstellen, verwenden Sie den Befehl ‘File -> New Project’.

Es öffnet sich dann ein Fenster, in dem Sie entscheiden können, ob Sie das Projekt in einem neuen Verzeichnis, einem schon existierenden Verzeichnis oder mit Versionskontrolle (zum Beispiel mit git) anlegen wollen. Letzteres ist sehr sinnvoll, würde hier aber den Rahmen sprengen.

Durch das Anlegen eines neuen Projekts in RStudio wird eine Projektdatei (mit der Erweiterung .Rproj) im Projektverzeichnis erstellt. Diese Datei kann später als Verknüpfung zum direkten Öffnen des Projekts verwendet werden.

3.2 Datentypen und -transformationen

Das Einlesen der Daten sollte jetzt kein Problem mehr sein, allerdings haben die Daten manchmal nicht den Datentyp der für bestimmte statistische Auswertungen oder Abbildungen gebraucht wird. Zuerst eine Übersicht über die gängigsten Typen:

Mit der Funktion str()kann man testen, als welche Datentypen R die Daten interpretiert und bekommt eine Übersicht der Struktur des Datensatzes:

str(kirschen)

Wie oben erwähnt, werden für einige Auswertungen/Abbildungen manchmal andere Datentypen benötigt, als R erkannt hat. So werden zum Beispiel manchmal die Bezeichnungen für Behandlungsgruppen als ‘character’ eingelesen, für eine Varianzanalyse werden aber Faktoren benötigt. In diesem Fall kann der Datentyp mit der Funktion as.factor() umgewandelt werden.

kirschen$Inokulat <- as.factor(kirschen$Inokulat)

Das geht natürlich nur, wenn der neue Datentyp zu der umzuwandelnden Variablen passt: ‘Gi_12’ kann nicht in eine numerische Variable umgewandelt werden, wohl aber die Zahlen 1, 2, 3, .. in Faktoren (die dann als “1”, “2”, “3” ausgegeben werden). Andere häufig genutzte Funktionen sind as.numeric() , as.character() und as.Date().

3.3 Umgang mit NAs

Auch NAs können manchmal zu Fehlermeldungen oder unvollständigen Abbildungen führen. Da NA nicht Null bedeutet, sondern eher ‘wissen wir nicht’, kann zum Beispiel nicht ohne Weiteres ein Mittelwert berechnet werden.

Daten <- c(5, 3, 9, NA, 4, NA)
mean(Daten)

Die meisten Funktionen haben aber die Option, NA-Werte bei der Berechnung auszuschließen:

mean(Daten, na.rm = TRUE)

Schauen Sie also am besten zuerst in der Hilfe zu der Funktion , die Sie verwenden möchten (?mean), ob es ein Funktionsargument gibt, das den Umgang mit NA-Werten festlegt und ändern Sie die Einstellung entsprechend. Wenn es gar nicht anders geht, kann man als Alternative auch die Zeilen in den NAs vorkommen vor der Analyse herausfiltern. Wie solche Daten-Operationen funktionieren, lernen Sie im folgenden Kapitel.

3.4 Datensätze umstrukturieren

Manchmal ist es nötig, den eingelesenen Datensatz umzustrukturieren, Teildatensätze herauszufiltern oder nue Variablen zu berechnen. Dazu sollte man am besten das Paketbündel ‘tidyverse’ installieren. Hat man das getan, steht auch eine weitere Konvertierungsfunktion zur Verfügung, die aus einem normalen Datensatz (data.frame) ein ‘tibble’ macht. Tibble haben ein paar Vorteile, zum Beispiel, dass standardmäßig nur die ersten 10 Zeilen ausgegeben werden und zu jeder Spalte der Datentyp gleich dabei steht (wobei Fließkommazahlen hier nicht als numeric sondern double (dbl) bezeichnet werden. ‘fct’ steht für Faktor, ‘int’ für Integer).

kirschen <- as_tibble(kirschen)
kirschen

Daten <- c(5, 9, 2, 5, 3, 9, 5, 3, 5)
Ergebnis <- round(sqrt(sum(Daten^2)), 3)
Ergebnis

Daten_quad <- Daten^2
Daten_quad_sum <- sum(Daten_quad)
Daten_quad_sum_wurzel <- sqrt(Daten_quad_sum)
Ergebnis <- round(Daten_quad_sum_wurzel, 3)

Ergebnis <- Daten^2 %>%   
            sum() %>%     
            sqrt() %>%    
            round(3)      

library(tidyverse)
# mit der Funktion 'tibble()' legt man einen Datensatz an
Daten1 <- tibble(spalte1 = c(1, NA, 5, 7), spalte2 = c("a", "b", NA, "d"))
# wir schauen uns den Datensatz an:
Daten1
# wir löschen die Zeilen mit NA heraus
Daten2 <- Daten1 %>% drop_na()
# uns sehen uns das Ergebnis an:
Daten2

kirschen %>% select(Trockengewicht, Blattfläche) %>% head()

kirschen %>% select(-Trockengewicht, -Blattfläche) %>% head()

kirschen_neu <- kirschen %>% select(-Trockengewicht, -Blattfläche)
head(kirschen_neu)

kirschen_nurPDV <- kirschen %>% filter(Inokulat == "PDV")
summary(kirschen_nurPDV)

plot(kirschen_nurPDV$Inokulat, kirschen_nurPDV$Frischgewicht)

kirschen_nurPDV <- kirschen %>% 
                   filter(Inokulat == "PDV") %>%
                   droplevels()
                   
 plot(kirschen_nurPDV$Inokulat, kirschen_nurPDV$Frischgewicht)                  

kirschen %>% mutate(Gesamttrieblänge_m = Gesamttrieblänge/100)

Wachstum <- tibble(Baum = c("Nr1", "Nr2", "Nr3", "Nr4"), 
                   Zeitp1 = c(338,459,558,439), 
                   Zeitp2 = c(437,560,729,658), 
                   Zeitp3 = c(469,579,937,774))
                   
Wachstum_lang <- Wachstum %>% 
                 pivot_longer(
                 cols = c(Zeitp1, Zeitp2, Zeitp3),
                 names_to = "Zeitpunkt",
                 values_to = "Größe_cm")
                 
# cols: die Zeilen die zusammengefügt werden sollen
# names_to: neuer Spaltenname des Faktors
# values_to: neuer Spaltenname der Daten

# 'Zeitpunkt' und 'Baum' wurden noch nicht als Faktor erkannt, deshalb:
Wachstum_lang$Zeitpunkt <- as.factor(Wachstum_lang$Zeitpunkt)
Wachstum_lang$Baum <- as.factor(Wachstum_lang$Baum)

ggplot(Wachstum_lang, aes(x = Zeitpunkt, y = Größe_cm, color=Baum)) + 
   geom_point() +   labs(y="Größe (cm)", x="Zeitpunkt")

Wachstum_ww <- Wachstum_lang %>% 
               pivot_wider(
               names_from = "Zeitpunkt",
               values_from = "Größe_cm"
               )

4.1 Das Grundprinzip der Varianzanalyse

Wir nehmen an, wir haben einen sehr einfachen Versuch mit Tomaten durchgeführt, in dem die Hälfte aller 40 Pflanzen gedüngt wird und die andere Hälfte als Kontrollgruppe ungedüngt bleibt. Am Ende bestimmen wir das Gesamtgewicht der Tomaten pro Pflanze als ein Maß für den Ertrag. Jetzt interessiert uns, ob die Düngung einen signifikanten Effekt auf den Ertrag hat. Unsere Hypothese ist: Die Düngung hat einen Effekt auf den Ertrag. Statistisch ausgedrückt bedeutet das, dass die Ertragsdaten aus zwei Normalverteilungen mit unterschiedlichen Mittelwerten stammen. Diese Hypothese ist auf folgender Abbildung auf der linken Seite dargestellt. Das zugehörige statistische Modell hat zwei Parameter: a ist der Mittelwert der Normalverteilung ohne Dünger, b die Differenz zwischen a und dem Mittelwert für die Normalverteilung mit Dünger (die Effektgröße). Die Variable Dünger nimmt entweder den Wert 0 (ohne Dünger) oder 1 (mit Dünger) an. Die Streuung der Daten (sie liegen nicht alle genau auf dem Mittelwert) wird durch den Fehlerterm aufgefangen.

Hypothesen, statistische Modelle und Umsetzung in R. Links: Düngung hat einen Effekt auf den Ertrag -> die Daten stammen aus zwei Normalverteilungen mit unterschiedlichen Mittelwerten; rechts (Nullhypothese): Düngung hat keinen Effekt auf den Ertrag -> die Daten stammen aus einer einzigen Normalverteilung.

Die Null-Hypothese zu diesem Modell ist: Die Düngung hat keinen Effekt auf den Ertrag -> die Daten stammen alle aus einer einzigen Normalverteilung mit dem Mittelwert a (rechte Seite der Abbildung). Entsprechend einfacher ist das statistische Modell.

In R formulieren wir diese beiden konkurrierenden Modelle mit der Funktion lm(). lm steht für ‘linear model’ und es schätzt Mittelwert (oder in anderen Fällen auch Achsenabschnitt und Steigung einer Regressionsgeraden), indem es die Werte so wählt, dass die Fehlerquadrate (quadrierter Abstand der Messungen zu den Mittelwerten) möglichst klein sind.

# hier lesen wir keine Daten aus Excel ein sondern geben sie einfach händisch direkt 
# in R ein. Zuerst Ertragsdaten von 10 ungedüngten Pflanzen
ohne <- c(417.1,558.8,519.1,611.7,450.3,461.9,517.3,453.2,533.1,514.7)
# dann von 10 gedüngten Pflanzen
mit <- c(581.9,517.4,641.3,659.8,587.8,483.3,703.0,589.5,532.4,569.3)
# hier werden die beiden Datensätze mit c() wie combine
#  zu einem Vektor (= Spalte) zusammengeführt 
Ertrag <- c(ohne, mit)
# jetzt generieren wir noch einen Vektor der die Behandlung anzeigt
Düngung <- gl(2, 10, 20, labels = c("ohne","mit"))
# Mit diesen beiden Vektoren können wir das statistische Modell formulieren, 
# dass der Ertrag von der Düngung abhängt. In R wird das mit der Tilde ~ ausgedrückt
#  (oben links auf der Tastatur). 
# Wir nennen dieses Modell 'model1'
model1 <- lm(Ertrag ~ Düngung)
# mit der folgenden Zeile schauen wir uns eine summary (Zusammenfassung)
#  des gefitteten Modells an

summary(model1)

Interpretation der Ausgabe einer model summary

In der Zeile ‘Intercept’ finden wir unter ‘Estimate’ den geschätzten Mittelwert für die Gruppe ohne Düngung, also das a in der Abbildung oben, linke Seite. Darunter, in der Zeile ‘Düngungmit’ ist die Differenz (das b aus der Abbildung) das auf das a aufaddiert werden muss, um die Schätzung für den Mittelwert der gedüngten Pflanzen zu erhalten.

Als nächstes formulieren wir die Nullhypothese

null_model <- lm(Ertrag ~ 1)
summary(null_model)

## 
## Call:
## lm(formula = Ertrag ~ 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -128.04  -38.30  -12.39   43.08  157.85 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   545.15      16.68   32.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 74.59 on 19 degrees of freedom

Wie erwartet, erhalten wir jetzt nur einen Mittelwert über alle Tomatenpflanzen.

Welches der beiden Modelle ist nun besser? Das Modell mit Düngung als erklärender Variable passt sich besser an die Daten an -> kleinere Summe der Fehlerquadrate, ist aber komplizierter.

4.2 Das Prinzip der Parsimonie

Die Grundlagen für das Prinzip der Parsimonie, auch bezeichnet als Occam’s Rasor. Verändert nach einer Folie von Frank Schurr

Um zu entscheiden, ob eine erklärende Variable einen signifikanten Effekt hat, stellen wir die Frage: Wie wahrscheinlich ist es, die Daten so zu beobachten, wie wir es tatsächlich getan haben, wenn in Wahrheit die erklärende Variable keinen Effekt hat (= die Nullhypothese wahr ist bzw. alle Daten aus einer einzigen Normalverteilung stammen). Um diese Wahrscheinlichkeit ausrechnen zu können, müssen die Daten innerhalb der Behandlungsgruppen normalverteilt sein und die Streuung (Varianz) muss in etwa gleich sein. Das sind die wichtigsten Vorraussetzungen für eine Varianzanalyse, weil sonst evtl. die berechnete Wahrscheinlichkeit nicht stimmt. In R nutzten wir die Funktion anova() um die beiden Modelle zu vergleichen und die oben genannte Wahrscheinlichkeit auszurechnen:

anova(model1, null_model)

Interpretation der Ausgabe einer ANOVA

Die Wahrscheinlichkeit nach der wir suchen, finden wir unter Pr(>F) und ist der berühmte p-Wert. Die Wahrscheinlichkeit zufällig solche (oder noch extremere) Daten zu finden, wie wir sie für die Tomaten gefunden haben, obwohl die Düngung in Wahrheit keinen Einfluss auf den Ertrag hat, ist nur 0,86% und damit so unwahrscheinlich, dass wir den Effekt der Düngung als eindeutig signifikant interpretieren können. Damit ist model1 das bessere Model. Eine allgemeine, wenn auch etwas arbiträre Übereinkunft, ist, dass alle Variablen mit einem p-Wert von unter 5% bzw p < 0.05 als signifikant gelten.

Kritik an der Verwendung des p-Wertes
In den letzten Jahren gibt es mehr und mehr Kritik an dem starken Fokus auf den p-Wert. Zum einen, weil zumeist einfach das Signifikanzniveau verwendet wird, was fast immer verwendet wird, ohne weiter darüber nachzudenken: 5%. Dieses Niveau stellt einen Kompromiss zwischen der Fähigkeit eine Entdeckung zu machen und unserer Bereitschaft ein fehlerhaftes Testergebnis zu akzeptieren dar. Es sollte eigentlich einleuchtend sein, dass das Signifikanzniveau für verschiedene Fragestellungen auch unterschiedlich angesetzt werden sollte (was der Begründer der frequentistischen Statistik, Ronald Fisher auch schon 1956 selbst so geraten hat). Zum anderen ist der Nachweis eines signifikanten Effekts außer von der Effektgröße stark vom Informationsgehalt der Daten abhängig (wie groß ist die Stichproble). So kann es sein, dass, wenn man hunderte von Proben hat, eine Variable als signifikant nachgewiesen werden kann, ihr Effekt und die wissenschaftliche Relevanz aber nur minimal ist. Es ist also immer, immer wichtig, den Bezug zur eigentlichen Frage und zum untersuchten System zu behalten und die Ergebnisse mit gesundem Menschen- und Sachverstand zu interpretieren. Im Kapitel zur bayesischen Statistik werden wir eine andere Herangehensweise kennenlernen.

4.3 Annahmen der ANOVA überprüfen

Wie oben schon kurz angesprochen, müssen für glaubwürdige Ergebnisse einer Varianzanalyse einige Vorraussetzungen erfüllt sein

1. Unabhänigkeit der Messungen - das werden wir im Kapitel ‘Experimentelles Design’ genauer besprechen
2. Keine großen Ausreißer in den Daten - das geht am einfachsten durch das Plotten der Daten und evtl. Ausreißer herausnehmen. Bleibt noch:
3. Die abhängige Variable ist innerhalb der Behandlungsgruppen normalverteilt
4. Die Varianz in jeder Gruppe sollte annähernd gleich sein. Das wird auch mit dem Begriff Homoskedastizität beschrieben

Um die Normalverteilung der Daten innerhalb der Behandlungsgruppen zu prüfen, ist es einfacher, sich die Daten nicht gruppenweise anzuschauen, sondern einfach die Residuen aller Daten (also den Abstand jedes Datenpunktes zum geschätzten Gruppen-Mittelwert) zu nehmen. Es gibt numerische Tests, um die Normalverteilung der Residuen zu prüfen, in den letzten Jahren setzt sich aber mehr und mehr eine visuelle Überprüfung durch. Am besten eignet sich dafür ein qq-Plot (Quantile-Quantile-Plot). Wenn die Punkte annähernd auf einer Geraden liegen, sind die Residuen normalverteilt. Da ‘annähernd’ ein dehnbarer Begriff ist, liefert die Funktion qqPlot aus dem Package ‘car’ ein Konfidenzintervall. Liegen die Punkte innerhalb dieses Intervalls, kann man von einer Normalverteilung ausgehen.

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

qqPlot(resid(model1))

## [1] 17  4

Außer dem Plot gibt die Funktion die ‘Namen’ der Punkte aus, die im Bezug auf die x-Achse outlier darstellen (also Punkte, die besonders stark abweidchen,hier Punkt 4 und Punkt 17). Es lohnt sich, diese Punkte genauer anzuschauen (waren die Versuchsbedingungen hier aus irgendeinem Grund anders?) und evtl. auszuschließen, insbesondere, wenn sie außerdem außerhalb der Konfidenzlinien liegen.

Eine gute Erklärung wie qq-plots Funktionieren finden Sie hier: https://www.youtube.com/watch?v=okjYjClSjOg … auf Englisch und mit toller Intro :)

Um die Homoskedastizität visuell zu überprüfen, plotten Sie die Residuen gegen die von der Funktion lm() geschätzten Mittelwerte. Wenn die Daten der Gruppen etwa einen gleichen Bereich auf der y-Achse umfassen, sind die Varianzen homogen. Hierfür gibt es leider derzeit keine Konfidenzintervalle, eine Daumenregel ist aber, dass die Varianzen als homogen gelten, wenn die größte Varianz nicht mehr als 1,5fach größer ist als die kleinste Varianz. In unserem Beispiel reichen die Residuen beider Gruppen von etwa -100 bis 100, die Varianzen sind also annähernd gleich groß. Formal kann Varianzhomogenität mit dem Leven’s Test überprüft werden. R bietet dazu die Funktion leveneTest() aus dem package car. Wer interessiert daran ist, kann sich mit ?leveneTest die Hilfeseite mit Erklärungen zur Anwendung der Funktion anschauen.

plot(fitted(model1), resid(model1))
abline(h=0, lty = 3, col = "red")

Für unsere Analyse sind alle Annahmen erfüllt. Wir können also davon ausgehen, dass die Düngung in unserem hypothetischen Experiment tatsächlich einen signifikanten Effekt auf den Ertrag hat. In einer Arbeit würden wir den F-Wert (aus summary(model1)), die Anzahl der Freiheitsgrade und den p-Wert angeben, wobei p-Werte, die kleiner als 0.05 sind, üblicherweise mit p < 0.05 angegeben werden.

Und hier der Link zu Quiz 3: https://ecampus.uni-bonn.de/goto_ecampus_tst_2816432.html

5.1 Nicht-normalverteilte Daten

Im letzten Kapitel haben wir gesehen, dass eine Voraussetzung für die Durchführung einer ANOVA die Normalverteilung der Residuen (Fehler, Abstände der Messpunkte zum Mittelpunkt) ist. Bei kontinuierlichen Daten wie Größe, Gewicht etc. ist das normalerweise der Fall, wenn das Modell richtig formuliert ist (d.h. wenn alle Faktoren, die einen Einfluss haben, enthalten sind und deren Beziehung - additiv, multitplikativ, exponentiell - korrekt dargestellt ist). Sind die Daten nicht kontinuierlich, können wir erstmal nicht von einer Normalverteilung der Residuen ausgehen. Häufig sind auch die Varianzen (= Streuung der Daten) nicht homogen über die Behandlungsgruppen. Das ist zum Beispiel bei Zähldaten der Fall, zumindest, wenn es sich um eher kleine Zahlen handelt (bis etwa 7).

5.2 Daten transformieren?

In so einem Fall wurde früher oft dazu geraten, die Daten zu transformieren, um sie ANOVA-fähig zu machen. Das birgt allerdings verschiedene Schwierigkeiten:

• die Interpretation der Ergebnisse wird schwieriger, weil die Ergebnisse und Signifikanzen nur für die transformierten Daten zutreffen. Damit verbunden ändert sich auch die Funktion des Zusammenhanges der erklärenden Parameter. Wenn die erklärenden Variablen vorher addiert wurden, um die abhängige Variable zu schätzen, müssen sie zum Beispiel nach einer Logarithmierung multipliziert werden.
• das Logarithmieren von Zähldaten, das häufig propagiert wird (wobei häufig eine 1 addiert wird, um Nullen zu vermeiden), führt häufig zu einer nicht gewollten Abnahme der Varianz der abhängigen Variable, die abhängig vom Mittelwert ist.
• in jedem Fall muss die Annahme der normalverteilten Residuen auch nach einer Transformation überprüft werden.

Die meisten Ratschläge, Daten zu transformieren, stammen aus der Zeit, als Varianzanalysen noch von Hand durchgeführt wurden und kompliziertere Modelle kaum zu rechnen waren. Zum Glück haben wir heute Computer und R, um auch viele nicht-normalverteilte Daten mittels einer ANOVA (Varianzanalyse) analysieren zu können. Der Trick hierbei ist, nicht die schon bekannten linearen Modelle (lm) zu verwenden, sondern generalisierte lineare Modelle (glm).

5.3 Generalisierte lineare Modelle

… können Poisson-verteilte Residuen (aus Zähldaten), binomial-verteilte Residuen (aus Anteilsdaten) und viele andere Residuenverteilungen beschreiben. Ich werde hier nicht auf die Charakteristika dieser Verteilungen eingehen, weil es für das Verständnis der Auswertung nicht essentiell ist. Wenn Sie Interesse haben, finden Sie aber gute Erklärungen über eine Internet-Suche.

Vergleich der Funktionsweisen eines linearen und eines generalisierten linearen Modells. Verändert nach einer Folie von Frank Schurr.

In einem normalen linearen Modell (oberer Teil der Abbildung) liefern die erklärenden Variablen (zum Beispiel ‘Düngung’)und die Modellgleichung direkt eine Schätzung der abhängigen Variable \(\mu\) (zum Beispiel ‘Ertrag). Die Abstände zu den tatsächlich beobachteten Daten \(y\) sind die Fehler oder Residuen und es wird angenommen, dass sie normalverteilt sind. Bei einem GLM wird zusätzlich eine ’link-function’ benötigt, die das Ergebnis der Modellgleichung auf den biologisch sinnvollen Bereich abbildet. So wird zum Beispiel verhindert, dass negative Zahlen geschätzt werden (ein Baum mit -5 Äpfeln). Zusätzlich ist die Annahme der Verteilung der Residuen/Fehler flexiber, wie oben bereits erwähnt.

Exkurs: Schätzung der Modell-Parameter
Grundsätzlich findet man die Modell-Parameter (bei einer Geradengleichung \(y = a * x + b\), wären das zum Beispiel \(a\) und \(b\)), für die die Daten am wahrscheinlichsten sind, indem man das Maximum-Likelihood Prinzip anwendet, sprich die deviance - also die Abweichung der geschätzten von den beobachteten Werten - minimiert. Bei normalverteilten Fehlern ist die deviance als die Summe der Fehlerquadrate RSS (residual sum of squares) definiert: \(\sum (\mu - y)^2\). Bei nicht-normalverteilten Fehlern funktioniert der Trick mit dem Minimieren der Summe der Fehlerquadrate nicht. Die deviance ist hier etwas komplizierter definiert. Für Poisson-verteilte Daten zum Beispiel \(2\sum(y_i log(y/\mu)-(y-\mu))\).

5.4 Beispiel: Tomaten

Leichter verständlich, wird es sicher, wenn wir uns ein konkretes Beispiel anschauen - hier können Sie einen Datensatz zu einem Tomatenversuch herunterladen: https://ecampus.uni-bonn.de/goto_ecampus_file_2786370_download.html Speichern Sie ihn in ihrem R-Projekt (wenn Sie eins angelegt haben) im Ordner ‘data’.

# Zuerst lesen wir den Tomaten-Datensatz ein
tomaten_daten <- read.csv("data/Tomaten.csv")

# Eine kurze Kontrolle, ob alles funktioniert hat
summary(tomaten_daten)

##     Datum             Zeitpunkt           Haus        Bedachung        
##  Length:865         Min.   : 1.000   Min.   :1.000   Length:865        
##  Class :character   1st Qu.: 3.000   1st Qu.:2.000   Class :character  
##  Mode  :character   Median : 5.000   Median :4.000   Mode  :character  
##                     Mean   : 5.089   Mean   :3.504                     
##                     3rd Qu.: 7.000   3rd Qu.:5.000                     
##                     Max.   :10.000   Max.   :6.000                     
##                                                                        
##   Pflanzen.Nr     Veredelung         Salzstress          Trieblänge   
##  Min.   : 1.00   Length:865         Length:865         Min.   :  0.0  
##  1st Qu.:25.00   Class :character   Class :character   1st Qu.: 80.0  
##  Median :49.00   Mode  :character   Mode  :character   Median :104.0  
##  Mean   :48.55                                         Mean   :106.4  
##  3rd Qu.:72.00                                         3rd Qu.:130.0  
##  Max.   :96.00                                         Max.   :196.0  
##                                                                       
##  Anzahl_Blätter  Anzahl_Blütenstände Anzahl_grüne_Früchte Anzahl_rote_Früchte
##  Min.   : 0.00   Min.   : 0.000      Min.   : 0.00        Min.   : 0.000     
##  1st Qu.:15.00   1st Qu.: 3.000      1st Qu.: 2.00        1st Qu.: 0.000     
##  Median :19.00   Median : 4.000      Median : 7.00        Median : 0.000     
##  Mean   :19.62   Mean   : 4.467      Mean   :11.14        Mean   : 1.029     
##  3rd Qu.:24.00   3rd Qu.: 6.000      3rd Qu.:17.00        3rd Qu.: 0.000     
##  Max.   :36.00   Max.   :36.000      Max.   :58.00        Max.   :15.000     
##                                                                              
##  Anzahl_Blütenendfäule Anzahl_geerntete_Früchte Gewicht_Trieb   Gewicht_Blätter
##  Min.   : 0.000        Min.   : 0.000           Min.   :  0.0   Min.   :   0   
##  1st Qu.: 0.000        1st Qu.: 0.000           1st Qu.:235.0   1st Qu.: 820   
##  Median : 1.000        Median : 0.000           Median :368.0   Median :1237   
##  Mean   : 3.042        Mean   : 0.341           Mean   :328.6   Mean   :1077   
##  3rd Qu.: 5.000        3rd Qu.: 0.000           3rd Qu.:428.0   3rd Qu.:1418   
##  Max.   :17.000        Max.   :10.000           Max.   :579.0   Max.   :1783   
##                                                 NA's   :768     NA's   :768    
##  Gewicht_grüne_Früchte Gewicht_rote_Früchte Gewicht_Blütenendfäule
##  Min.   :   0.0        Min.   :  0.0        Min.   :  0.00        
##  1st Qu.: 365.0        1st Qu.:  0.0        1st Qu.:  0.00        
##  Median : 688.0        Median :  0.0        Median :  0.00        
##  Mean   : 633.2        Mean   : 18.9        Mean   : 21.93        
##  3rd Qu.: 964.0        3rd Qu.:  0.0        3rd Qu.: 28.00        
##  Max.   :1222.0        Max.   :312.0        Max.   :221.00        
##  NA's   :768           NA's   :768          NA's   :768

# Es fällt auf, dass die beiden erklärenden Variablen 'Veredelung' und 'Salzstress', 
# die wir gleich verwenden wollen, as 'character' eingelesen worden sind. 
# Das kann bei den Abbildungen zu Problemen führen. Deshalb wandeln wir sie 
# zuerst in 'factors' um:
tomaten_daten$Veredelung <- as.factor(tomaten_daten$Veredelung)
tomaten_daten$Salzstress <-
  as.factor(tomaten_daten$Salzstress)

# Für die folgende Abbildung benötigen wir das package 'lattice' 
# (muss wie immer zuerst installiert und dann geladen werden).
library(lattice)

# Um die Daten ein bisschen kennenzulernen (und noch ein paar neue Arten von plots kennenzulernen), 
# bilden wir sie erstmal in einem xyplot ab. Insbesondere interessieren wir uns 
# jetzt für die Anzahl an Früchten, die Blütenendfäule haben in Abhängigkeit von 
# Salzstress und Veredelung.
# 
with(tomaten_daten, xyplot(Anzahl_Blütenendfäule ~ Salzstress, groups = Veredelung))

Hier sehen wir nicht besonders viel, außer dass die Varianz sehr hoch ist. Machen wir noch einen Interaktionsplot, in dem die Mittelwerte der Gruppen (mit/ohne Salzstress, mit/ohne Veredelung) berechnet und abgebildet werden.

# hier machen wir ein Interaktionsplot
with(tomaten_daten, interaction.plot(x.factor = Veredelung, 
   trace.factor = Salzstress, response = Anzahl_Blütenendfäule))

Jetzt erkennen wir, dass veredelte Tomaten mit Salzstress eine deutlich geringere Zahl von Früchten mit Blütenendfäule haben, als ohne Salzstress. Bei nicht-veredelten Pflanzen ist der Effekt genau umgekehrt, allerdings deutlich kleiner. Es sieht also so aus als gäbe es eine Interaktion zwischen den beiden erklärenden Variablen Salzstress und Veredelung: der Effekt der einen erklärenden Variablen auf die abhängige Variable hängt vom Wert der anderen erklärenden Variablen ab. Um diesen Effekte auf Signifikanz zu testen, fitten wir jetzt generalisierte lineare Modelle, die mit der Funktion glm() aufgerufen werden:

glm_model_1 <- glm(Anzahl_Blütenendfäule ~ Veredelung * Salzstress, 
    family = poisson, data = tomaten_daten)

Wie Sie sehen, funktioniert das ganz ähnlich wie mit der Funktion lm(). Es muss lediglich der zusätzliche Parameter ‘family’ angegeben werden. Hier wählen wir ‘poisson’, weil Zähldaten Poisson-verteilt sind.

5.5 Overdispersion

Aber Vorsicht! Bevor wir das Modell weiter analysieren können müssen wir uns noch das Verhältnis von Residual deviance (die Fehler, die es nach der Schätzung der Mittelwerte noch gibt) und den Freiheitsgraden/degrees of freedom anschauen: Ein GLM nimmt eine bestimmte Beziehung zwischen der Residual deviance und der Modellvorhersage an. Bei vielen agrarwissenschaftlichen Daten ist die Varianz aber größer als unter Poisson-Verteilung erwartet (z.B. weil nicht gemessene Variablen die abhängige Variable beeinflussen) -> das bezeichnet man als Overdispersion. Wenn die Residual deviance (zweitletzte Zeile in der summary) mehr als 1,5 mal höher ist als die Freiheitsgrade, liegt Overdispersion vor.

summary(glm_model_1)

## 
## Call:
## glm(formula = Anzahl_Blütenendfäule ~ Veredelung * Salzstress, 
##     family = poisson, data = tomaten_daten)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.702  -2.351  -1.222   1.225   5.770  
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    1.01664    0.04093  24.840  < 2e-16 ***
## Veredelungohne                 0.08507    0.05675   1.499    0.134    
## Salzstressohne                 0.27847    0.05414   5.143 2.70e-07 ***
## Veredelungohne:Salzstressohne -0.37364    0.07854  -4.757 1.96e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 4293.7  on 864  degrees of freedom
## Residual deviance: 4256.1  on 861  degrees of freedom
## AIC: 5833.3
## 
## Number of Fisher Scoring iterations: 6

4256/861 ist deutlich größer als 1,5, wir haben in den Daten also tatsächlich Overdispersion.

Lösungen:

• wenn möglich, weitere erklärende Variablen mit in das Modell einbeziehen
• Wenn das nicht hilft/möglich ist: Verwendung von quasipoisson-Fehlern:

glm_model_1 <- glm(Anzahl_Blütenendfäule ~ Veredelung * Salzstress, 
      family = quasipoisson, data = tomaten_daten)

In diesem Modell wird nun zusätzlich geschätzt, wie hoch die Varianz in den einzelnen Behandlungsgruppen ist. Hiermit können wir weiterarbeiten: um die Interaktion zwischen Veredelung und Salzstress auf Signifikanz zu testen, formulieren wir ein einfacheres, genestetes Modell, in dem diese Interaktion fehlt (wenn sich ein Modell B von einem Modell A ableiten lässt, dass heißt, ein oder mehrere erklärende Variablen oder Interaktionen herausgenommen worden sind, sagt man, das Modell B ist in Modell A genestet). Dazu ersetzen wir das * (das für Interaktion steht) mit einem +. So haben beide erklärenden Variablen einen unabhängigen Einfluss auf die Anzahl der Früchte mit Blütenendfäule aber keine Interaktion mehr.

glm_model_2 <- glm(Anzahl_Blütenendfäule ~ Veredelung + Salzstress, 
        family = quasipoisson, data = tomaten_daten)

Mit der Funktion anova() vergleichen Sie die Modelle und testen, ob die Interkation signifikant ist. Zusätzlich müssen Sie hier die Statistik angeben, mit der die Modelle verglichen werden sollen, weil das bei glms mit eindeutig ist. Für quasipoisson-verteilte Residuen, wird der F-test benötigt.

anova(glm_model_1, glm_model_2, test = "F") 

## Analysis of Deviance Table
## 
## Model 1: Anzahl_Blütenendfäule ~ Veredelung * Salzstress
## Model 2: Anzahl_Blütenendfäule ~ Veredelung + Salzstress
##   Resid. Df Resid. Dev Df Deviance      F  Pr(>F)  
## 1       861     4256.1                             
## 2       862     4278.8 -1  -22.723 4.8598 0.02775 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wie schon vermutet, ist die Interaktion signifikant. Wenn eine Interaktion signifikant ist, bleiben auch die beiden erklärenden Variablen, die interagieren auf jeden Fall im Modell. Die Varianzanalyse ist also jetzt abgeschlossen.

Wir könnten uns jetzt mit summary(*Modellname*) noch die geschätzten Parameter-Werte anschauen. Hier müssen Sie aber beachten, dass sie auf der link-Skala angegeben sind (siehe oben, die Parameter werden so angepasst, dass das Modell sinnvolle Ergebnisse liefert) und zurück-transformiert werden müssen, um den eigentlichen Wert zu erhalten. Am einfachsten ist es, die Funktion predict() zu verwenden: Als erstes Argument nennen wir den Namen des glm-Modells, dann erstellen wir einen Datensatz mit den Variablen und Werten, für die wir die Vorhersagen sehen möchten, also mit/ohne Salzstress und mit/ohne Veredelung in allen 4 möglichen Kombinationen. Als type wählen wir ‘response’ aus, weil wir die Schätzungen nicht auf der link-Skala sondern auf der response-Skala (das sind die zurücktransformierten, interpretierbaren Werte, siehe oben) haben wollen.

predict(glm_model_1, data.frame(Salzstress = c("mit", "mit", "ohne", "ohne"), 
                    Veredelung = c("mit", "ohne", "mit", "ohne")), type = "response")

##        1        2        3        4 
## 2.763889 3.009302 3.651376 2.736111

Der geschätzte Wert für die Anzahl an Früchten mit Blütenendfäule für Pflanzen mit Salzstress und mit Veredelung ist also 2,76, für Pflanzen mit Salzstress und ohne Veredelung 3.01 usw.

Link zu Quiz 4:
https://ecampus.uni-bonn.de/goto_ecampus_tst_2823115.html

6.1 Anteilsdaten auswerten

Im Keimungsversuch mit Radis und Bohnen auf zwei verschiedenen Substraten haben Sie Anteilsdaten aufgenommen: wie viele der 80 ausgesäten Samen sind gekeimt? Die csv-Dateien finden Sie hier: https://ecampus.uni-bonn.de/goto_ecampus_file_2786369_download.html

Wie Sie sich schon denken können, erfüllt auch dieser Datentyp nicht die Annahmen, die für eine ‘normale’ ANOVa gegeben sein müssen. Zum Glück können wir wieder auf die generalisierten linearen Modelle zurückgreifen, allerdings mit zwei kleinen Änderungen. Aber fangen wir von vorne an.

Zuerst laden wir das Paket ‘tidyverse’, das uns zum Beispiel das Umstrukturieren von Datensätzen erleichtert, lesen die Daten ein (bei mir sind sie wieder im Ordner ‘data’ innerhalb des Ordners, in dem sich das r-Skript befindet gespeichert) und kontrollieren, ob alles geklappt hat:

library(tidyverse)

Ansaaten <- read.csv("data/Ansaaten.csv", header = T) 
head(Ansaaten)
str(Ansaaten)
summary(Ansaaten)

Um uns die Keimungsraten genauer anzuschauen, erstellen wir einen neuen Datensatz ‘Keimung’, indem wir mit dem pipe-Operator %>% aus dem Datensatz ‘Ansaaten’ die entsprechenden Spalten auswählen:

Keimung <- Ansaaten %>% select(Art, ID, Substrat, Ausfaelle, Angelaufen)
head(Keimung)

##     Art ID  Substrat Ausfaelle Angelaufen
## 1 Bohne  1 Presstopf         5         79
## 2 Bohne  2 Presstopf         3         81
## 3 Bohne  3 Presstopf         6         78
## 4 Bohne  4 Presstopf         6         78
## 5 Bohne  5 Presstopf        21         63
## 6 Bohne  6 Presstopf        23         61

Um einen Säulendiagramm zu erstellen, in dem die Anzahl der angelaufenen Samen und der Ausfälle gestapelt dargestellt werden, müssen wir den Datensatz umstrukturieren, um für jede Anzahl ‘Ausfaelle’ oder ‘Angelaufen’ eine eigene Spalte zu bekommen. Dazu nutzen wir die Funktion ‘reshape’. Wichtig sind hierbei die Parameter

• idvar = die Kombination aus Plot, ID und Behandlungslevels, die eine Zeile in der neuen Datenstruktur definiert
• timevar = Bezeichnung für neue Spalte
• times = das, was in diese Spalte eingetragen werden soll

Keimung_w <- reshape(Keimung, 
        direction = "long",
        varying = list(names(Keimung)[4:5]),
        v.names = "Anzahl",
        idvar = c("Art", "ID", "Substrat"),
        timevar = "Beobachtung",
        times = c("Ausfaelle", "Angelaufen"))
head(Keimung_w)

##                               Art ID  Substrat Beobachtung Anzahl
## Bohne.1.Presstopf.Ausfaelle Bohne  1 Presstopf   Ausfaelle      5
## Bohne.2.Presstopf.Ausfaelle Bohne  2 Presstopf   Ausfaelle      3
## Bohne.3.Presstopf.Ausfaelle Bohne  3 Presstopf   Ausfaelle      6
## Bohne.4.Presstopf.Ausfaelle Bohne  4 Presstopf   Ausfaelle      6
## Bohne.5.Presstopf.Ausfaelle Bohne  5 Presstopf   Ausfaelle     21
## Bohne.6.Presstopf.Ausfaelle Bohne  6 Presstopf   Ausfaelle     23

Jetzt können wir den Plot erstellen

ggplot(Keimung_w, aes(x=factor(ID), y=Anzahl, fill=Beobachtung)) +
  facet_grid(. ~ Art * Substrat) +
geom_bar(position="stack", stat="identity") +
  scale_fill_manual("legend", 
      values = c("Angelaufen" = "darkgreen", "Ausfaelle" = "lightgrey")) +
geom_col(position = position_stack(reverse = TRUE))

Wir sehen, dass die Erfolgsrate bei den Radis unabhängig vom Substrat sehr hoch ist. Bei den Bohnen gibt es mehr Ausfälle, insbesondere bei den auf Steinwolle ausgesäten.

Sind diese Unterschiede signifikant? Um das zu testen, können wir wieder die generalisierten linearen Modelle fitten und mittels der anova-Funktion vergleichen. Bei Anteilsdaten wie diesen erfolgt davor allerdings noch ein Vorbereitungsschritt: Wir generieren aus den beiden Spalten ‘Angelaufen’ und ‘Ausfaelle’ eine einzige abhängige Variable, die wir ‘Keimungsrate’ nennen. Es ist wichtig, die Rate nicht selber auszurechnen: ‘1 von 10 Samen gekeimt’ hat im Modell weniger Gewicht als ‘10 von 100’ Samen gekeimt. Stattdessen verwenden wir die Funktion ‘cbind()’ :

Ansaaten$Keimungsrate <- cbind(Ansaaten$Angelaufen, Ansaaten$Ausfaelle)

Dann können wir die Funktion glm() wie bekannt nutzen. Bei Anteilsdaten nutzen wir als family ‘binomial’

Modell1 <- glm(Keimungsrate ~ Substrat * Art, family = binomial, data = Ansaaten)
summary(Modell1)

## 
## Call:
## glm(formula = Keimungsrate ~ Substrat * Art, family = binomial, 
##     data = Ansaaten)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -7.6302  -1.6393   0.6342   2.5570   7.4933  
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   1.61803    0.09285  17.427  < 2e-16 ***
## SubstratSteinwolle           -1.57041    0.11570 -13.574  < 2e-16 ***
## ArtRadies                     1.61079    0.20275   7.945 1.95e-15 ***
## SubstratSteinwolle:ArtRadies  1.33731    0.26856   4.980 6.37e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1266.76  on 39  degrees of freedom
## Residual deviance:  561.64  on 36  degrees of freedom
## AIC: 705.33
## 
## Number of Fisher Scoring iterations: 5

Die summary des Modells zeigt, dass die Daten overdispersed sind. Auch hier gibt es für diesen Fall die Möglichkeit als family ‘quasibinomial’ anzugeben:

Modell1 <- glm(Keimungsrate ~ Substrat * Art, family = quasibinomial, data = Ansaaten)

Zunächst lassen wir die Interaktion zwischen Substrat und Art weg, indem wir das Multiplikationszeichen durch ein Plus ersetzen

Modell2 <- glm(Keimungsrate ~ Substrat + Art, family = quasibinomial, data = Ansaaten)

Jetzt vergleichen wir die Modelle mittels anaova(). Als Test muss jetzt der F-Test angegeben werden

anova(Modell1, Modell2, test = "F")

## Analysis of Deviance Table
## 
## Model 1: Keimungsrate ~ Substrat * Art
## Model 2: Keimungsrate ~ Substrat + Art
##   Resid. Df Resid. Dev Df Deviance     F Pr(>F)
## 1        36     561.64                         
## 2        37     584.97 -1  -23.326 1.603 0.2136

Wie wir sehen, ist die Interaktion nicht signifikant (p > 0.05) Also vereinfachen wir das Modell weiter und nehmen jetzt das Substrat als erklärende Variable raus

Modell3 <- glm(Keimungsrate ~ Art, family = quasibinomial, data = Ansaaten)
# Wieder vergleichen:
anova(Modell2, Modell3, test = "F")

## Analysis of Deviance Table
## 
## Model 1: Keimungsrate ~ Substrat + Art
## Model 2: Keimungsrate ~ Art
##   Resid. Df Resid. Dev Df Deviance      F   Pr(>F)   
## 1        37     584.97                               
## 2        38     767.96 -1  -182.99 11.644 0.001573 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Das Substrat hat einen hochsignifikanten Effekt auf die Keimungsrate. Als letztes testen wir noch die Art

Modell4 <- glm(Keimungsrate ~ Substrat, family = quasibinomial, data = Ansaaten)

Jetzt müssen wir Modell4 natürlich gegen Modell2 testen. Die Modelle müssen immer ‘genestet’ sein, sonst kann kein Signifikanzniveau berechnet werden. Genestet bedeutet, dass ein Modell eine einfachere Variante (eine oder mehrere Variable/n oder Interaktion/en weniger) des anderen Modells ist.

anova(Modell2, Modell4, test = "F")

## Analysis of Deviance Table
## 
## Model 1: Keimungsrate ~ Substrat + Art
## Model 2: Keimungsrate ~ Substrat
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1        37     584.97                                 
## 2        38    1108.34 -1  -523.38 33.303 1.283e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Auch die Art hat einen hochsignifikanten Einfluss auf die Keimungsrate. Damit ist die Analyse der Keimungsraten abgeschlossen.

6.2 Bonitur-Daten

Bonitur-Daten sind recht speziell, weil es keine metrischen Daten sind, sondern Ordinal-Daten und darüber hinaus auch einigermaßen subjektiv. Trotzdem sind hierfür Methoden zur Auswertung entwickelt worden. Sie stammen aus der Psychologie, in der häufig Untersuchungen stattfinden, bei denen Personen zum Beispiel ihre Zufriedenheit auf einer Skala von 1 bis 5 bewerten sollen.

Wir schauen uns dazu exemplarisch die Bonitur-Daten zur Entwicklung der Bohnen an. Dazu filtern wir den Datensatz Ansaaten, indem wir nur die Zeilen nehmen, in denen das Argument Art == “Bohne” zutrifft (das doppelte Gleichzeichen wird immer benutzt, wenn eine Abfrage gemacht wird, ein einfaches Gleichzeichen ist in R eine Zuweisung also synonym mit dem Pfeil <- ). Mit dem Pipe-Operator %>% verbinden wir diese Auswahl mit der Auswahl der Spalten, die wir mit der Funktion select() machen. Hier können wieder einfach die header (=Spaltenüberschriften) der benötigten Spalten angegeben werden. Vorher laden wir noch zwei Pakete, die für die Auswertung von Ordinal-Daten entwickelt wurden

library(likert)

## Loading required package: xtable

## 
## Attaching package: 'likert'

## The following object is masked from 'package:car':
## 
##     recode

## The following object is masked from 'package:dplyr':
## 
##     recode

library(ordinal)

## 
## Attaching package: 'ordinal'

## The following object is masked from 'package:dplyr':
## 
##     slice

Bohne_Entw_l <- filter(Ansaaten, Art == "Bohne") %>% select(Substrat, ID, Entwicklung)
head(Bohne_Entw_l)

##    Substrat ID Entwicklung
## 1 Presstopf  1           8
## 2 Presstopf  2           5
## 3 Presstopf  3           8
## 4 Presstopf  4           8
## 5 Presstopf  5           4
## 6 Presstopf  6           5

Als erstes sollten die Daten wieder geplottet werden. Bei Bonitur-Daten eignen sich Histogramme am besten. Um das zu tun, müssen die Bonitur-Noten in Faktoren umgewandelt werden. Für die weitere Analyse geben wir auch direkt an, wieviele mögliche Levels (verschiedene Noten) es gibt und dass sie eine aussagekräftige Reihenfolge haben (ordered = TRUE):

Bohne_Entw_l$Entwicklung <- factor(Bohne_Entw_l$Entwicklung, ordered = TRUE, levels = 1:9)

Jetzt können die Histogramme erstellt werden:

ggplot(Bohne_Entw_l, aes(Entwicklung)) +
  geom_bar() +
  facet_wrap(~Substrat, ncol=1) +
  xlab("Bonitur Note") +
  ylab("Anzahl")

Das package ‘likert’ stellt einige Plot-Funktionen zur Verfügung, die speziell für Bonitur- und andere Typen von likert-Daten geeignet sind. Dazu brauchen wir die Daten allerdings im ‘wide’ Format. Im Moment haben wir die Daten im ‘long’ Format: es gibt eine eigene Spalte für das Substrat, sodass ‘Presstopf’ und ‘Steinwolle’ je 10 Mal genannt werden. Wir nutzen wieder die Funktion ‘reshape’ wie oben, jetzt allerdings in umgekehrter Richtung, vom ‘long’ zum ‘wide’ Format:

Bohne_Entw_w <- reshape(Bohne_Entw_l, v.names="Entwicklung", idvar = "ID", 
    timevar = "Substrat", direction="wide")
head(Bohne_Entw_w)

##   ID Entwicklung.Presstopf Entwicklung.Steinwolle
## 1  1                     8                      5
## 2  2                     5                      5
## 3  3                     8                      5
## 4  4                     8                      3
## 5  5                     4                      2
## 6  6                     5                      1

‘idvar’ sind jetzt also die Zeilenbezeichnungen, ‘timevar’ die Spaltenbezeichnungen, ‘v.names’ die eigentlichen Werte.

Als nächstes müssen wir ein sogenanntes ‘likert-Objekt’ aus den beiden Spalten mit den Bonitur-Noten erstellen, dann kann die plot-Funktion des likert-packages genutzt werden:

Bohne_Entw_likert = likert(Bohne_Entw_w[2:3])

plot(Bohne_Entw_likert)

Dies ist eine andere - etwas kompaktere - Darstellung der gleichen Daten. Die %-Werte geben an, wieviel % im schlechten Bereich (1-3), im mittleren Bereich (4-6) und im guten Bereich (7 - 9) waren.

Bonitur-Daten sind Ordinal-Daten und haben damit andere Eigenschaften als metrische Daten. Wie oben erwähnt, erfüllen sie häufig nicht die Annahmen, die für eine normale Varianzanalyse zutreffen müssen. Wir greifen hier auf ‘cumulative link models’ (CLM) aus der library ‘ordered’ zurück. CLMs sind im Prinzip das gleiche wie proportional odds model, falls Ihnen dieser Begriff mal über den Weg läuft. Ähnlich wie in den vorherigen beiden Kapiteln formulieren wir wieder ein komplexeres Modell mit Substrat als erklärender Variable und das Null-Modell, indem angenommen wird, dass alle Bohnen sich ungeachtet der Behandlung gleich Entwickeln:

Model1 <- clm(Entwicklung ~ Substrat, data = Bohne_Entw_l, threshold = "equidistant")

Model2 <- clm(Entwicklung ~ 1, data = Bohne_Entw_l, threshold = "equidistant")

Als threshhold geben wir ‘equidistant’ an, weil wir davon ausgehen, dass zum Beispiel der Unterschied in der Entwicklung zwischen den Bonitur-Noten 2 und 3 genauso groß ist, wie der Unterschied zwischen den Noten 3 und 4 usw.

Auch für CLMs gibt es die Funktion anova(). Hier wird allerdings kein t- oder F-Test durchgeführt, sondern ein Chi-square Test, der uns einen p-Wert liefert:

anova(Model1, Model2)

## Likelihood ratio tests of cumulative link models:
##  
##        formula:               link: threshold: 
## Model2 Entwicklung ~ 1        logit equidistant
## Model1 Entwicklung ~ Substrat logit equidistant
## 
##        no.par    AIC  logLik LR.stat df Pr(>Chisq)   
## Model2      2 78.001 -37.000                         
## Model1      3 73.330 -33.665  6.6711  1   0.009799 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wie sie sehen, ist der p-Wert etwa 0.01, das heißt, das Substrat hat einen signifikanten Effekt auf die Entwicklung der Bohnen. Hiermit ist auch eine grundlegende Analyse der Bonitur-Daten abgeschlossen.

Ein interessanter Link zum Thema, inklusive Auswertungsmöglichkeiten mit R: ‘Do not use averages with Likert scale data’

Link zu Quiz 5: https://ecampus.uni-bonn.de/goto_ecampus_tst_2843090.html

7.1 Zeitreihen plotten

Zunächst laden wir die library ‘ggplot2’, lesen die Wachstumsdaten aus der Datei ‘Braeburn.csv’ ein (im Order ‘data’ auf eCampus) und kontrollieren den Import. ID und Zeitpunkt wandeln wir in Faktoren um. Wie Sie sehen, haben wir Daten von 4 Früchten, 2 von Ästen mit niedrigen Behangdichten und 2 von Ästen mit hohen Behangdichten, deren Umfang stündlich gemessen wurde.

library(ggplot2)

Braeburn <- read.csv("data/Braeburn.csv", header = TRUE)
summary(Braeburn)

##     Datum               Umfang            ID          Behang         
##  Length:2508        Min.   : 1434   Min.   :1.00   Length:2508       
##  Class :character   1st Qu.: 5049   1st Qu.:1.75   Class :character  
##  Mode  :character   Median : 7294   Median :2.50   Mode  :character  
##                     Mean   : 7443   Mean   :2.50                     
##                     3rd Qu.: 9583   3rd Qu.:3.25                     
##                     Max.   :14609   Max.   :4.00                     
##    Zeitpunkt  
##  Min.   :  1  
##  1st Qu.:157  
##  Median :314  
##  Mean   :314  
##  3rd Qu.:471  
##  Max.   :627

Braeburn$ID <- as.factor(Braeburn$ID)
Braeburn$Zeitpunkt <- as.factor(Braeburn$Zeitpunkt)

Um das Wachstum dieser 4 Früchte zu visualisieren, nutzen wir die ggplot Funktion und geben bei den Aestetics (ein Parameter der Funktion) den Zeitpunkt (x-Achse), den Umfang dividiert durch 1000 um von Mikro- auf Millimeter zu kommen (y-Achse) und den Behang (niedrig/hoch) für die Gruppierung durch Farbe an:

ggplot(Braeburn, aes(Zeitpunkt, Umfang/1000, color=Behang)) + 
   geom_point() +   labs(y="Umfang (mm)", x="Zeit") +
   theme_bw()

Im besten Fall, insbesondere wenn wir Unterschiede im Wachstum durch verschiedene Behandlung statistisch nachweisen wollen, haben wir nicht nur 2 Wiederholungen wie hier sondern deutlich mehr. Dies ist im Datensatz ‘Hühner’ der Fall, in dem das Gewicht von Junghennen, die mit unterschiedlichem Futter aufgezogen wurden, über 21 Tage protokolliert wurde.

Auch diesen Datensatz lesen wir ein (auch Ordner ‘data) und wandeln ’Futter’ und ‘ID’ in Faktoren um:

Hühner <- read.csv("data/Huehner.csv", header = TRUE)
head(Hühner)

##   X Gewicht Tag ID Futter
## 1 1      42   0  1      1
## 2 2      51   2  1      1
## 3 3      59   4  1      1
## 4 4      64   6  1      1
## 5 5      76   8  1      1
## 6 6      93  10  1      1

summary(Hühner)

##        X            Gewicht           Tag              ID       
##  Min.   :  1.0   Min.   : 35.0   Min.   : 0.00   Min.   : 1.00  
##  1st Qu.:145.2   1st Qu.: 63.0   1st Qu.: 4.00   1st Qu.:13.00  
##  Median :289.5   Median :103.0   Median :10.00   Median :26.00  
##  Mean   :289.5   Mean   :121.8   Mean   :10.72   Mean   :25.75  
##  3rd Qu.:433.8   3rd Qu.:163.8   3rd Qu.:16.00   3rd Qu.:38.00  
##  Max.   :578.0   Max.   :373.0   Max.   :21.00   Max.   :50.00  
##      Futter     
##  Min.   :1.000  
##  1st Qu.:1.000  
##  Median :2.000  
##  Mean   :2.235  
##  3rd Qu.:3.000  
##  Max.   :4.000

Hühner$Futter <- as.factor(Hühner$Futter)
Hühner$ID <- as.factor(Hühner$ID)

Das entsprechende Diagramm sieht folgendermaßen aus:

ggplot(Hühner, aes(Tag, Gewicht, color=Futter)) + 
   geom_point() +   labs(y="Gewicht (g)", x="Tage seit Geburt") +
   theme_bw()

Für eine bessere Übersicht wäre es wünschenswert, nicht alle Datenpunkte zu plotten¸ sondern den Mittelwert des Gewichts der Hühner einer Behandlung und ein Maß entweder für die Varianz in den Daten oder die Genauigkeit des Mittelwertes. Mit der ggplot-Funktion ‘stat_summary’ können wir genau das machen: der Wert ‘mean_se’ im Parameter ‘fun.data’ bewirkt, dass Mittelwert und Standardfehler des Mittelwertes berechnet und geplottet werden:

ggplot(Hühner, aes(Tag, Gewicht, color=Futter)) + 
  stat_summary(fun.data=mean_se, geom="pointrange") + 
  labs(y="Gewicht (g)", x="Zeit (Tage seit Geburt)")

Wir sehen in diesem Plot einen deutlichen Trend, nämlich dass das Futter des Typs 3 zu schnellerer Gewichtszunahme führt, als die übrigen. Die nächsten beiden Abschnitte erklären, wie wir testen, ob dieser Trend tatsächlich signifikant ist.

7.2 Zufällige Effekte

Im Prinzip sollte man meinen, dass sich Gewichts- oder Wachstumsdaten gut eignen sollten, um eine Varianzanalyse durchführen zu können: sie sind kontinuierlich und wir erwarten im Prinzip eine Normalverteilung der Residuen. Allerdings ist eine grundlegende Annahme verletzt: die Unabhängigkeit der Daten. Wenn wir immer wieder den gleichen Apfel oder das gleiche Huhn messen, können wir davon ausgehen, dass die Werte einander ähnlicher sind, als wenn wir immer andere Indiviuen nehmen würden. Wir haben also keine echten unabhängigen Wiederholungen sondern Pseudoreplikationen. Andere Beispiele die zu Pseudoreplikationen führen sind Versuchsblöcke, Standorte, unterschiedliche Probennehmer, etc.

Wie können solche Faktoren adäquat im Modell berücksichtigt werden? Der klassische Ansatz war lange, sie als ‘normale’ Variablen, also als festen Effekt mit ins Modell aufzunehmen. Das verbraucht aber viele Freiheitsgrade, weil für jeden Level des Zufallseffekts (also jedes Huhn) ein Parameter geschätzt wird (z.B. das mittlere Gewicht jedes Huhns, siehe Abbildung unten). Hier gibt es eine recht gute Erklärung für Freiheitsgrade: https://welt-der-bwl.de/Freiheitsgrade. Je weniger Freiheitsgrade es gibt, desto schwieriger ist es in einer Varianzanalyse, tatsächlich existierende Unterschiede in den Behandlungsgruppen auch als solche zu erkennen (Fähigkeit, Unterschiede zu erkennen = Power der Analyse). Deshalb ist der bessere Ansatz, die Variable als Zufallseffekt einzubeziehen. Das beruht auf der Annahme, dass die Basis-Parameter für jedes Level dieser Variable (jedes Huhn) aus einer gemeinsamen Verteilung stammen. Geschätzt werden muss für das Modell nur die Breite, also Varianz dieser Verteilung. Auf diese Weise wird nur 1 Freiheitsgrad verwendet und die Power der Analyse bleibt höher.

Linke Seite: Klassischer Ansatz - für jedes Huhn wird ein Parameter für den Achsenabschnitt Gewicht geschätzt. Rechte Seite: Nur die breite der Verteilung, aus der die Parameter stammen, wird geschätzt, es wird also nur 1 Freiheitsgrad ‘verbraucht’. Verändert nach einer Folie von Frank Schurr

Feste Effekte

• die geschätzten Paremeter sind von Interesse, es sind zum Beispiel die unterschiedlichen Behandlungslevels (Behangdichte, Futter, …)
• die Namen der Levels haben eine (agrarwissenschaftliche) Bedeutung

Zufällige Effekte

• die geschätzten Parameter interessieren uns nicht
• die Namen der Level haben keine agrarwissenschafliche Bedeutung
• nur einbezogen, um Pseudoreplikation zu vermeiden (ID des Huhn, des Apfels, Block, Ort, Probennehmer)

7.3 Zeitreihen analysieren

Um zufällige Faktoren zu berücksichtigen, benötigen wir gemischte Modelle.

Vergleich der Struktur eines linearen Modells mit einem gemischten linearen Modell, verändert nach einer Folie von Frank Schurr

Ein Paket, das Funktionen für gemischte Modelle zur Verfügung stellt ist ‘lme4’. Dieses Paket laden wir jetzt (nachdem Sie es installiert haben).

library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## 
## Attaching package: 'lme4'

## The following objects are masked from 'package:ordinal':
## 
##     ranef, VarCorr

Die Funktion lmer() fitted lineare gemischte Modelle. Die abhängige Variable und die festen erklärenden Variablen werden wie bekannt in Zusammenhang gebracht: In unserem Beispiel steht das Gewicht links von der Tilde und Tag und Futter werden mit einem Multiplikationszeichen verbunden, um auszudrücken, dass in dem Modell beide Faktoren und ihre Interaktion einen Einfluss auf das Gewicht haben können. Dahinter kommen in runden Klammern die Zufallseffekte: Hier gehen wir davon aus, dass es einen Zufallseffekt geben könnte (individuelle Unterschiede im Gewicht der Hühner), der sich auf den y-Achsenabschnitt der Zeitreihen auswirkt. Das wird ausgedrückt, indem die ID hinter den senkrechten Strich geschrieben wird. Um die Modelle im Anschluss vergleichen zu können, müssen wir noch die Standardmethode ‘REML’ auf ‘FALSE’ setzen.

model1 <- lmer(Gewicht ~ Tag * Futter + (1| ID), data=Hühner, REML = FALSE)

Danach geht es wie gewohnt weiter: wir vereinfachen das Modell, zuerst durch Ersetzen des Multiplikationszeichens mit einem Additionszeichen, womit wir die Interaktion zwischen Tag und Futter rauswerfen. Ist die Interaktion nicht signifikant, nehmen wir als nächstes das Futter ganz raus. Der Zeitpunkt (=Tag) sollte bei solchen Zeitreihen-Analysen allerdings immer im Modell bleiben, es sei denn, wir gehen davon aus, dass es möglicherweise gar keine Veränderung der abhängigen Variable über die Zeit gibt.

model2 <- lmer(Gewicht ~ Tag + Futter + (1 | ID), data=Hühner, REML = F)
model3 <- lmer(Gewicht ~ Tag + (1 | ID), data=Hühner, REML = F)

Dann kann man die Modelle mittels anova() vergleichen. Wie Sie sehen, können auch mehrere Modelle gleichzeitig verglichen werden.

anova(model1, model2, model3)

## Data: Hühner
## Models:
## model3: Gewicht ~ Tag + (1 | ID)
## model2: Gewicht ~ Tag + Futter + (1 | ID)
## model1: Gewicht ~ Tag * Futter + (1 | ID)
##        npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## model3    4 5630.3 5647.8 -2811.2   5622.3                          
## model2    7 5619.2 5649.7 -2802.6   5605.2  17.143  3  0.0006603 ***
## model1   10 5508.0 5551.6 -2744.0   5488.0 117.184  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sie sehen, dass sowohl der Einfluss des Futters auf den Gewichtsverlauf als auch die Interaktion zwischen Futter und Tag hochsignifikant sind. Die Vorhersagen des besten Modells (also model1) können wir nun noch zu dem Plot hinzufügen:

last_plot() + stat_summary(aes(y=fitted(model1)), fun=mean, geom="line")

Und Quiz 6: https://ecampus.uni-bonn.de/goto_ecampus_tst_2853026.html

8.1 Post-hoc Tests - die Problematik

Natürlich können wir jetzt jedes Substrat nochmal einzeln mit einer ANOVA gegen jedes andere Substrat testen und schauen, ob das Ergebnis siginifkant ist. Allerdings gibt es dabei ein Problem: jeder Test hat per Definition eine Irrtumswahrscheinlichkeit von 5% (wir nennen eine erklärende Variable signifikant, wenn der p-Wert, also die Wahrscheinlichkeit solche oder noch extremere Daten zu finden, obwohl die Null-Hypothese wahr ist, kleiner als 0.05 = 5% ist). Je mehr Vergleiche wir machen, desto höher ist also auch die Wahrscheinlichkeit, dass wir in der gesamten Auswertung eine Signifikanz finden, obwohl gar kein wahrer Effekt vorhanden ist - das Problem der multiplen Tests.

source: https://xkcd.com/882/

Einige Statistiker sagen deshalb, dass post-hoc Tests grundsätzlich vermieden werden sollten und auch nicht nötig sind, wenn das experimentelle Design des zugrundeliegenden Versuchs nur gut genug angelegt wurde. Die agrarwissenschaftlichen Realität sieht aber oft anders aus: es wird häufig eine Vielzahl von Behandlungslevels gegeneinander getestet und es interessiert uns nicht nur, OB die Wahl des Substrats einen Effekt auf die Keimungsrate einer Art hat, sondern auch bei welchem davon die Keimungsrate signifikant höher ist. Wahr ist aber, dass man gut durchdenken sollte, wann post-hoc Tests wirklich gebraucht werden und wie genau man sie durchführt.

8.2 Die Herangehensweise

Die generelle Herangehensweise ist, dass man für die Anzahl an Vergleichen, die man macht, korrigiert, indem die einzelnen Signifikanz-Niveaus so angepasst werden, dass man INSGESAMT auf eine Irrtumswahrscheinlichkeit von 5% kommt. Gängige Ansätze sind die Bonferroni-Korrektur und der Tukey’s HSD Test. Hier nutzen wir eine angepasste Bonferroni Korrektur, die sich ‘Holm’ nennt. Wer sich für die methodischen Details interessiert findet hier eine recht gute Erklärung.

Bitte laden Sie die Daten zum Ansaaten-Versuch von diesem Jahr herunter: https://ecampus.uni-bonn.de/goto_ecampus_file_2859385_download.html
und installieren Sie das R-Paket ‘multcomp’.

# Laden der benötigten Pakete
library(multcomp)

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

library(tidyverse)

# Einlesen der Daten und Umwandeln der erklärenden Variablen in Faktoren
ansaaten <- read.csv("data/Ansaaten_neu.csv", header = TRUE)
ansaaten$Art <- as.factor(ansaaten$Art)
ansaaten$Substrat <- as.factor(ansaaten$Substrat)

# Kombinieren der beiden Spalten mit Anteilsdaten
ansaaten$keimung <- cbind(ansaaten$Angelaufen, ansaaten$Ausfälle)

# Umformen des Datensatzes, um einen Plot zu machen
Keimung_w <- reshape(ansaaten, 
                     direction = "long",
                     varying = list(names(ansaaten)[4:5]),
                     v.names = "Anzahl",
                     idvar = c("Art", "Substrat", "Wdh"),
                     timevar = "Beobachtung",
                     times = c("Ausfaelle", "Angelaufen"))
head(Keimung_w)

##                                Art   Substrat Wdh keimung.1 keimung.2
## Bohne.Jiffy.1.Ausfaelle      Bohne      Jiffy   1        63        21
## Bohne.Jiffy.2.Ausfaelle      Bohne      Jiffy   2        34        50
## Bohne.Presstopf.1.Ausfaelle  Bohne  Presstopf   1        76         8
## Bohne.Presstopf.2.Ausfaelle  Bohne  Presstopf   2        76         8
## Bohne.Steinwolle.1.Ausfaelle Bohne Steinwolle   1        64        20
## Bohne.Steinwolle.2.Ausfaelle Bohne Steinwolle   2        70        14
##                              Beobachtung Anzahl
## Bohne.Jiffy.1.Ausfaelle        Ausfaelle     63
## Bohne.Jiffy.2.Ausfaelle        Ausfaelle     34
## Bohne.Presstopf.1.Ausfaelle    Ausfaelle     76
## Bohne.Presstopf.2.Ausfaelle    Ausfaelle     76
## Bohne.Steinwolle.1.Ausfaelle   Ausfaelle     64
## Bohne.Steinwolle.2.Ausfaelle   Ausfaelle     70

# Plot
ggplot(Keimung_w, aes(x=factor(Wdh), y=Anzahl, fill=Beobachtung)) +
  facet_grid(. ~ Art * Substrat) +
  geom_bar(position="stack", stat="identity") +
  scale_fill_manual("legend", 
                    values = c("Angelaufen" = "darkgreen", "Ausfaelle" = "lightgrey")) +
  geom_col(position = position_stack(reverse = TRUE))

# Generalisiertes lineares Modell unserer Hypothese: 
# Substrat und Art haben einen Effekt auf die Keimungsrate
model1 <- glm(keimung ~ Substrat * Art, family=quasibinomial, data = ansaaten)

# Genestetes Modell ohne den Effekt von 'Substrat', 
# um die Signifikanz dieser Variable testen zu können
model2 <- glm(keimung ~ Art, family = quasibinomial, data = ansaaten)

# Varianzanalyse der beiden Modelle
anova(model1, model2, test = "F")

## Analysis of Deviance Table
## 
## Model 1: keimung ~ Substrat * Art
## Model 2: keimung ~ Art
##   Resid. Df Resid. Dev Df Deviance      F    Pr(>F)    
## 1         8      31.78                                 
## 2        14     326.94 -6  -295.16 12.977 0.0009685 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mit der Funktion interaction() generieren wir jetzt eine neue Spalte ‘gruppe’, in der die erklärenden Variablen kombiniert werden. Mit diesen Behandlungsgruppen fitten wir ein neues glm, um sie dann in einem nächsten Schritt in einer post-hoc Analyse gegeneinander testen zu können.

ansaaten$gruppe<-interaction(ansaaten$Art, ansaaten$Substrat)
head(ansaaten)

##     Art   Substrat Wdh Angelaufen Ausfälle keimung.1 keimung.2           gruppe
## 1 Bohne      Jiffy   1         63       21        63        21      Bohne.Jiffy
## 2 Bohne      Jiffy   2         34       50        34        50      Bohne.Jiffy
## 3 Bohne  Presstopf   1         76        8        76         8  Bohne.Presstopf
## 4 Bohne  Presstopf   2         76        8        76         8  Bohne.Presstopf
## 5 Bohne Steinwolle   1         64       20        64        20 Bohne.Steinwolle
## 6 Bohne Steinwolle   2         70       14        70        14 Bohne.Steinwolle

model_posthoc <- glm(keimung ~ gruppe, family=quasibinomial, data = ansaaten)

Bevor wir das tun, sollten wir überlegen, welche Gruppen wir tatsächlich vergleichen wollen. Die Abbildung oben legt nahe, dass Jiffy stripes zu einer höheren Keimungsrate führen als Steinwolle und Sand und wir wollen vielleicht testen, ob diese Unterschiede signifikant sind. In der Funktion glht , die den post-hoc Test für uns durchführt können wir festlegen, welche Gruppen wir gegeneinander testen wollen:

summary(glht(model_posthoc, 
             linfct = mcp(gruppe =
#Ist die Differenz zwischen diesen Gruppen ungleich 0?
            c("(Bohne.Jiffy)-(Bohne.Sand)=0",
              "(Bohne.Steinwolle)-(Bohne.Sand)=0"))),test = adjusted("holm"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: glm(formula = keimung ~ gruppe, family = quasibinomial, data = ansaaten)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error z value Pr(>|z|)  
## (Bohne.Jiffy) - (Bohne.Sand) == 0       -1.4319     0.5202  -2.753   0.0118 *
## (Bohne.Steinwolle) - (Bohne.Sand) == 0  -0.3725     0.5639  -0.661   0.5089  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- holm method)

Jetzt sehen Sie, dass der Unterschied zwischen Jiffy stripes und Sand bei der Bohne signifikant ist, der Unterschied Steinwolle - Sand aber nicht.

Signifikante Unterschiede, die zwischen den Behandlungsgruppen gefunden wurden, werden in Abbildungen häufig durch verschiedene Buchstaben gekennzeichnet. Zum Beispiel so (hier liegen andere Daten zugrunde):

https://r-graph-gallery.com/84-tukey-test.html

Allerdings bietet R keine guten Standardlösungen dafür an. Wenn Sie den TukeyHSD Test durchführen, bietet das Packet ‘multcompView’ eine Möglichkeit. R-Code um die Abbildung oben zu erstellen finden Sie hier.

Test 7 finden Sie hier: https://ecampus.uni-bonn.de/goto_ecampus_tst_2861802.html

9.1 Der Satz von Bayes

Keine Einführung in die Bayesische Statistik kommt ohne den Satz von Bayes aus: Thomas Bayes war ein englischer Geistlicher, Philosoph und Statistiker von dem 1763 - posthum - ein Manuskript mit dem Titel ‘An essay towards solving a problem in the doctrine of chances’ veröffentlicht wurde. In diesem beschreibt er, wie basierend auf vorhandenem Wissen und zusätzlichen Beobachtungen X die Wahrscheinlichkeit dafür berechnet werden kann, dass eine Hypothese H zutrifft:

\[P(H|X) = \frac {P(X|H)×P(H)}{P(X)}\]

In Worten ausgedrückt: Die Wahrscheinlichkeit (P = probability) dass die Hypothese zutrifft, nachdem wir die Daten betrachtet haben (P(H|X), a-posteriori Wissen), ist gleich der Wahrscheinlichkeit, die Daten so zu beobachten, angenommen die Hypothese trifft zu (P(X|H), Likelihood), mal der Wahrscheinlichkeit, dass die Hypothese zutrifft, bevor wir die Daten betrachtet haben (P(H), a-priori Wissen oder Prior), geteilt durch die Wahrscheinlichkeit die Daten so zu beobachten ohne irgendeine Annahme zur Hypothese (P(X), Normalisierungskonstante).
Wichtig ist noch anzumerken, dass Hypothesen dabei als Parameterwerte ausgedrückt werden. Ist unsere Hypothese also, dass unter Behandlung A der Mittelwert der abhängigen Variable größer ist als unter Behandlung B, berechnen wir die Wahrscheinlichkeit für den Parameter Mittelwert(A) - Mittelwert(B) und bekommen als Ergebnis, mit welcher Wahrscheinlichkeit er > 0 ist.

9.2 Das Prinzip des Bayes-Statistik

In der Bayes-Statistik wird diese Idee verwendet, um bereits vorhandenes Wissen - das A-priori-Wissen (P(HP)) - mit den Information, die in den neu erhobenen Daten X stecken zu kombinieren und so das ‘geupdatete- A-posteriori-Wissen zu generieren. Wir sehen außerdem, dass wir als Ergebnis die Wahrscheinlichkeit unserer Hypothese bekommen. Das ist wesentlich leichter bekömmlich und näher an unserem ’normalen’ Denken als die Interpretation eines p-Wertes: ‘Die Wahrscheinlichkeit solche oder extremere Daten zu finden, wenn die Null-Hypothese wahr ist’.

Exkurs: Wahrscheinlichkeitsverteilungen
Um Wahrscheinlichkeiten darzustellen, zum Beispiel das a-priori Wissen, werden Wahrscheinlichkeitsverteilungen verwendet. In der foldenden Abbildung sehen Sie eine solche Wahrscheinlichkeitsverteilung: Je höher die Dichte, desto wahrscheinlicher der Parameterwert, der auf der x-Achse abgetragen ist. Die Gesamtfläche unter der Kurve beträgt immer 1 (genau irgendein Parameterwert trifft immer zu).

Wie berechnen wir also diese Wahrscheinlichkeit? Leider ist das meistens nicht ganz so trivial, wie es aussieht. Zwar können wir meistens die Likelihood berechnen (das wird auch in der frequentistischen Statistik zur Bestimmung der Parameterwerte unserer Modelle gemacht) und unser a-priori Wissen durch eine Wahrscheinlichkeitsverteilung festlegen. Schwierig wird es aber, zumindest sobald wir keine ganz simplen Modelle mehr haben, bei dem Teil P(X), der Wahrscheinlichkeit der Daten. Hierzu müsste man die Wahrscheinlichkeit der Daten über alle überhaupt möglichen Parameterwerte integrieren und das ist oft nicht möglich. Zum Glück kann dieses Problem aber mit einer Simulationstechnik, die in den 80er und 90er Jahren entwickelt wurde, umgangen werden:

9.2.1 MCMC - Markov Chain Monte Carlo

MCMC-Verfahren können Wahrscheinlichkeitsverteilung annähern (in diesem Fall P(H|X)), die man nicht analystisch berechnen kann. Die intuitivste und ziemlich geniale Erklärung, die ich dazu gefunden habe ist die von Michael Franke auf seiner github Seite https://michael-franke.github.io/intro-data-analysis/Ch-03-03-estimation-algorithms.html Ich habe sie hier in Teilen übersetzt:

Image by brgfx on Freepik

Wie die Äpfel an die Bäume kommen
Jedes Jahr im Frühling schickt Mutter Natur die Kinder los, um Äpfel in den Apfelbäumen zu verteilen. Für jeden Baum soll die Anzahl der Äpfel proportional zur Anzahl der Blätter sein: Riese Ralf mit doppelt so vielen Blättern wie Dünner Daniel soll also am Ende auch doppelt so viele Äpfel haben. Wenn es also insgesamt \(n_a\) Äpfel gibt, und \(L(t)\) die Anzahl der Blätter von Baum \(t\) ist, soll jeder Baum \(A(t)\) Äpfel bekommen.

\[A(t) = \frac {L(t)}{\sum L(t')}n_a\]

Das Problem ist nur, dass Mutter Natur \(\sum L(t')\) (die Normalisierungskonstante) nicht ausrechnen kann, also nicht weiß, wie viele Blätter alle Bäume zusammengenommen haben.
Die Kinder (Markov-Ketten) können aber zählen. Sie können zwar nicht alle Bäume besuchen und sich die Zahlen auch nicht so lange merken, aber immerhin für je zwei Bäume. Was sie jetzt tun, ist folgendes: sie starten je an einem zufälligen Baum (Parameterwert), hängen schonmal einen Apfel rein, zählen dort die Blätter \(L(t_1)\) und suchen dann nach einem Baum in der Nähe, von dem sie auch die Blätter zählen \(L(t_2)\) (der Zähler unserer Verteilung). Ist die Zahl der Blätter im zweiten Baum höher, gehen sie dort auf jeden Fall hin und hängen einen Apfel hinr ein. Ist sie niedriger gehen sie nur mit der Wahrscheinlichkeit \(L(t_2)/L(t_1)\) dorthin und hängen einen Apfel auf. So laufen sie weiter zwischen den Bäumen hin und her und am Ende sind die Äpfel ungefähr richtig verteilt. Je mehr Äpfel es sind, desto besser das Ergebnis. Die Besuchshäufigkeit der Kinder hat sich also der gewünschten Verteilung (die Mutter Natur nicht ausrechnen konnte) angenähert!

MCMC Verfahren machen das gleiche: eine MCMC Kette zieht zufällig Werte für die Modell-Parameter und berechnet damit Ergebnis1 = Likelihood der Daten * Prior. Dann zieht sie zufällig Werte, die um die ersten Werte herum liegen, und berechnet dafür wiederum Likelihood der Daten * Prior = Ergebnis2. Wenn Ergebnis2 höher ist als Ergebnis1, springt sie dorthin und zieht von dort aus neue Parameter-Werte. Wenn das Ergebnis2 niedriger ist, springt sie nur mit der Wahrscheinlichkeit Ergebnis2/Ergebnis1 dorthin. Jetzt werden wieder zufällig Werte gezogen usw.
In der folgenden Abbildung sind die erfolgreichen Sprünge durch blaue Pfeile symbolisiert, die abgelehnten Sprünge durch grüne Pfeile.

Visualisierung einer MCMC-Simulation. Sie müssen sich vorstellen, dass Sie die lachsfarbene Parameterlandschaft, die die Werte für Likelihood * Prior darstellt, nicht sehen können, sondern nur die Werte der einzelnen Punkte berechnen können. Weitere Erklärungen im Text. Quelle: https://www.turing.ac.uk/research/research-projects/adaptive-multilevel-mcmc-sampling

Wenn dieser Prozess lang genug fortgesetzt wird, nähert sich die Verteilung der blauen Punkte der a-posteriori Wahrscheinlichkeitsverteilung der Parameter an, die wir haben wollen: Wir haben a-priori Wissen und die Information, die in den neu gesammelten Daten steckt, erfolgreich zum a-posteriori Wissen kombiniert!

9.3 Ein (etwas konstruiertes) Beispiel

Wir müssen uns entscheiden, welche Hühnerrasse wir in größerem Maßstab auf unserer Streuobstwiese halten möchten. Wir haben selbst schon 3 Kraienköppe und 6 Niederrheiner. Die Niederreihner sind uns ein bisschen lieber, weil sie weniger aggressiv sind, aber die Kraienköppe scheinen eine etwas höhere Legeleistung zu haben. Ist es deshalb wirklich sinnvoll, eher Kraienköppe zu wählen?

Schritt 1: Definition der A-priori-Wahrscheinlichkeitsverteilung

Wir erkundigen uns über die Legeleistung der beiden Rassen und erfahren, dass Niederrheiner zwischen 190 und 210 Eier pro Jahr legen, Kraienköppe zwischen 200 und 220. Entsprechend formulieren wir unsere a-priori Wahrscheinlichkeitsverteilungen:

library(brms, warn.conflicts=F, quietly=T)

## Loading 'brms' package (version 2.18.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

# mit der Funktion set_prior definieren wir 2 Normalverteilungen, eine mit Mittelwert 200 
# und Standardabweichung 1,5, und eine mit Mittelwert 210 und Standardabweichung 1,5. 
# Unter 'coef' legen wir fest, zu welchen Parameter im Modell, das wir später 
# fromulieren, diese Priors gehören.
priors <- c(set_prior("normal(200, 5)", coef = "rasseN"), set_prior("normal(210, 5)", coef = "rasseK"))

#So können wir sie plotten, um ein Gefühl dafür zu bekommen
curve(dnorm(x, 200, 5), from = 170, to = 230, xlab="Legeleistung", ylab="Dichte")

2. Erhebung von neuen Daten

Jetzt geben wir unsere eigenen Beobachtungen der Legeleistung der 3 Kraienköppe und 6 Niederrheinern aus dem vorigen Jahr an:

rasse <- c("K", "K", "K", "N", "N", "N", "N", "N", "N")
ll <- c(225, 222, 221, 219, 219, 216, 221, 218, 217)
daten <- data.frame(rasse=rasse, ll=ll)

3. Kombination der a-priori-Verteilung mit den Daten, um die a-posteriori-Verteilung zu erhalten

Um die a-posteriori Verteilung zu schätzen nutzen wir die Funktion ‘brm’ des Packets ‘brms’. Wie wir es schon von anderen Auswertungen kennen, formulieren wir zuerst unser Modell, nämlich das die Legeleistung ll von der Rasse abhängt. Die -1 in der Modell-Formel bewirkt, dass das Modell die Mittelwerte für die Legeleistung schätzt und nicht nur den Mittelwert für die erste Rasse und den Unterschied dazu. Unter ‘data’ geben wir unsere selbst erhobenen Daten ein und unter ‘prior’ die oben definierten Priors. Das argument ‘silent’ bestimmt, wie viele Informationen auf der Konsole ausgegeben werden, wenn die MCMC Ketten loslaufen. Diese Simulationen laufen je nach Modell-Komplexizität einige Sekunden oder Minuten.

legeleistung <- brm(ll ~ rasse -1, data = daten, prior = priors, silent = 2)

4. Interpretation

Der obige Code hat die MCMC Simulation durchgeführt und wir können uns jetzt die a-posteriori Verteilungen für die Legeleistungen der Hühnerrassen anschauen:

# So kann man sich die Zusammenfassung anschauen
summary(legeleistung)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: ll ~ rasse - 1 
##    Data: daten (Number of observations: 9) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rasseK   221.71      1.53   218.12   224.21 1.00     1770     1514
## rasseN   217.61      1.15   214.91   219.46 1.00     1802     1249
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     2.37      0.90     1.32     4.69 1.00     1190     1212
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Die Zusammenfassung zeigt uns zuerst nochmal unsere Eingaben und ein paar Informationen zu den Markov-Ketten. Am interessantesten für uns sind natürlich die Schätzer für rasseK und rasseN und ihre Kredibilitätsintervalle. Wie Sie sehen, überschneiden sich die Konfidenzintervalle: der untere Wert von rasseK (l-95%) ist mit etwa 218 kleiner als der obere Wert (u-95%) von etwa 219 von rasseN, die legeschwächer ist.) Sigma ist die geschätzte Standardabweichung der angenommenen Normalverteilung und interessiert uns hier weniger.

plot(legeleistung)

In der Abbildung sind die Ergebnisse noch leichter zu interpretieren. Links sehen Sie die a-posteriori Verteilungen für beide Rassen (beachten Sie, dass die x-Achsen der beiden oberen Verteilungen nicht ausgerichtet sind). Rechts daneben die Dynamik der Markov-Ketten. Aus den a-posteriori Verteilungen ließe sich jetzt auch die genaue Wahrscheinlichkeit berechnen, dass die Legeleistung wirklich unterschiedlich ist. Da wir aber schon gesehen haben, dass die Kredibilitätsintervalle sich überschneiden, bleiben wir bei unserer Vorliebe für die Niederrheiner, auch wenn es nicht ausgeschlossen ist, dass sie ein paar Eier weniger legen.

9.4 Frequentistisch versus Bayesisch

Die sinnvollere Interpretierbarkeit der Ergebnisse von bayesischen Analysen wird am besten durch folgene Abbildung dargestellt:

Quelle: Korner-Nievergelt und Hüppop (2016). Fünf mögliche Resultate der Schätzung eines Effektes, zum Beispiel der Unterschied im Ertrag bei anderer Düngung. Die Punkte geben die geschätzten Differenzen, die vertikalen Balken die 95 % Unsicherheitsintervalle (Vertrauensintervall bzw. Kredibilitätsintervall) an. Die Resultate des zugehörigen Nullhypothesentest sind in der ersten Zeile aufgeführt. In der zweiten Zeile finden sich die posterior Wahrscheinlichkeit für die Hypothese, dass der Effekt „wirtschaftlich relevant“ ist. Die Hintergrundfarbe unterscheidet schematisch „wirtschaftlich relevante“ (orange) von „wirtschaftlich nicht relevanten“ (hellblau) Effektgrößen.

Bei einem unkritischen frequentistischen Ansatz würde das Ergebnis lauten: “Der Effekt bei Gruppe/Behandlung A und E ist signifikant, also Düngen wir jetzt mit dem getesteten Mittel. Bei den anderen Gruppen gibt es anscheinend keinen Einfluss.”
Natürlich gibt es auch in B einen wichtigen Effekt und auch wenn der Effekt in E statistisch signifikant ist, ist die Effektgröße auf den Ertrag so klein, dass es von wirtschaftlichem Standpunkt aus wahrscheinlich keinen Sinn macht, zusätzlich zu düngen. Die Ergebnisse der bayesischen Analyse (also die Wahrscheinlichkeiten für die Hypothese, dass es einen Effekt gibt), spiegeln diese Interpretationen viel direkter wider.

Auch in diesem Kapitel gibt es eine Quiz-Fragen.

10.1 Ein Thema wählen

Die Bachelorarbeit ist eine fachliche Weichenstellung in Ihrer Ausbildung und Sie werden einige Monate mit dem Anfertigen der Arbeit verbringen. Deshalb ist es wichtig, sich vor der Wahl eines Themas ein paar Gedanken zu machen. Das allerwichtigste Kriterium ist: Was interessiert mich? Wenn Sie ein ernsthaftes Interesse an dem Thema Ihrer Arbeit haben, vielleicht sogar begeistert davon sind, ist das die beste Voraussetzung für einen guten Verlauf der Arbeit und ein überzeugendes Ergebnis.

Weitere Entscheidungskriterien sind:

• Welche Methoden will ich lernen?
• In welchem Bereich habe ich schon Vorwissen und kann es in die Arbeit einbringen?
• Was ist meine berufliche Zielsetzung? Möchten ich die Arbeit in Zusammenarbeit mit einem Betrieb schreiben? (In diesem Fall werden häufig Themen vorgeschlagen)
• Gibt es Arbeitsgruppen / Betreuuer*innen, die Experten auf dem Gebiet sind, und mich deshalb mit Enthusiasmus unterstützen?

Wo Sie Themen finden:

• Uni-Homepage: Forschungsthemen und Publikationen der Arbeitsgruppen anschauen,
• Uni-Homepage: explizite Auflistung der Bachelor-Themen, zum Beispiel:

https://www.gartenbauwissenschaft.uni-bonn.de/teaching/

https://www.aol.uni-bonn.de/de/lehre/angebot-von-bachelor-und-masterarbeitsthemen-2018

• Mit Lehrenden sprechen, natürlich insbesondere wenn die Module Sie interessiert haben
• Mit Mitstudierenden sprechen

10.2 Literatur recherchieren

Wenn Sie ein Thema und eineN BetreuerIn gefunden haben, sollten Sie sich genügend Zeit nehmen, das Wissen, das zu dem Thema schon vorhanden ist, zu recherchieren. Nur so können Sie sicherstellen, dass Sie nicht stumpf wiederholen, was schon gemacht worden ist, sondern etwas Neues herausfinden, das das Forschungsfeld oder die praktische Nutzung voran bringt. Leider fehlt auch den engagiertesten BetreuerInnen häufig die Zeit, immer über alle neuen Publikationen zu den angebotenen Themen informiert zu sein, deshalb ist es wichtig hier gründlich zu arbeiten. Nichts ist ärgerlicher als am Ende einer Arbeit festzustellen, dass genau das gleiche schon gemacht wurde, insbesondere wenn Sie die Arbeit publizieren wollen.

Üblicherweise wird der/die BetreuerIn Ihnen als Einstieg ein paar Referenzen zu wichtigen Publikationen aus dem Bereich geben. Diese Publikationen können Sie als ‘Anker’ verwenden: Welche Literatur wird zitiert (Referenzliste anschauen)? Welche neueren Publikation zitieren dieses Paper (Kann in Suchmaschinen wie dem Web of Science angezeigt werden)?

Wenn Sie einige Publikationen gelesen haben, kennen Sie die relevanten Fachbegriffe und den verwendeten Wissenschaftsjargon. Jetzt ist es sinnvoll, sich ein Liste mit keywords anzulegen, die Sie für die Suche mit den Suchmaschinen nutzen. Häufig werden bei den Publikationen auch explizit keywords angegeben, die Sie verwenden können. Machen Sie sich mit den Such-Operatoren der Suchmaschienen vertraut (z.B. AND, OR, NOT, siehe http://images.webofknowledge.com/images/help/WOS/hs_search_operators.html). Hier eine Auflistung von Suchmaschinen:

• Web of Science (Clarivate): http://www.webofknowledge.com
• Google Scholar (https://scholar.google.com/, ggf. Advanced Search nutzen)
• ULB Bonn, Bonnus: http://www.ulb.uni-bonn.de
• Science Direct (Elsevier): http://www.sciencedirect.com
• Researchgate: http://www.researchgate.net
• Statistische Datensätze: Statistisches Bundesamt Deutschland (http://www.destatis.de)
• JSTOR (Journal STORage): http://www.jstor.org
• Scopus (Elsevier): http://www.scopus.com
• FAO: http://www.fao.org

Je nach Thema, werden Sie kaum Publikationen finden oder damit überschwemmt werden. Deshalb hier ein paar Such-Taktiken um die Suche auszuweiten oder einzugrenzen:

• mehrere Suchmaschinen mit denselben Suchbegriffen verwenden
• Allgemeine und spezifische Suchbegriffe vergleichen, z.B. Fruit trees ↔︎ apple trees ↔︎ Malus sylvestris
• bei guten Publikationen gezielt nach anderen Veröffentlichungen der Autoren suchen, sowie nach zitierten und zitierenden Publikationen
• möglichst aktuelle Publikationen suchen, Publikationszeitraum ggf. eingrenzen
• qualitativ hochwertige Journals bevorzugen (-> Impact factor)
• gibt es Publikationen des Betreuers/der Betreuerin oder der Arbeitsgruppe. Die sollten Sie natürlich nicht vergessen :)

Impact Factor IF = SC / SA
SC: Summe der Zitationen in dem Jahr von Interesse, die Artikel der letzten 2 Jahre im selben Journal zitieren.
SA: Summe der Artikel, die in den letzten 2 Jahren in dem Journal publiziert wurden

In den Suchmaschinen werden Sie in manchen Fällen direkt die PDF-Datei mit der gesamten Veröffentlichung finden (wenn die Artikel frei verfügbar sind, oder die Uni Bonn ein Abo für das Journal hat). Manchmal finden Sie aber auch nur den Titel und ein paar weitere Meta-Daten. Was dann?

• Google suche nach Titel in Anführungszeichen

  • direkt zum PDF
  • Online Repository (z.B. Researchgate, BioRxiv)

• Kontakt zu Autoren (corresponding author, z.B. über Research Gate oder per E-Mail)

• KollegInnen fragen, die Zugang haben

• Andere Wege…

Etwas Werbung für Preprint-Server
Preprint-Server beherbergen Datenbanken mit Manuskripten zu einem bestimmten Themenkomplex, die ohne Review Prozess und ohne spezielle Formatierung dort abgelegt und veröffentlicht werden. Eine eigene kritische Begutachtung ist deshalb zwingend notwendig. Organisiert werden solche Server von non-profit Organisation, die der Idee des Rechts auf freien Zugang zu Wissen folgen. Die Manuskripte sind zitierfähig und können trotzdem bei den meisten wissenschaftlichen Journalen (nochmal) veröffentlicht werden. Vorbild ist der Preprint-Server arxiv.org (englisch ausgesprochen archive-dot-org) der hauptsächlich in der Mathematik und Physik genutzt wird. Mittlerweile gibt es aber auch BioRxiv für die Biowissenschaften und Agrirxiv für die Agrarwissenschaften https://agrirxiv.org/about/

Egal wie sorgfältig Sie suchen, Sie werden nicht alles vorhandene Wissen zu dem Thema finden können, schlicht weil es aus unterschiedlichen Gründen nicht verfügbar ist, wie das folgende Bild veranschaulicht. Zu ändern ist das kaum, aber es ist gut, sich diser Einschränkung bewusst zu sein.

10.3 Literatur verwalten

All die Literatur, die Sie gefunden und gelesen haben, muss verwaltet werden, hauptsächlich deshalb, weil Sie sie in Ihrer Arbeit zitieren werden. Dazu ist es sehr, sehr sinnvoll, ein Literaturverwaltungsprogramm zu verwenden. Am besten wirklich gleich von Anfang an, das erspart Ihnen viel Arbeit!

Mit solchen Programmen können Sie Literatur sammeln (einschließlich der PDFs), Literatur teilen (gemeinsame Ordner), Referenzen in den Text einfügen, automatisch ein Literaturverzeichnis erstellen und auch automatisch den Reference Style auswählen und verändern. Häufig verwendete Programme sind

• Mendeley (https://www.mendeley.com, kostenlos)
• Zotero (https://www.zotero.org, kostenlos)
• Citavi (https://www.citavi.com/de, kostenlose Basisversion)
• Endnote (https://endnote.com, kostenlose Testversion)

In unserer Arbeitsgruppe haben wir uns für Zotero entschieden, weil es eine gute Funktionalität besitzt und kostenfrei ist. Hier ein Video zur Benutzung von Zotero (noch eine Anmerkung zum Video: Zotero kann nicht nur mit Chrome, wie hier gezeigt, sondern auch mit anderen Internet-Browsern verwendet werden, es gibt mittlerweile ein GoogleDocs PlugIn und eine Kompatibilität mit dem Editor R-Studio):

Youtube Video zu Zotero

10.4 Eine gute wissenschaftliche Frage formulieren

Wenn Sie nach der Literaturrecherche eine guten Überblick über das Thema und vor allem auch die noch vorhandenen Wissenslücken haben, sollten Sie in der Lage sein, zusammen mit Ihrem Betreuer/Ihrer Betreuerin eine gute wissenschaftliche Frage formulieren. Warum ist das so wichtig? Die beste Antwort darauf gibt Douglas Adams in seinem Buch ‘The Hitchhiker’s guide to the Galaxy’:

Der Computer Deep Thought hatte siebeneinhalb Millionen Jahre gerechnet, um jetzt eine Antwort auf die Große Frage zu geben…

“Both of the men had been trained for this moment, their lives had been a preparation for it, they had been selected at birth as those who would witness the answer, but even so they found themselves gasping and squirming like excited children.
”And you’re ready to give it to us?” urged Loonsuawl.
“I am.”
“Now?”
“Now,” said Deep Thought.
They both licked their dry lips.
“Though I don’t think,” added Deep Thought. “that you’re going to like it.”
“Doesn’t matter!” said Phouchg. “We must know it! Now!”
“Now?” inquired Deep Thought.
“Yes! Now…”
“All right,” said the computer, and settled into silence again. The two men fidgeted. The tension was unbearable.
“You’re really not going to like it,” observed Deep Thought.
“Tell us!”
“All right,” said Deep Thought. “The Answer to the Great Question…”
“Yes..!”
“Of Life, the Universe and Everything…” said Deep Thought.
“Yes…!”
“Is…” said Deep Thought, and paused.
“Yes…!”
“Is…”
“Yes…!!!…?”
“Forty-two,” said Deep Thought, with infinite majesty and calm.”
“Forty-two!” yelled Loonquawl. “Is that all you’ve got to show for seven and a half million years’ work?”
“I checked it very thoroughly,” said the computer, “and that quite definitely is the answer. I think the problem, to be quite honest with you, is that you’ve never actually known what the question is.”

— Douglas Adams

Ohne konkrete Frage - keine hilfreiche Antwort! Die Frage muss also inhaltlich so konkret wie möglich und sprachlich so präzise wie möglich formuliert werden. Aus der Frage wird das experimentelle Design abgeleitet (manchmal auch erst eine Hypothese), das der Arbeit zugrunde liegt und zu neuen Erkenntnissen führen soll. Darüber hinaus sollte die Beantwortung der Frage relevant für die Wissenschaft und/oder Praxis sein (das sollte der/die BetreuerIn einschätzen) und zudem sollte sie die Arbeit klar eingrenzen, so dass die Beantwortung zeitlich zu schaffen ist.

Beispiel für eine viel zu unspezifische Frage:
Ist es gut, Apfelbäume zu düngen?

Besser:

• Hat die Düngung mit x einen positiven Effekt auf den Gesamtertrag (kg Früchte) der Apfel-Sorte a?
• Welche Düngungsverfahren xyz zeigt den stärksten positiven Effekt auf den Gesamtertrag (kg Früchte) der Apfel-Sorte a?
• Hat die Düngung mittels x einen positiven Effekt auf den Gesamtertrag (kg Früchte) auf alle untersuchten Apfel-Sorten abc? Bei welcher Sorte ist der stärkste Effekt zu beobachten?
• Ist der Effekt von Düngung mit unterschiedlichen Mitteln xyz abhängig von der untersuchten Apfel-Sorte abc?

Hier der Link zu Quiz 8: https://ecampus.uni-bonn.de/goto_ecampus_tst_2890306.html

11.1 Prinzipien für ein gutes Design

Bevor Sie einen Versuch im Freiland oder Gewächshaus starten, um die wissenschaftliche Frage, die Sie nun formuliert haben, zu beantworten, ist es sehr wichtig, genau zu überlegen, wie der Versuch aufgebaut sein soll. Dabei bekommen Sie natürlich Unterstützung von Ihrer Betreuerin/Ihrem Betreuer. Je nachdem wie kompliziert das Design aussehen wird und wie gut sich der Betreuer/die Betreuerin selbst mit statistischen Auswertungen auskennt, kann es auch sinnvoll sein, sich schon zu diesem Zeitpunkt Rat von Statistik-Experten einzuholen. Das ist viel besser als später mit einem löchrigen Datensatz und phantasievollem Design zu kommen und zu hoffen, dass die Analyse durch irgendwelche Tricks noch gerettet werden kann. Die Devise ist: “Erst denken, dann handeln’. Wenn Sie ganz sicher gehen wollen, ist es ein weiterer guter - wenn auch etwas seltsam klingender - Rat: erfinden Sie erstmal Ihre Daten so wie Sie sie in dem angedachten Versuchsdesign erheben könnten. Bei der Auswertung dieser Pseudo-Daten wird häufig schnell klar, wo Schwierigkeiten auftreten könnten und wie Sie gegensteuern können. Natürlich müssen Sie die erfundenen Daten nachher sehr gewissenhaft durch die echten ersetzen!

Das Wichtigste ist, dass Sie mit Ihrem Versuch die Fragen, die Sie formuliert haben, klar beantworten können. Beim Versuchsdesign und auch später bei der Auswertung der Daten ist es oft sinnvoll, sich nochmal auf diese Fragen zurück zu besinnen. Anhand der Fragen entscheiden Sie, welche Pflanzen Sie verwenden, welche Behandlungsmethoden Sie anwenden und welche Daten Sie erheben. Sind Sie nur an der Zahl marktfähiger Früchte pro Pflanze interessiert, ist es normalerweise unnötig, über die gesamt Saison Wachstumsparameter der Pflanzen zu erheben. Einige zusätzliche Messungen können allerdings sinnvoll sein, um die gefundenen Ergebnisse interpretieren und diskutieren zu können.

Außer den oben genannten Aspekten gibt es eine ganze Reihe statistischer Grundprinzipien für ein gelungenes Design, die eine unkomplizierte und belastbare Analyse ermöglichen:

1. Repräsentativität der Stichprobe Mit unserem Experiment wollen wir eine allgemeine Aussage über eine bestimmte Grundgesamtheit, z.B. Apfelbäume der Sorte Pinova und deren Raktion auf eine bestimmte Behandlung machen. Da wir nicht alle existierenden Apfelbäume der Sorte Pinova beproben können, müssen wir mit einer Stichprobe, also einer Teilmenge aller vorhandenen Apfelbäume, zurecht kommen. Diese Stichprobe sollte repräsentativ sein, das heißt zum einen die gesamte Bandbreite an möglichen Werten (z.B. Ertrag) abdecken und zum anderen nicht verzerrt sein. Das können wir erreichen, indem wir eine ausreichend große Stichprobe erheben (d.h. genügend Wiederholungen, Replikationen, im statistischen Jargon meistens mit N bezeichnet) und sicherstellen, dass diese Wiederholungen zufällig gewählt werden. Beprobt man nur Apfelbäume der Sorte Pinova mit der Unterlage M18, dann ist die Stichprobe nur für diese Unterlage repräsentativ, nicht aber für alle Apfelbäume der Sorte Pinova.
Wichtig ist, dass Sie für jede Wiederholung auch einen eigenen Wert erheben:

Wenn Sie zum Beispiel die Früchte von 10 Erdbeerpflanzen, die einer bestimmten Behandlung ausgesetzt sind, zusammenwerfen und den Fruchtzuckergehalt bestimmen, bekommen Sie zwar eine bessere Schätzung vom mittleren Fruchtzuckergehalt der Früchte unter dieser Behandlung, als wenn Sie dazu nur die Früchte einer einzelnen Pflanze verwenden. Was aber dabei verloren geht ist das Wissen über die Variabilität, also die Varianz zwischen den Pflanzen, und die benötigen wir für die statistische Auswertung der Ergebnisse. Ideal wäre, den Fruchtzuckergehalt pro Pflanzen zu bestimmen. Wenn das zu aufwändig/kostenintensiv ist, ist es immernoch besser, zum Beispiel nur 4 Pflanzen anstelle von 10 zu verwenden, diese dann aber getrennt voneinander auf den Fruchtzuckergehalt zu analysieren, um 4 echte Wiederholungen zu generieren.

2: Unabhängigkeit Die einzelnen Wiederholungen müssen unabhängig voneinander sein, was bedeutet, dass die Werte einer Wiederholung keine Informationen über die Werte anderer Wiederholungen enthalten, es also keine Verbindung zwischen den Stichproben gibt.

Wenn Sie zum Beispiel 3 Apfelbäume je einer Behandlung unterziehen und danach pro Baum die Fruchtsäure von 10 Früchten messen, sind die Werte der 10 Äpfel keine unabhängigen Wiederholungen, weil die gemessenen Werte von Früchten desselben Baumes auch ohne die Behandlung ähnlicher wären, als die Werte von Früchten verschiedener Bäume. Der Fruchtsäuregehalt des ersten gemessenen Apfels gibt also Informationen darüber, wie hoch der Fruchtsäuregehalt des 2. Apfels desselben Baumes in etwa ist. Man spricht in diesem Fall von auch Pseudo-Replikationen. Wenn solch ein Design nicht zu vermeiden ist (weil keine riesige Zahl von Apfelbäumen zur Vergügung steht), gibt es aber Möglichkeiten, in der statistischen Analyse für solche Abhänigkeiten zu korrigieren. Die Analyse hat dann allerdings weniger ‘Kraft’ echte Unterschiede zwischen den Behandlungen aufzudecken.

3: Kontrolle Um den Effekt einer Behandlung quantifizieren zu können, müssen wir wissen, wie das Ergebnis ohne die Behandlung ausgesehen hätte. Dazu brauchen wir eine Kontrollgruppe, also Versuchseinheiten, die nicht mit der zu untersuchenden Methode behandelt wurden. Damit wir den eigentlichen Effekt der Behandlung messen können, müssen wir mögliche Nebeneffekte der Behandlung bedenken und die Kontrolle entsprechend anpassen.

Wenn wir zum Beispiel Tomaten wöchentlich mit einem in Wasser gelösten Pestizid besprühen, kann nicht nur das Pestizid an sich, sondern auch das Besprühen einen Effekt auf die Fruchtqualität haben. Deshalb sollte die Kontrollgruppe ebenfalls besprüht werden, in diesem Fall mit reinem Wasser.

4: Orthogonalität Häufig untersuchen wir nicht nur den Effekt einer Behandlung sondern interessieren uns für die Effekte mehrerer Behandlungen und deren Kombination. Zum Beispiel möchten wir die Reaktion von Getreide auf Phosphatdüngung und Stickstoffdüngung untersuchen. Wenn wir nur 3 Typen von Untersuchungsflächen haben (Kontrolle, Phosphor, Phospor + Stickstoff), können wir den Effekt von Phosphor + Stickstoff nicht eindeutig interpretieren, weil wir nicht wissen, wie die Reaktion auf reine reine Stickstoffdüngung gewesen wäre. Nur ein faktorielles Design in dem alle Level der einen Behandlungen mit allen Leveln der anderen Behandlung kombiniert werden, erlaubt eine vollständige Interpretation. Wir brauchen bei diesem Experiment also die Behandlungen Kontrolle, nur Phosphor, nur Stickstoff, Phosphor und Stickstoff. So sind die Behandlungen Phosphor und Stickstoff in der Analyse unabhängig voneinander (= orthogonal) und wir können auf die einzelnen Effekte von Phosphor und Stickstoff und auf mögliche Interaktionen (Hat die Verfügbarkeit von Phosphat einen Einfluss auf die Reaktion der Pflanze auf Stickstoff?) testen.

5: Konstante Bedingungen Außer den Effekten, auf die wir testen wollen, sollten die Umweltbedingungen für alle Versuchseinheiten möglichst konstant sein. Wir müssen zum Beispiel davon ausgehen, dass Sonneneinstrahlung in nur einem Teil des Gewächshauses einen Effekt auf verschiedene Pflanzenparameter (Chlorophyll, Wachstum, Zuckergehalt) hat. Ebenso kann die Bodenqualität über die Versuchsfläche variieren und ebenfalls einen Effekt auf Wachstum der Pflanzen etc. haben. Wenn möglich sollten solche Effekte vermieden werden. In der Praxis ist das allerdings häufig nicht möglich. Das wichtigste ist zu verhindern, dass nur bestimmte Level einer Behandlung (z.B. nur die mit Stickstoff gedüngten Pflanzen im Beispiel oben) im helleren Teil des Gewächshauses stehen. So kann man später nicht mehr herausrechnen, ob der beobachtete Effekt eine Reaktion auf die Düngung oder die stärkere Sonneneinstrahlung ist. Deshalb sollten die Proben immer zufällig durchmischt stehen, sodass es in allen Ecken des Gewächshauses Pflanzen mit allen Behandlungsleveln gibt. Denken Sie auch daran, dass es einen Einfluss auf die Werte haben kann, wer die Beprobung vorgenommen hat. Wenn Sie mit mehreren Personen Daten erheben, sollte nicht eine Person alle Kontrollen aufnehmen und die andere Person alle anderen Behandlungen. Schreiben Sie im Zweifel auf, wer welche Daten erhoben hat. Auch dieser Faktor kann bei der statistischen Analyse einbezogen werden.

Eine weitere Möglichkeit, Umweltheterogenität herauszurechen besteht in der Bildung von Blöcken. Wenn wir zum Beispiel in 8 Parzellen unsere Behandlungslevels wiederholen, stellt die Parzelle einen Block dar und wir nutzen Parzellennummer, um mögliche (unbekannte) Unterschiede zwischen den Parzellen aus den Ergebnissen herauszurechnen. Die für das Versuchsziel wichtigen Vergleiche müssen dafür innerhalb der Blöcke vorgenommen werden. Zweck der Blockbildung ist es, die Genauigkeit blockinterner Vergleiche zu erhöhen.

Regel 6: Balanciertes Design Ein balanciertes Desging bedeutet, dass in allen Behandlungsgruppen die gleiche Anzahl von Wiederholungen sind. Viele statistische Tests setzen balancierte Designs voraus. Einzelne fehlende Werte (z.B. wenn Pflanzen gestorben sind) sind aber kein Problem.

Das ideale Versuchsdesign wird häufig aufgrund von mangelndem Platz, mangelnder Zeit, mangeldem Geld… nicht möglich sein. Das ist auch allen Betreuern und auch Reviewern und Editoren von wissenschaftlichen Journalen auch klar. Wichtig ist, diese Beschränkungen in der Arbeit zu nennen und die möglichen Effekte zu diskutieren.

11.2 Wichtige Designs für manipulative Experimente

Es gibt einige typische Designs für manipulative Experimente, welche die obigen Regeln beachten und vernünftig ausgewertet werden können. Hier die drei wichtigsten:

11.3 Stichprobengröße abschätzen

Grundsätzlich gilt, dass man mit mehr Wiederholungen die Effekte einer Behandlung sicherer nachweisen kann. Natürlich ist die Anzahl an Wiederholungen meistens durch den zeitlichen, räumlichen und finanziellen Rahmen begrenzt. Wenn man etwas Vorwissen zu der Variabilität der Parameter, die gemessen werden sollen (Ertrag, Größe, etc) hat, kann man abschätzen, wie viele Wiederholungen benötigt werden, um einen Effekt einer bestimmten Größe nachzuweisen.

Dabei gilt: je größer die Reaktion auf die Behandlung ist, je weniger variable die gemessenen Parameter sind und je größer die Stichprobe, desto wahrscheinlicher lässt sich der Effekt nachweisen. Dies ist in folgender Abbildung veranschaulicht: die blauen und grünen Kurven stellen gemessene Werte (entlang der x-Achse) mit zwei verschiedenen Behandlungen dar. Die grüne Behandlung führt zu höheren Werten (z.B. Ertrag), sie liegen auf der x-Achste weiter rechts. Liegen die Kurven nah beieinander, wie im Beispiel auf der linken Seite, ist der Effekt weniger leicht nachzuweisen als wenn sie weiter voneinander entfernt liegen, wie im Beispiel rechts oben. Der Effekt ist auch schwerer nachzuweisen, wenn die Kurven breit sind, es also viel Variabilität in den Daten gibt. Wenn sich die aufgenommenen Werte innerhabl der Behandlungen stark ähneln, also nicht so variabel sind, sind die Kurven schmaler und der Effekt der Behandlung kann leichter nachgewiesen werden. Der Effekt kann auch leichter nachgewiesen werden, wenn es viele Wiederholungen gibt, die Lage und Form der Kurven also sicherer bestimmt werden kann, wie im Beispiel unten rechts im Gegensatz zum Beispiel unten links.

In R gibt es die Möglichkeit, die minimal nötige Anzahl an Wiederholungen zu schätzen, um einen Effekt einer bestimmten Größe nachweisen zu können. Dafür benötigt man

• d = Unterschied/Effekt, den man mindestens erkennen können möchte (detection level)
• α = Signifikanzniveau, das festlegt, ab wann man ein Ergebnis als signifikant einstuft (üblicherweise 0.05)
• power = gewünschte Trennschärfe des Testes, üblicherweise 0.9
• Varianz innerhalb der Behandlungen

Frage: Unterscheidet sich der Ertrag bei Himbeeren unter 4 verschiedenen Behandlungen?
Signifikanzniveau: α = 0.05
Trennschärfe: power = 0.9
Vorwissen: Varianz innerhalb der Gruppen ist in etwa 280.
Wir wollen einen Effekt nachweisen können wenn er mindestens 400 g beträgt.

In wie vielen Parzellen mit Himbeeren muss der Ertrag gemessen werden, um den Effekt der Behandlung statistisch nachweisen zu können?

power.anova.test(groups=4, n=NULL, between.var=400,
within.var=280, sig.level=0.05, power=0.9)

## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 4
##               n = 4.429616
##     between.var = 400
##      within.var = 280
##       sig.level = 0.05
##           power = 0.9
## 
## NOTE: n is number in each group

R hat berechnet, dass es mindestens n = 4.42… Wiederholungen braucht, um einen Effekt von mindestens 400g nachweisen zu können. Da wir nur mit ganzen Parzellen arbeiten, sollten Sie mindestens 5 Wiederholungen pro Behandlung einplanen.

11.4 Weitere Überlegungen zum Versuchsdesign

Außer den genannten statistischen Erwägungen sollte natürlich ein gartenbauliches Verständnis die Planung des Versuches leiten. Wichtig hierbei sind zum Beispiel Effekte die durch die Lage einer Wiederholung im Versuchsaufbau zustande kommt (Randeffekte) sowie die Praktikabilität der Datenerhebung.

Und der Link zum Quiz: https://ecampus.uni-bonn.de/goto_ecampus_tst_2896163.html