Chapter 4 Zähl- und Boniturdaten auswerten
In der letzten Woche haben Sie gelernt, wie man Experimente mit kontinuierlichen Daten der abhängigen Variable auswertet. In den Agrarwissenschaften kommen aber auch häufig andere Datentypen vor. In diesem Kapitel beschäftigen wir uns speziell mit Zähldaten (wieviele Äpfel trägt der Baum? wieviele Läuse sitzen auf der Pflanze? wieviele Eier legt das Huhn?)
Am Ende dieses Kapitels
- wissen Sie, warum Zähldaten häufig nicht mit einer normalen Varianzanalyse ausgewertet werden können.
- kennen Sie eine Methode, um solche Daten trotzdem zu analysieren.
- sind Sie eine Fallstudie nochmal im Detail durchgegangen.
4.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).
4.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).
4.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.
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))\).
4.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.
4.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.
##
## Call:
## glm(formula = Anzahl_Blütenendfäule ~ Veredelung * Salzstress,
## family = poisson, data = tomaten_daten)
##
## 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: 64256/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.
## 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 ' ' 1Wie 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.736111Der 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.
4.6 Anteils- und Boniturdaten auswerten
In dieser Woche geht es gleich um zwei Datentypen: Anteils- und Bonitur-Daten. Anteilsdaten können mit generalisierten linearen Modellen (GLMs) ausgewertet werden, die Sie beim letzten Mal schon kennengelernt haben - allerdings mit zwei kleinen Änderungen. Bonitur-Daten sind im Gartenbau ziemlich häufig, sind allerdings in Bezug auf die Auswertung ein eher undankbarer Datentyp: oft ist es schwierig bis unmöglich, sie durch Verteilungen wie der Normalverteilung zu charakterisieren. Deshalb wiedersetzen sie sich auch den Auswertungen, die auf solchen Verteilungen beruhen.
Am Ende dieses Kapitels wissen Sie
- wie die abhängige Variable in einem GLM zu Anteilsdaten generiert wird
- mit welchem Modell und welchem Test Sie Anteilsdaten auswerten, um signifikante Unterschiede zwischen den Behandlungen zu identifizieren
- wie man die etwas störrischen Bonitur-Daten zähmen, plotten und auswerten kann.
4.7 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:
## 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 61Um 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 23Jetzt 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()’ :
Dann können wir die Funktion glm() wie bekannt nutzen. Bei Anteilsdaten nutzen wir als family ‘binomial’
##
## Call:
## glm(formula = Keimungsrate ~ Substrat * Art, family = binomial,
## data = Ansaaten)
##
## 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: 5Die summary des Modells zeigt, dass die Daten overdispersed sind. Auch hier gibt es für diesen Fall die Möglichkeit als family ‘quasibinomial’ anzugeben:
Zunächst lassen wir die Interaktion zwischen Substrat und Art weg, indem wir das Multiplikationszeichen durch ein Plus ersetzen
Jetzt vergleichen wir die Modelle mittels anaova(). Als Test muss jetzt der F-Test angegeben werden
## 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.2136Wie 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 ' ' 1Das Substrat hat einen hochsignifikanten Effekt auf die Keimungsrate. Als letztes testen wir noch die Art
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.
## 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 ' ' 1Auch die Art hat einen hochsignifikanten Einfluss auf die Keimungsrate. Damit ist die Analyse der Keimungsraten abgeschlossen.
4.8 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
## 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##
## Attaching package: 'ordinal'## The following object is masked from 'package:dplyr':
##
## sliceBohne_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 5Als 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):
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:
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:
## 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 ' ' 1Wie 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’