Zweifaktorielle Varianzanalyse: Depressivität nach Herzinfarkt vs. Hirnschlag
Neele Engelmann & Simon Stephan
31 Mai 2019
1 Szenario
In einer medizinisch-psychologischen Studie soll die Hypothese geprüft werden, dass sich Depressivität nach schweren Erkrankungen in Abhängigkeit davon unterscheidet, ob Patient*innen einen Herzinfarkt oder einen Hirnschlag erlitten haben. Es wird vermutet, dass ein Hirnschlag zu mehr Depressivität führt als ein Herzinfarkt. Weiterhin soll geprüft werden, ob die Schwere der Erkrankung einen Einfluss hat. Hier lautet die Hypothese, dass ein schwerer Herzinfarkt/Hirnschlag zu mehr krankheitsbedingter Depressivität führt als ein durchschnittlich schwerer.
Zuletzt soll noch die Interaktion der beiden Faktoren geprüft werden. Es wird vermutet, dass die beiden Faktoren interagieren, aber ein konkretes Muster für die Interaktion wurde nicht vorhergesagt.
2 Daten einlesen, Pakete laden
data <- read.csv("Daten_Krankheitsstudie.csv")
library(Rcmdr)
library(RcmdrMisc)3 Statistische Hypothesen über lineares Modell formuliert:
Allgemeine Form aus Formelsammlung zur zweifaktoriellen Anova:
knitr::include_graphics("SH.png")
Für dieses Beispiel ergibt sich also:
Haupteffekt “Art des Ereignisses”:
- H1,1: QS(Modell mit Ereignis)/df(Modell mit Ereignis) > QS(Modell ohne Ereignis)/df(Modell ohne Ereignis)
- H0,1: QS(Modell mit Ereignis)/df(Modell mit Ereignis) <= QS(Modell ohne Ereignis)/df(Modell ohne Ereignis)
Haupteffekt “Schwere des Ereignisses”:
- H1,2: QS(Modell mit Schwere)/df(Modell mit Schwere) > QS(Modell ohne Schwere)/df(Modell ohne Schwere)
- H0,2: QS(Modell mit Schwere)/df(Modell mit Schwere) <= QS(Modell ohne Schwere)/df(Modell ohne Schwere)
**Interaktion "Art des Ereignisses*Schwere des Ereignisses"**
- H1,3: QS(Modell mit Schwere x Ereignis)/df(Modell mit Schwere x Ereignis) > QS(Modell ohne Schwere x Ereignis)/df(Modell ohne Schwere x Ereignis)
- H0,3: QS(Modell mit Schwere x Ereignis)/df(Modell mit Schwere x Ereignis) <= QS(Modell ohne Schwere x Ereignis)/df(Modell ohne Schwere x Ereignis)
Anmerkung: Diese Hypothesen sind streng genommen nicht gerichtet (obwohl größer/kleiner-Zeichen drinstehen), denn sie testen nur, ob der jeweilige Faktor einen Einfluss hat. Die Richtung der Effekte lesen wir in diesem Fall aus den Graphen ab. In der nächsten Sitzung (Kontraste) werden wir sehen, wie man auch gerichtete Mittelwertshypothesen in mehrfaktoriellen Design testen kann.
4 Grafik: Mittelwerte und Konfidenzintervalle
Commander: Grafiken, Plot für arithmetische Mittel, beide Prädiktoren auswählen
with(data, plotMeans(Depressivitaet, Ereignis, Schwere, error.bars="conf.int",
level=0.95, connect=TRUE, legend.pos="farright"))
5 Zweifaktorielle Varianzanalyse, Variante 1: klassisches Anova-Vorgehen über Commander
Statistik, Mittelwerte vergleichen, Mehrfaktorielle Varianzanalyse:
AnovaModel.1 <- lm(Depressivitaet ~ Ereignis*Schwere, data=data, contrasts=list(Ereignis ="contr.Sum", Schwere ="contr.Sum"))
Anova(AnovaModel.1,type="II")## Anova Table (Type II tests)
##
## Response: Depressivitaet
## Sum Sq Df F value Pr(>F)
## Ereignis 537.63 1 44.7364 8.344e-10 ***
## Schwere 546.13 1 45.4436 6.434e-10 ***
## Ereignis:Schwere 56.03 1 4.6625 0.03289 *
## Residuals 1394.07 116
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1with(data, (tapply(Depressivitaet, list(Ereignis, Schwere), mean, na.rm=TRUE))) # means## Durchschnitt ueberdurchschnitt
## Herzinfarkt 11.30000 14.2
## Hirnschlag 14.16667 19.8with(data, (tapply(Depressivitaet, list(Ereignis, Schwere), sd, na.rm=TRUE))) # std. deviations## Durchschnitt ueberdurchschnitt
## Herzinfarkt 3.425563 3.478010
## Hirnschlag 3.553370 3.407902xtabs(~ Ereignis + Schwere, data=data) # counts## Schwere
## Ereignis Durchschnitt ueberdurchschnitt
## Herzinfarkt 30 30
## Hirnschlag 30 30Interpretation des Outputs:
Signifikanter Haupteffekt für “Ereignis”: Ein Modell, das zusätzlich “Art des Ereignisses” als Prädiktor enthält, sagt Depressivität besser vorher als ein Modell, das nur “Schwere des Ereignisses” enthält
Signifikanter Haupteffekt für “Schwere”: Ein Modell, das zusätzlich “Schwere des Ereignisses” als Prädiktor enthält, sagt Depressivität besser vorher als ein Modell, das nur “Art des Ereignisses” enthält.
Signifikante Interaktion: Ein Modell, das auch die Interaktion der beiden Faktoren enthält, sagt Depressivität besser vorher als ein Modell, das nur die beiden Haupteffekte enthält.
Reminder:
6 Varianzanalyse, Variante 2: Modellvergleiche per Hand (QS Typ 2)
Aus der Vorlesung (VL 7, Folie 11):
Diese drei Modellvergleiche können wir natürlich auch per Hand berechnen!
Modellvergleich 1:
Haupteffekt für Art des Ereignisses: Modell 1 nur mit “Schwere” als Prädiktor" (ohne “Ereignis”“, ohne Interaktion) gegen Modell 2 mit”Schwere"" und “Ereignis”" (ohne Interaktion)
M1 <- lm(Depressivitaet ~ Schwere, data=data)
M2 <- lm(Depressivitaet ~ Schwere + Ereignis, data=data)
anova(M1,M2) #Achtung: auf Kleinschreibung des anova-Befehls achten!## Analysis of Variance Table
##
## Model 1: Depressivitaet ~ Schwere
## Model 2: Depressivitaet ~ Schwere + Ereignis
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 118 1987.7
## 2 117 1450.1 1 537.63 43.379 1.347e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Signifikanter Haupteffekt für Art des Ereignisses!
Modellvergleich 2:
Haupteffekt für Schwere des Ereignisses: Modell 3 nur mit “Ereignis” als Prädiktor (ohne “Schwere” und ohne Interaktion) gegen Modell 2 mit “Schwere” und “Ereignis” (ohne Interaktion). Letzteres wurde im vorherigen Chunk bereits erstellt.
M3 <- lm(Depressivitaet ~ Ereignis, data=data)
anova(M3,M2)## Analysis of Variance Table
##
## Model 1: Depressivitaet ~ Ereignis
## Model 2: Depressivitaet ~ Schwere + Ereignis
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 118 1996.2
## 2 117 1450.1 1 546.13 44.064 1.044e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Signifikanter Haupteffekt für Schwere des Ereignisses!
Interaktion: Modell 2 (ohne Interaktion) gegen neues Modell 4 (mit Interaktion)
M4 <- lm(Depressivitaet ~ Schwere*Ereignis, data=data)
anova(M2,M4)## Analysis of Variance Table
##
## Model 1: Depressivitaet ~ Schwere + Ereignis
## Model 2: Depressivitaet ~ Schwere * Ereignis
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 117 1450.1
## 2 116 1394.1 1 56.033 4.6625 0.03289 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Signifikante Interaktion zwischen Art und Schwere des Ereignisses! Woher wissen wir, in welche Richtung die Interaktion verläuft? Graphische Inspektion oder Kontrastanalyse (nächste Sitzung).
7 Warum sind die F-Werte für die Haupteffekte leicht unterschiedlich?
Wenn man den Output vom Anova-Befehl (Set von drei Modellvergleichen) mit dem Output der drei manuellen Modellvergleiche vergleicht, fällt auf, dass die Quadratsummen für jeden der drei Prädiktoren (Ereignis, Schwere und Interaktion) identisch sind, die F-Werte für die Haupteffekte allerdings minimal unterschiedlich. Woran liegt das?
Da die ersten beiden manuellen Modellvergleiche die Interaktion noch nicht beinhalten, unterscheiden sich hier Quadratsumme der Residuen sowie die Freiheitsgrade der Residuen leicht (beides zusammen bildet ja den Nenner des F-Bruchs). Die Quadratsumme der Residuen ist etwas höher, da sie noch Varianz beinhaltet, die später von der Interaktion aufgeklärt werden kann. Die Freiheitsgrade der Residuen sind etwas höher (117 statt 116), weil ja ohne die Interaktion ein Prädiktor weniger im Modell ist - das Modell nur mit Schwere und Ereignis (ohne Interaktion) enthält zwei Dummy-Variablen und ein Intercept, daher ergibt sich hier 117 als Freiheitsgrade.
8 Was ist der Unterschied zwischen dem “anova”-Befehl und dem “Anova”-Befehl aus dem Commander?
“anova” führt genau einen spezifizierten Modellvergleich durch (den man manuell eingibt, dafür muss man die zu vergleichenden Modelle vorher erstellt haben). “Anova” führt die oben durchgeführten drei Modellvergleiche vom Typ II (default) als Set durch (basierend auf den Komponenten des linearen Modells, das man in der Klammer des Befehls bzw. durch Klicken im Commander spezifiziert - man muss die drei Modelle dafür nicht vorher manuell anlegen).
9 Effektstärken
EtaSquared für das Anova-Modell:
library(lsr)
etaSquared(AnovaModel.1)## eta.sq eta.sq.part
## Ereignis 0.21217901 0.27832134
## Schwere 0.21553357 0.28148301
## Ereignis:Schwere 0.02211377 0.03864101Reminder partial vs. generalised eta squared (VL 7, Folie 18):
knitr::include_graphics("Etasquared.png")
Partial Eta Squared ist also der Anteil von Varianz, den ein Faktor noch aufklärt, wenn die durch andere Faktoren bereits aufgeklärte Varianz ignoriert wird.
10 Output eines einfachen linearen Modells mit zwei Haupteffekten und Interaktion
Statt einer Serie von geplanten Modellvergleichen können wir uns auch den Output von M4 per summary-Befehl anschauen, also des Modells, in das beide Haupteffekte und die Interaktion aufgenommen wurden:
summary(M4)##
## Call:
## lm(formula = Depressivitaet ~ Schwere * Ereignis, 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) 11.3000 0.6329 17.854
## Schwereueberdurchschnitt 2.9000 0.8951 3.240
## EreignisHirnschlag 2.8667 0.8951 3.203
## Schwereueberdurchschnitt:EreignisHirnschlag 2.7333 1.2658 2.159
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Schwereueberdurchschnitt 0.00156 **
## EreignisHirnschlag 0.00176 **
## Schwereueberdurchschnitt:EreignisHirnschlag 0.03289 *
## ---
## 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-15Hier wird das Modell mit allen drei Dummy-Variablen gegen das Nullmodell verglichen. Wie sind nun die Koeffizienten zu interpretieren (siehe VL, Folie 18)? Bei mehr als einer Dummy-Variable verändert sich die Interpretation!
- Intercept = Wert der AV, wenn alle Prädiktoren = 0, also für eine Person mit durchschnittlich schwerem Herzinfarkt
- Schwere = Der Unterschied zwischen überdurchschnittllichem und durchschnittlichem Verlauf bei Herzinfarkt
- Ereignis = Der Unterschied zwischen Hirnschlag und Herzinfarkt bei durchschnittlichem Verlauf
- Interaktion = Das “Extra”, wenn sowohl Hirnschlag als auch überdurchschnittliche Schwere vorliegen. Der vorhergesagte Mittelwert für diese Gruppe ergäbe sich also auf beta0 + beta1 + beta2 + beta3
Zum Abgleich können wir die durch das Modell vorhergesagten Werte dem Datensatz hinzufügen:
data<- within(data, {
fitted.AnovaModel.1 <- fitted(AnovaModel.1)
})wir sehen:
- Für Herzinfarkt & durchschnittlich wird 11.3 vorhergesagt (Intercept von 11.3)
- Für Herzinfarkt & überdurchschnittlich 14.2 vorhergesagt (= Intercept + beta1 = 11.3 + 2.9)
- Für Hirnschlag & durchschnittlich wird 14.16667 vorhergesagt (= Intercept + beta2 = 11.3 + 2.8667)
- Für Hirnschlag & überdurchschnittlich wird 19.8 vorhergesagt (Intercept + alle betas)