Datensatz

Hier wurde ein Experiment mit Mohrrüben durchgeführt ( Jansen und Hoekstra, 1993). In einer randomisierten vollständigen Blockanlage wurden 16 Sorten bei zwei verschiedenen Insektizidbehandlungen untersucht und der Anteil der mit Rübenfliegen befallenen Pflanzen ermittelt. Der Datensatz ist im agridat package unter dem Namen jansen.carrot hinterlegt.

library(agridat)
data(jansen.carrot)
library(data.table)
# kürzerer Name und data.table Format
carrot <- as.data.table(jansen.carrot) 
head(carrot) # erste 6 Zeilen
##    trt gen block  n  y
## 1:  T1 G01    B1 53 44
## 2:  T1 G01    B2 48 42
## 3:  T1 G01    B3 51 27
## 4:  T1 G02    B1 48 24
## 5:  T1 G02    B2 42 35
## 6:  T1 G02    B3 52 45
str(carrot)
## Classes 'data.table' and 'data.frame':   96 obs. of  5 variables:
##  $ trt  : Factor w/ 2 levels "T1","T2"..
##  $ gen  : Factor w/ 16 levels "G01",""..
##  $ block: Factor w/ 3 levels "B1","B2"..
##  $ n    : int  53 48 51 48 42 52 49 49..
##  $ y    : int  44 42 27 24 35 45 8 16 ..
##  - attr(*, ".internal.selfref")=<exter..

Je Parzelle gibt y die Anzahl schadhafter Mohrrüben und n die Gesamtzahl von Mohrrüben an. Der prozentuale Befall kann also je Parzelle als y/n angegeben werden. Wir können eine neue Spalte mit eben diesem Anteil berechnen und genauer anschauen.

# Erstelle Spalte "anteil"
carrot[, anteil := y/n] 
boxplot(anteil ~ gen, data=carrot, las=2)

boxplot(anteil ~ trt, data=carrot, las=2)

Beim Betrachten der Boxplots fällt auf, dass die Genotypen 1 und 2 einen höheren Anteil befallener Möhren aufweisen und somit anfälliger scheinen.

Problemfall Prozentzahlen in einem linearen Modell

Wäre in diesem Experiment ein anderes Merkmal wie z.B. der Ertrag gemessen worden, würde die Auswertung jetzt wie bei einer “normalen” zweifaktoriellen Varianzanalyse in einem randomisierten vollständigen Blockdesign laufen (siehe Kapitel 4). Die Besonderheit an der Auswertung von Prozentzahlen ist, dass oft mindestens eine der folgenden Annahmen einer Varianzanalyse verletzt wird (Die angegebene Reihenfolge der Voraussetzungen entspricht dabei ihrer Wichtigkeit):

  1. Additivität der systematischen Effekte
  2. Varianzhomogenität
  3. Normalverteilung der Fehler

Warum das so ist kann man sich z.B. vorstellen indem man sich vergegenwärtigt, dass Werte einer Sorte weniger schwanken, wenn sie gegen 0% oder 100% gehen. Die Werte einer mittelmäßig anfälligen Sorte hingegen schwanken wahrscheinlich mehr - sowohl über als auch unter 50%. Demnach ist die Varianz bei mittlerem Befall am größten, während sie bei kleinem oder hohem Befall sinkt. Anders ausgedrückt besteht ein direkter Zusammenhang zwischen dem Erwartungswert (einer Sorte) und der Varianz, was der Annahme der Varianzhomogenität widerspricht, weil laut dieser die Varianzen immer gleich sein müssten, egal wie hoch der Erwartungswert (der Sorte).

Um dies zu prüfen können wir ein einfaches, lineares Modell wie in Kapitel 4 aufstellen und uns die Residuenplots anschauen. Tatsächlich deutet der linke Plot ganz und gar nicht auf Varianzhomogenität hin:

# Modell für zweifaktorielle Varianzanalyse in rand. vollst. Blockanlage
mod <- lm(anteil ~ block + gen + trt + gen:trt, data=carrot)
# Residuenplots mit autoplot Funktion aus ggfortify package
library(ggfortify)
autoplot(mod, which=1:2)

Es muss klar sein, das die Varianzanalyse auch dann (in R) durchgeführt werden kann und somit Ergebnisse erzeugt werden, wenn die Annahmen verletzt sind. Diese Ergebnisse sind dann aber nicht zulässig/aussagekräftig. Es liegt daher in der Verantwortung des Auswertenden sicherzustellen, dass die Annahmen erfüllt sind.

Sieht das jetzt gut genug aus?
Häufig wird bemängelt, dass eine Beurteilung von Residuenplots subjektiv ist: Gerade mit wenig Erfahrung fällt die Entscheidung schwer ob der jeweilige Plot nun gut aussieht oder nicht. Um Orientierungspunkte zu geben, sind hier ein paar Residuenplots zum Verlgeich dargestellt: [UNDER CONSTRUCTION].
Es soll auch darauf hingewiesen sein, dass es Tests gibt, die testen ob alle Voraussetzungen erfüllt sind. Allerdings bergen solche Tests das Problem, dass sie das Gesamtbild nur “schwarz-weiß” darstellen, also auf ein “ja” oder “nein” reduzieren, was zu einem Informationsverlust führt. Residuenplots können eben nicht nur perfekt oder vollkommen daneben aussehen und liefern somit eine genauere Aussage darüber ob bzw. wie gut die Voraussetzungen erfüllt sind.

Lösung 1: Datentransformation

Wenn man nicht die gemessenen Werte selbst auswertet, sondern sie vor der Auswertung transformiert, kann dies dazu führen, dass die Annahmen der Varianzanalyse (besser als vorher) erfüllt sind und man doch eine “normale” Varianzanalyse durchführen darf - nur eben mit den transformierten Werten.

Transformieren heißt hier durch eine Funktion umrechnen wie z.B. die Wurzel ziehen oder logarithmieren. Wenn also ein Fall wie der oben auftritt würde man also tatsächlich versuchen nicht den prozentualen Befall selbst (anteil) auszuwerten, sondern eine neue Spalte erstellen in der man beispielsweise log(anteil) berechnet und diese Werte dann auswertet. Auf den ersten Blick mag dies nicht logisch bzw. zielführend erscheinen (zumindest ging es mir so), da man so nicht mehr die natürliche Skala analysiert. Doch tatsächlich können dieselben biologisch relevanten Schlussfolgerungen, also z.B. “Sorte 1 ist anfälliger als Sorte 3” genauso gut auf Basis von logarithmierten Werten wie auf nicht logarithmierten Werten stattfinden - wobei “genauso gut” eben nicht stimmt, wenn die untransformierten Werte die Varianzanalyseannahmen nicht erfüllen, die transformierten Werte aber schon. In solch einem Fall sollte man also die transformierten Werte auswerten.

Ein weiteres Gedankenspiel zur Sinnigkeit von Transformationen
In einem Versuch wurden sowohl Blattlänge als auch Blattfläche bestimmt. Keins der beiden Merkmale kann als wichtiger oder richtiger gewertet werden als das andere. Da Blätter beim Wachsen ihre Form kaum ändern, ergibt sich allerdings kein linearer Zusammenhang zwischen den beiden Merkmalen. Das liegt daran, dass die Fläche etwa der quadrierten Länge entspricht und demnach die Blattfläche schneller größer wird als die Blattlänge. Ein Längenunterschied zwischen zwei Blättern von 1 cm führt bei kleinen Blättern zu einem kleineren Flächenunterschied als bei zwei großen Blättern. Wenn also die Annahmen z.B. für Blattlänge stimmen, dürften sie das automatisch nicht für Blattfläche tun.
Würde man allerdings die Blattflächenwerte transformieren indem man die Wurzel zieht, so ist es naheliegend, dass diese Werte die Annahmen einer Varianzanalyse erfüllen. Die natürliche Skala ist nicht immer die einfachste für eine statistische Analyse und in diesem Fall wäre es von Vorteil die Wurzel aus der Blattfläche anstatt der Blattfläche selbst zu analysieren.

Ziel einer Transformation der Daten ist es also die Voraussetzungen, die für eine Varianzanalyse erforderlich sind, annähernd zu erfüllen. Es ist gängige Praxis für bestimmte Datentypen verschiedene Transformationen anzuwenden wie in der Übersichtstabelle beschrieben (Quelle: Biometrieksript Prof. H.-P. Piepho)


Anwendung

Laut der Tabelle sollen wir also probieren unseren prozentualen Mohrrübenbefall mit der Arcus-Sinus-Wurzel-Transformation zu transformieren und prüfen ob die Annahmen dann erfüllt sind:

# Erstelle Spalte "arcsin.anteil"
carrot[, arcsin.anteil := asin(sqrt(y/n))]
boxplot(arcsin.anteil ~ gen, data=carrot, las=2)

boxplot(arcsin.anteil ~ trt, data=carrot, las=2)

# Ersetze Zielvariable durch "arcsin.anteil"
mod.arcsin <- lm(arcsin.anteil ~ block + gen + trt + gen:trt, data=carrot) 
# Residuenplots
autoplot(mod.arcsin, which=1:2)

Tatsächlich sieht die Varianz homogener als noch im vorigen Residuenplot aus. Hinsichtlich der Annahme der Varianzhomogenität könnte man nun schon eher mit der Varianzanalyse fortfahren und ggf. Mittelwertvergleiche durchführen wie in Kapitel 4. Es soll noch angemerkt sein, dass man am Ende der Auswertung die adjustierten Mittelwerte auch wieder zurücktransformieren kann, sodass sie wieder der natürlichen Skala entsprechen und somit (sowohl in einer Tabelle als auch in Graphen) intuitiver verständlich sind. Es sollte dann aber deutlich vermerkt werden, dass die Varianzanalyse und post-hoc Tests auf einer anderen, nämlich der transformierten Skala, stattgefunden haben.

Oft verbessert eine Transformation die Daten übrigens nicht nur hinsichtlich einer Annahme, sondern aller Annahmen. Das kann man sich so erklären, dass wenn man einen fiktiven Datensatz mit perfekter Einhaltung der Annahmen transformieren würde, er wahrscheinlich nicht nur hinsichtlich einer Annahme Abstriche machen muss.

Für weiterführende, aber dennoch eingängig formulierte Informationen mit direktem Bezug zu R empfehle ich Kapitel 3.3 von Doormann & Kühn.

Lösung 2: Generalisierte Modelle

Hin und wieder müssen Daten analysiert werden, bei denen selbst verschiedene Transformationen zu keiner zufriedenstellenden Verbesserung hinsichtlich der Annahmen einer Varianzanalyse führen. In solch einem Fall können Generalisierte lineare Modelle (GLM) angewendet werden. Bisher haben wir nur einfache, oder auch “generelle lineare Modelle (gLM)” genutzt.

Generalisiert meint hier, dass zwei der Annahmen, nämlich die der Varianzhomogenität und der Normalverteilung der Fehler wegfallen. Es können nämlich andere Verteilungen wie z.B. die Poisson-Verteilung, die Binomial-Verteilung oder die Gammaverteilung angenommen werden. Und ja, die Varianzen müssen nicht länger homogen sein, aber es muss die Beziehung zwischen Erwartungswert und Varianz bekannt sein bzw. geraten werden, d.h. es muss eine Varianzfunktion angegeben werden.

Außerdem bleibt natürlich die übrige Annahme der Additivität der systematischen Effekte bestehen. In die Klasse der GLMs fallen eine ganze Reihe gängiger Modelle, z.B. die Logit- und Probit-Regression für binomialverteilte Daten oder die loglinearen Modelle für poissonverteilte Daten. Die für uns “normalen” Modelle sind also ein Sonderfall der generalisierten Modelle, die normalverteilte Fehler mit homogener Varianz annehmen. Es soll gesagt sein, dass in der Praxis die oben beschrieben Schritte der Datentransformation manchmal auch übersprungen werden und einfach direkt ein GLM für bestimmte Datentypen angepasst wird - so z.B. für Prozentwerte (binomialverteilt) oder nach oben offene Zähldaten (poissonverteilt).

Was ist eigentlich mit Additivität gemeint?
Nehmen wir an in einem Versuch wie dem Mohrrübenversuch hier, gibt es prozentuale Befallswerte von zwei Sorten (A und B) in mehreren Blöcken. In Block 1 hat Sorte A 50% Befall und Sorte B 70% Befall. In Block 2 herrschen deutlich günstigere Infektionsbedingungen, sodass Sorte A einen Befall von 90% aufweist. Das sind 40% mehr als in Block 1. Es ist allerdings unmöglich, dass auch Sorte B 40% mehr Befall in Block 2 aufweist, da der Befall nicht höher als 100% sein kann. Eine Additivität der Effekte ist nicht gegeben, sodass ein additives Modell nicht sinnvoll ist.

Übrigens muss klar sein, dass die Verwendung von GLM nicht automatisch bedeutet, dass keine Datentransformation mehr nötig ist. Stattdessen kann sich mit der eventuellen Datentransformation nun darauf konzentriert werden die letzte verbleibende Annahme der Additivität zu erzielen. Die Logit-Transformation ist in solch einem Fall eine gängige Methode, könnte also auch für das Beispiel im Anmerkungskasten hierüber angewendet werden.

Maximum Likelihood
Einfache lineare Modelle (z.B. Regression oder Varianzanalyse) können mit der Kleinstquadratlösung (auch Methode der minimalen Abweichungsquadrate, ordinary least squares (OLS)) analytisch mittels Matrixalgebra gelöst werden (siehe z.B. Kapitel 6.6 von Doormann & Kühn). Generalisierte lineare Modelle müssen stattdessen mit dem maximum likelihood-Ansatz gelöst werden. Es kann gezeigt werden, dass für normalverteilte Daten bei einfachen linearen Modellen der maximum likelihood-Ansatz zum Kleinsquadratansatz führt. Mittels des maximum likelihood-Ansatzes können demach auch Poisson-, binomial- oder anders verteilte Datensätze analysiert werden. Allerdings ist dieser eben nicht mehr analytisch lösbar. Stattdessen werden Parameter im Modell iterativ numerisch gelöst. Dieser Prozess ist rechenaufwendig und wurde deshalb erst mit der Entwicklung von leistungsfähigen Computern sowie der entsprechenden Statistiksoftware populär.

Anwendung

Um also nun den prozentualen Befall unseres Mohrrübendatensatzes mit einem GLM auszuwerten müssen wir die Funktion glm() nutzen und zusätzlich Argumente zur Verteilung machen.

mod.glm <- glm(y/n ~ block + gen + trt + gen:trt, 
               weights= n,
               family = binomial(link='logit'),
               data   = carrot)

Wie man sehen kann, wäre übrigens das vorige Berechnen der anteil-Spalte für das Aufstellen des Modells nicht nötig gewesen, da wir auch direkt y/n als Zielvariable schreiben können. Wichtig ist dabei aber, dass wir bei Angabe eines solchen Bruchs auch immer noch weights=n dazuschreiben müssen! (Zur Info: Alternativ kann man statt y/n auch cbind(y, n-y) schreiben und das weights=n weglassen)

Als erwartete Verteilung schreiben wir hier vorerst family=binomial, sodass eine Binomialverteilung für die Fehler angenommen wird. Der Zusatz (link="logit") ist hier eigentlich überflüssig, weil dies die default Linkfunktion für binomialverteilte Daten ist - ich schreibe es aber der Klarheit gerne dazu. Da glm() eine Basis-Funktion von R ist - also kein extra package benötigt - können Residuenplots wieder einfach mit der autoplot() Funktion des ggfortify packages betrachtet werden:

autoplot(mod.glm, which=1:2)

Tatsächlich müssen wir die Residuenplots ja aber eben nicht betrachten um Normalverteilung und Varianzhomogenität zu beurteilen, da diese bei einem GLM ja nicht mehr angenommen werden müssen.

Devianzanalyse

Sobald wir mit unserem Modell zufrieden sind, würden wir gerne wie in Kapitel 4 als nächstes mittels F-Test einer zweifaktoriellen Varianzanalyse prüfen wollen welche Behandlungseffekte signifikant sind. Bei GLM werden Varianzanalyse-Tabellen allerdings durch Devianzanalyse-Tabellen ersetzt. Praktischerweise macht dies für den Anwender keinen großen Unterschied in der Codierung und Interpretation der Ergebnisse.

library(car)
anova(mod.glm, test="Chisq") # Option 1
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: y/n
## 
## Terms added sequentially (first to last)
## 
## 
##         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       95     901.70              
## block    2    16.40        93     885.30  0.000274 ***
## gen     15   500.52        78     384.78 < 2.2e-16 ***
## trt      1   107.72        77     277.05 < 2.2e-16 ***
## gen:trt 15    68.91        62     208.14 6.971e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod.glm, test.statistic="Wald") # Option 2
## Analysis of Deviance Table (Type II tests)
## 
## Response: y/n
##         Df   Chisq Pr(>Chisq)    
## block    2  17.065  0.0001969 ***
## gen     15 413.459  < 2.2e-16 ***
## trt      1  94.423  < 2.2e-16 ***
## gen:trt 15  66.875    1.6e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Laut der Devianzanalyse-Tabelle gibt es also eine signifikante Wechselwirkung zwischen den beiden Behandlungseffekten. Demnach müsste man im Anschluss die adjustierten Mittelwerte aller Genotype-Treatment-Kombinationen berechnen und miteinander (via post-hoc Tests) vergleichen wie z.B. in Kapitel 4 oder Kapitel 5.

Über- und Unterdispersion

Wir wollen uns allerdings noch eines weiteren Aspektes der GLM zuwenden - dem Dispersionsparameter. Wie schon gesagt gibt es für jede Verteilung eine zugehörige Varianzfunktion und somit eine erwartete Varianz. Tatsächlich sind die beobachteten Varianzen aber manchmal größer oder kleiner als ihr Erwartungswert laut Varianzfunktion. Dieses Phänomen wird dann entsprechend als Über- und Unterdispersion bezeichnet.

Sollte solch ein Fall vorliegen, kann im Modell ein zusätzlicher Dispersionsparameter \(\phi\) angepasst werden, der die Über- oder Unterdispersion sozusagen auffängt. Läge weder Über- noch Unterdispersion vor, so wäre \(\phi=1\). Bei Unterdispersion ist \(\phi<1\) und bei Überdispersion ist \(\phi>1\).

Als Kontrolle lassen wir uns den Dispersionsparameter des eben angepassten Modells ausgeben:

summary(mod.glm)$dispersion 
## [1] 1

Es muss aber klar sein, dass \(\phi\) hier nicht auf 1 geschätzt wurde, sondern auf 1 fixiert war, da wir explizit nur nach eine Binomialverteilung gefragt hatten. Um einen von 1 abweichenden Dispersionsparameter zu erlauben, müssen wir in glm() das binomial durch quasibinomial ersetzen:

mod.glm2 <- glm(y/n ~ block + gen + trt + gen:trt, 
               weights= n,
               family = quasibinomial(link='logit'),
               data   = carrot)

summary(mod.glm2)$dispersion 
## [1] 3.215226

Und tatsächlich liegt eine relativ starke Überdispersion vor. Da das Modell sich nun verändert hat, ergeben sich auch andere Testergebnisse in der Devianzanalyse-Tabelle:

Anova(mod.glm2, test.statistic="Wald")
## Analysis of Deviance Table (Type II tests)
## 
## Response: y/n
##         Df    Chisq Pr(>Chisq)    
## block    2   5.3077    0.07038 .  
## gen     15 128.5940  < 2.2e-16 ***
## trt      1  29.3673  5.988e-08 ***
## gen:trt 15  20.7994    0.14336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Die Wechselwirkung ist nun nicht mehr signifikant.

 

Bei Fragen kannst du mir gerne schreiben!

schmidtpaul@hotmail.de