In dieser Sitzung befassen wir uns mit dem Thema Kontrastanalyse bei mehrfaktoriellen Versuchdesigns. Wir schauen uns dazu wieder das Depressionsbeispiel an, das Sie schon aus dem Video zur 2x2 ANOVA kennengelernt haben. Die Hypothesen damals waren: (Haupeteffekt 1) Dass Patienten nach einem Hirnschlag im Schnitt zu stärkerer Depressivität neigen als Patienten nach einem Herzinfarkt. (Haupteffekt 2) Dass Patienten, bei denen der Vorfall (Herzinfarkt bzw. Schlaganfall) überdurchschnittlich stark war, zu mehr Depressivität neigen als Patienten, bei denen die Vorfälle durschnittlich stark waren. (Interaktionseffekt) Dass der Unterschied in der Depressivität zwischen über- und durchschnittlichen Ereignissen bei Patienten mit Schlaganfall größer ist als bei Patienten mit Herzinfarkt.

Wie Sie schon wissen erlaubt eine 2x2 ANOVA die Testung von Haupteffekten und Interkationen aber eben keine gerichtete Testung. Die Richtung muss mit Hilfe von Abbildungen beurteilt werden. Eine Kontrastanalyse dagegen erlaubt eine gerichtete Testung. Die einfaktorielle Kontrastanalyse kennen Sie schon aus Quanti I. Heute schauen wir uns an, wie Kontrastanalysen auch für mehrfaktorielle Designs gerechnet werden können.

Wir laden zunächst die Daten und erstellen Abbildungen.

0.0.1 Daten einlesen und Pakete laden

data <- read.csv("Uebung9_Kontraste_Krankheitsstudie.csv")

library(Rcmdr)
library(RcmdrMisc)

# define factors
data$Ereignis <- factor(data$Ereignis, levels = c("Herzinfarkt", "Hirnschlag"))
data$Schwere <- factor(data$Schwere, levels = c("Durchschnitt", "ueberdurchschnitt"))

data$Kombiniert <- factor(data$Kombiniert, levels = c("Herz_durchschnitt", "Herz_ueberdurchschnitt", 
                                                        "Hirn_durchschnitt", "Hirn_ueberdurchschnitt"))

0.0.2 Grafik erstellen

with(data, plotMeans(Depressivitaet, Ereignis, Schwere, 
  error.bars="conf.int", level=0.95, connect=TRUE, legend.pos="farright"))

library(plyr)
ddply(data, .(Ereignis, Schwere), summarize,  Rate1=mean(Depressivitaet))
##      Ereignis           Schwere    Rate1
## 1 Herzinfarkt      Durchschnitt 11.30000
## 2 Herzinfarkt ueberdurchschnitt 14.20000
## 3  Hirnschlag      Durchschnitt 14.16667
## 4  Hirnschlag ueberdurchschnitt 19.80000
with(data, plotMeans(Depressivitaet, Ereignis, 
  error.bars="conf.int", level=0.95, connect=TRUE, legend.pos="farright"))

ddply(data, .(Ereignis), summarize,  Rate1=mean(Depressivitaet))
##      Ereignis    Rate1
## 1 Herzinfarkt 12.75000
## 2  Hirnschlag 16.98333
with(data, plotMeans(Depressivitaet, Schwere, 
  error.bars="conf.int", level=0.95, connect=TRUE, legend.pos="farright"))

ddply(data, .(Schwere), summarize,  Rate1=mean(Depressivitaet))
##             Schwere    Rate1
## 1      Durchschnitt 12.73333
## 2 ueberdurchschnitt 17.00000

0.0.3 Aufstellen der statistischen Hypothesen über Konstraste

Was sind Kontraste? Eine Kontrastkodierung ermöglicht das Testen von gerichteten Hypothesen in ein- und mehrfaktoriellen Designs. Das Prinzip ist dabei dasselbe wie bei einer Dummykodierung (siehe VL 6), nur dass gezielt andere (meist inhaltlich interessantere) Vergleiche (sog. Kontraste) kodiert werden können. Dabei sind verschiedene Vergleiche möglich, allerdings sollte jeder Mittelwert nur einmal in einen Vergleich eingehen. Dann ist gewährleistet, dass die Kontraste unabhängig voneinander (orthogonal) sind. Das vermeidet Fehlerkumulation (jeder Kontrast entspricht nur einem statistischen Test, auch wenn mehrere Gruppen in ihn eingehen), und es wird insgesamt genau dieselbe Varianz aufgeklärt wie bei einer Varianzanalyse - mit dem Vorteil, dass wir flexibler und vor allem gerichtet testen können.

0.0.3.1 Statistische Hypothesen (siehe Formelsammlung):

H1,1: Kontrast “Art des Ereignisses” > 0; H0,1: Kontrast “Art des Ereignisses” <=0

H1,2: Kontrast “Schwere des Ereignisses” > 0; H0,2: Kontrast “Schwere des Ereignisses” <= 0

H1,3: Kontrast “Art x Schwere des Ereignisses” > 0; H0,3: Kontrast “Art x Schwere des Ereignisses” <= 0

Kontraste werden auch mit “D” (groß geschrieben, nicht zu verwechseln mit Cohen’s d). Wir könnten die Kontraste oben also auch D1 (HE Art), D2 (HE Schwere) und D3 (Interaktion Art*Schwere) nennen.

0.0.4 Kontraste definieren per Hand (zur Demonstration)

Kontraste werden definiert, indem man den Gruppen des Designs bestimmte Kontrastkoeffizienten zuweist. Wie wählt man diese? Folgende Prinzipien sind relevant für das Vergeben von Kontrastkoeffizienten für Haupteffekte:

  1. Die Gruppen, die für den Haupteffekt miteinander verglichen werden sollen, erhalten Kontrastkoeffizienten mit unterschiedlichen Vorzeichen

  2. Die Kontrastkoeffizienten müssen in der Summe Null ergeben

  3. Für einige Arten der Kontraste ist es nötig, dass die Summe der Beträge der Koeffizienten 2 ergibt (in einem Design mit 4 Gruppen). (Für einige Arten funktioniert es auch mit anderen Beträgen, aber man macht zumindest nichts falsch, wenn sich an diese Regel hält.)

Welche Kontrastkoeffizienten sollten wir also für die erwarteten Haupteffekte in unserem Design vergeben?

Kontrast “Art des Ereignisses”

  • Die Hypothese ist, dass Hirn-Ereignisse zu mehr Depressivität führen als Herz-Ereignisse. Die Hirn-Gruppen sollten also andere Vorzeichen erhalten als die Herz-Gruppen. Wenn wir den Hirn-Gruppen positive Vorzeichen und den Herz-Gruppen negative Vorzeichen geben, entspricht der Kontrast der Richtung unserer Hypothese (andersherum geht es auch, aber dann müssen wir das später bedenken, wenn wir die Ergebnisse interpretieren).

  • Wenn wir uns an die oben aufgelisteten Prinzipien 2 und 3 halten wollen, sollten wir also folgende Kontrastkoeffizienten für die vier Gruppen vergeben:

    • 0.5 für “Hirn, überdurchschnittlich”
    • 0.5 für “Hirn, durchschnittlich”
    • -0.5 für “Herz, überdurchschnittlich”
    • -0.5 für "Herz, durchschnittlich

Kontrast “Schwere des Ereignisses”

  • Die Hypothese ist, dass schwere Ereignisse mehr Depressivität hervorrufen als durchschnittliche Ereignisse. Analog zum Vorgehen oben ergeben sich also folgende Kontrastkoeffizienten:
    • 0.5 für “Hirn, überdurchschnittlich”
    • 0.5 für “Herz, überdurchschnittlich”
    • -0.5 für “Hirn, durchschnittlich”
    • -0.5 für “Herz, durchschnittlich”

Interaktionskontrast

Wie müssen die Kontrastkoeffizienten vergeben werden, wenn wir die gerichtete Interaktionshypothese testen wollen, dass die Schwere des Ereignisses bei Hirnschlag die Depressivität stärker erhöht als bei einem Herzinfarkt? Am leichtesten lässt sich dies herleiten, wenn wir uns die statistische Hypothese noch einmal als Gleichung über Mittelwerte aufschreiben:

M(Hirn,ü) - M(Hirn,d) > M(Herz,ü) - M(Herz,d), also:

[M(Hirn,ü) - M(Hirn,d)] - [M(Herz,ü) - M(Herz,d)] > 0

Wenn wir die eckigen Klammern auflösen, erhalten wir:

M(Hirn,ü) - M(Hirn,d) - M(Herz,ü) + M(Herz,d) > 0

In diese Gleichung setzen wir einfach die Kontrastkoeffizienten ein (jeder Mittelwert wird mit 0.5 gewichtet) und lesen anschließend die Vorzeichen ab:

0.5 x M(Hirn, ü) - 0.5 x M(Hirn,d) - 0.5 x M(Herz,ü) + 0.5 x M(Herz,d) > 0

Es ergeben sich also folgende Kontrastkoeffizienten für unseren gerichteten Interaktionskontrast:

  • 0.5 für “Hirn, überdurchschnittlich”
  • -0.5 für “Hirn, durchschnittlich”
  • -0.5 für “Herz, überdurchschnittlich”
  • 0.5 für “Herz, durchschnittlich”

Hier zeigt sich schonmal ein deutlicher Unterschied zum unspezifischen Interaktionstest der Varianzanalyse: Wenn wir solche Kontrastkoeffizienten vergeben, können wir gezielt testen, obb dieses Interaktionsmuster vorliegt.

0.0.5 Anlegen der Kontraste im R-Commander

Es gibt (mindestens) zwei Wege, diese Kontrastkoeffizienten in R einzuspeisen. Erstens könnten wir manuell drei Hilfsvariablen anlegen (wie beim manuellen Anlegen der Dummy-Variablen), eine für jeden Kontrast. In diesen Hilfsvariablen würde jede Person den Koeffizienten ihrer Gruppe erhalten (Beispiel: In der Spalte für den Haupteffektkontrast “Art des Ereignisses” würden alle Personen, die in der Hirn-Gruppe waren, eine 0.5 erhalten und alle in der Herzgruppe eine -0.5. Analog würde für die anderen beiden Kontrastspalten vorgegangen). Diese drei Hilfsvariablen würden wir in einem linearem Modell dann als Prädiktoren für unser Kriterium “Depressivität” verwenden - der Output dieses Modells liefert uns die Ergebnisse für unsere Kontraste.(Wer will, kann das ja mal zur Übung machen und prüfen, ob die gleichen Ergebnisse wie weiter unten erhalten werden - das sollte so sein.)

Etwas schneller geht es, wenn man die Kontrastkoeffizienten über den R-commander definiert und anschließend ein lineares Modell berechnet mit einer Variable als Prädiktor, zu der die Kontraste hinzugefügt wurden.

Dies muss eine Variable sein, die alle im Design vorhandenen Faktoren enthält, deshalb wurde in diesem Beispiel zusätzlich die Variable “Kombiniert” angelegt, in der sowohl Art als auch Schwere des Ereignisses kodiert wurden.

Wir benennen die Kontraste geben für jeden die Koeffizienten ein, die wir weiter oben ermittelt haben:

Die Kontraste wurden der Variable “Kombiniert” hinzugefügt. Resultierender code:

.Contrasts <- matrix(c(-0.5,-0.5,0.5,0.5,-0.5,0.5,-0.5,0.5,0.5,-0.5,-0.5,
  0.5), 4, 3)
colnames(.Contrasts) <- c('HE_Art', 'HE_Schwere', 'Interaktion')
contrasts(data$Kombiniert) <- .Contrasts
remove(.Contrasts)

Wir können jetzt ein lineares Modell mit der Variable “Kombiniert” als Prädiktor und der Variable “Depressivität” als Kriterium definieren und uns die Ergebnisse ausgeben lassen:

KontrastModell <- lm(Depressivitaet ~ Kombiniert, data=data)

summary(KontrastModell)
## 
## Call:
## lm(formula = Depressivitaet ~ Kombiniert, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.167 -2.225  0.200  2.700  6.700 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            14.8667     0.3165  46.978  < 2e-16 ***
## KombiniertHE_Art        4.2333     0.6329   6.689 8.34e-10 ***
## KombiniertHE_Schwere    4.2667     0.6329   6.741 6.43e-10 ***
## KombiniertInteraktion   1.3667     0.6329   2.159   0.0329 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.467 on 116 degrees of freedom
## Multiple R-squared:  0.4498, Adjusted R-squared:  0.4356 
## F-statistic: 31.61 on 3 and 116 DF,  p-value: 5.211e-15
Anova(KontrastModell)
## Anova Table (Type II tests)
## 
## Response: Depressivitaet
##            Sum Sq  Df F value    Pr(>F)    
## Kombiniert 1139.8   3  31.614 5.211e-15 ***
## Residuals  1394.1 116                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effektstärke für die einzelnen Effekte

  • Z.B. für den HE Ereignis
# Cohen's d = Kontrast / siB 
# siB = QS_res/ df_res

siB <- sqrt(1394.1/116) 
D_Ereignis <- 4.2333

d <- D_Ereignis/siB
d
## [1] 1.221128

Interpretation des Outputs

Der Anteil der aufgeklärten Varianz und der zugehörige F-Wert haben genau dieselbe Bedeutung wie immer.

Die “estimates” geben jetzt aber nicht mehr, wie bei der Dummy-Kodierung, die Veränderungen im Kriterium an, wenn sich der Prädiktor um eine Ausprägung ändert. Stattdessen sind sie der Wert der definierten Kontraste (Gewichtung der Gruppenmittelwerte mit den Koeffizienten). Zur Veranschaulichung wurden unten die entsprechenden Mittelwerte und anschließend die Werte der Kontraste nachgerechnet:

M_Herz_d <- mean(data$Depressivitaet[1:30])
M_Herz_u <- mean(data$Depressivitaet[31:60])
M_Hirn_d <- mean(data$Depressivitaet[61:90])
M_Hirn_u <- mean(data$Depressivitaet[91:120])

#Haupteffektkontrast Art
0.5*M_Hirn_d + 0.5*M_Hirn_u - 0.5*M_Herz_d - 0.5*M_Herz_u
## [1] 4.233333
#Das ist nichts anderes als das: Wir ziehen den Mittelwert der Herzgruppen vom Mittelwert der Hirngruppen ab:
mean(data$Depressivitaet[61:120]) - mean(data$Depressivitaet[1:60])
## [1] 4.233333
#Haupteffektkontrast Schwere
-0.5*M_Hirn_d + 0.5*M_Hirn_u - 0.5*M_Herz_d + 0.5*M_Herz_u
## [1] 4.266667
#Das ist nichts anderes als das: Wir ziehen den Mittelwert der "durchschnittlich"-Gruppen vom Mittelwert der "überdurchschnittlich"-Gruppen ab (siehe anderer Haupteffekt).

#Interaktionskontrast
-0.5*M_Hirn_d + 0.5*M_Hirn_u + 0.5*M_Herz_d - 0.5*M_Herz_u
## [1] 1.366667
#Das ist nichts anderes als ein Vergleich der jeweiligen Differenzen (siehe weiter oben bei der Herleitung der Kontrastkoeffizienten für die Interaktion). Genau genommen ist es die halbe Differenz der Differenzen (mit anderen Kontrastkoeffizienten könnte man auch direkt die Differenz der Differenzen erhalten - aber so oder so bedeutet ein höherer Wert dieses Estimates, das die Interaktion in die vorhergesagte Richtung geht.)

Je größer der Wert jedes Kontrasts, desto stärker gehen die Daten in die vorhergesagte Richtung (und desto eher wird der Kontrast natürlich auch im t-test signifikant).

Der t-test über die Kontraste wird zweiseitig gerechnet. Da wir die Kontraste aber ja gerichtet aufgestellt haben, dürfen wir die p-Werte noch halbieren.

Das Intercept ist der Gesamtmittelwert der Daten (also der Depressionsscore, von dem man ausgehen kann, wenn man keine Informationen über Art und Schwere des Ereignisses hat).

Achtung: Das manuelle Anlegen der orthogonalen Kontraste funktioniert nur auf diese Weisen nur, solange alle Gruppen im Design gleich groß sind (ansonsten müssen leicht andere Koeffizienten vergeben werden). Ist das nicht der Fall, sollte man auf voreingestellte Kontraste zurückgreifen (siehe nächster Abschnitt), denn R korrigiert das ggf. automatisch.

0.0.6 Kontrastanalyse mittels voreingestellter Kontraste (empfohlen)

Schneller als mit dem manuellen Anlegen geht die Kontrastanalyse mit voreingestellten Kontrasten. Diese bieten verschiedene Möglichkeiten, die Gruppen eines Designs miteinander zu vergleichen. Wir werden zwei Methoden kennenlernen, die verschiedende Arten von orthogonalen Kontraste ausgeben:

  1. “Differenz” (auch: “sum-to-zero”): Vergleicht die Mittelwerte von Faktorstufen mit dem Gesamtmittelwert der Daten.

  2. “Helmert”: Vergleicht den Mittelwert jeder Faktorstufe mit dem Mittelwert der vorhergehenden Faktorstufen. (siehe Vorlesung, Folie 13)

Bei nur zwei Ausprägungen pro Faktor bietet sich die Differenzmethode an, bei mehr als zwei Ausprägungen empfehlen sich Helmert-Kontraste.

Die Reihenfolge der Faktorstufen ist wichtig, wenn man voreingestellte Kontraste nutzt (besonders bei Helmert-Kontrasten). Wir können die Reihenfolge unter Datenmanagement, Variablen, Neuanordnung der Faktorstufen ändern. Beim Nutzen der Differenzmethode bietet es sich an, bei beiden Faktoren (Art und Schwere) die Gruppe mit dem höheren Mittelwert an die erste Stelle zu setzen. So erhalten wir am Ende einen positiven Wert des Kontrasts. (Das ist nicht zwingend nötig, man sollte nur wissen, in welcher Reihenfolge die Faktorstufen aktuell stehen und was das später für die Interpretation der Ergebnisse heißt.)

data$Ereignis <- with(data, factor(Ereignis, 
  levels=c('Hirnschlag','Herzinfarkt')))

data$Schwere <- with(data, factor(Schwere, 
  levels=c('ueberdurchschnitt','Durchschnitt')))

#Prüfen der korrekten Anordnung:
levels(data$Ereignis)
## [1] "Hirnschlag"  "Herzinfarkt"
levels(data$Schwere)
## [1] "ueberdurchschnitt" "Durchschnitt"

Wir können die Methoden im Commander im selben Menü auswählen wie vorher beim manuellen Anlegen der Kontraste. Die benötigten Koeffizienten werden für uns im Hintergrund automatisch angelegt.

Achtung: Die Kontraste müssen bei diesem Vorgehen nicht mehr der “Kombiniert”-Variable zugefügt werden, sondern den separaten Variablen (Art und Schwere). Es muss also zweimal eine Kontrastmethode ausgewählt werden, für jeden Faktor einzeln (und natürlich beide Male dieselbe Methode):

Resultierender Code:

contrasts(data$Ereignis) <- "contr.Sum"
contrasts(data$Schwere) <- "contr.Sum"

Wir können prüfen, welche Kontrastkoeffizienten R im Hintergrund vergeben hat:

contrasts(data$Ereignis)
##             [S.Hirnschlag]
## Hirnschlag               1
## Herzinfarkt             -1
contrasts(data$Schwere)
##                   [S.ueberdurchschnitt]
## ueberdurchschnitt                     1
## Durchschnitt                         -1

Als nächstes legen wir wieder ein lineares Modell an mit den Variablen Ereignis und Schwere als Prädiktoren und Depressivität als Kriterium. Die definierten Kontraste werden automatisch berücksichtigt. Wenn wir Ereignis und Schwere multiplikaktiv verknüpfen, wird auch automatisch ein Kontrast für eine Interaktion berechnet.

KontrastModell_Diff <- lm(Depressivitaet ~ Ereignis*Schwere, data=data)

summary(KontrastModell_Diff)
## 
## Call:
## lm(formula = Depressivitaet ~ Ereignis * Schwere, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.167 -2.225  0.200  2.700  6.700 
## 
## Coefficients:
##                                                     Estimate Std. Error t value
## (Intercept)                                          14.8667     0.3165  46.978
## Ereignis[S.Hirnschlag]                                2.1167     0.3165   6.689
## Schwere[S.ueberdurchschnitt]                          2.1333     0.3165   6.741
## Ereignis[S.Hirnschlag]:Schwere[S.ueberdurchschnitt]   0.6833     0.3165   2.159
##                                                     Pr(>|t|)    
## (Intercept)                                          < 2e-16 ***
## Ereignis[S.Hirnschlag]                              8.34e-10 ***
## Schwere[S.ueberdurchschnitt]                        6.43e-10 ***
## Ereignis[S.Hirnschlag]:Schwere[S.ueberdurchschnitt]   0.0329 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.467 on 116 degrees of freedom
## Multiple R-squared:  0.4498, Adjusted R-squared:  0.4356 
## F-statistic: 31.61 on 3 and 116 DF,  p-value: 5.211e-15

Der untere Teil des Outputs bleibt nach wie vor identisch. Was geben die Estimates jetzt an?

  • Intercept: Gesamtmittelwert
  • Haupteffekte: Abweichung des ersten Faktorlevels vom Gesamtmittelwert (halbiertes estimate im Vergleich zu den manuell angelegten Kontrasten).
  • Interaktion: Ebenfalls halbiertes estimate im Vergleich zu den manuell angelegten Kontrasten. Oben hat der Estimate des Interaktionsmodells die halbe Differenz der Mittelwertsdifferenzen angegeben. Jetzt wurde dieser Wert nocheinmal halbiert. Inhaltlich ist das nicht mehr sehr interessant. Da wir die Richtung der Interaktion auch nicht manuell definiert haben, bietet es sich hier eher an, diese über den Graphen zu prüfen.

Veranschaulichung Estimates bei Differenzmethode

#Hirn vs. Gesamtmittelwert

M_gesamt <- mean(data$Depressivitaet)

M_Hirn <- mean(data$Depressivitaet[data$Ereignis=="Hirnschlag"])    
M_Hirn - M_gesamt
## [1] 2.116667
#Ueberdurchschnittlich vs. Gesamtmittelwert

M_ueberdurchschnittlich <- mean(data$Depressivitaet[data$Schwere=="ueberdurchschnitt"])

M_ueberdurchschnittlich - M_gesamt
## [1] 2.133333

0.0.7 Beurteilung der statistischen Hypothesen:

  • H1,1 angenommen
  • H1,2 angenommen
  • H1,3 angenommen