Kapitel 9 Inferenzstatistik
In diesem Kapitel wird die Umsetzung der gängigsten statistischen Modelle in R gezeigt. Dabei liegt der Fokus klar auf der Programmierung und nicht auf dem für die Interpretation notwendigen Statistikwissen. Letzteres sollte aus einem der angegeben Statistikbüchern entnommen werden. Durch diese strikte Trennung können wir in den Unterkapiteln auch fortgeschrittenere Modelle betrachten, die in anderen Programmierbüchern leider häufig ausgelassen werden. Nach diesem Kapitel kannst du (fast) alle deiner Forschungsfragen mit R beantworten.
9.1 Einführung
Welches statistische Modell zur Beantwortung der eigenen Forschungsfragen in Form von Hypothesen geeignet ist, hängt im Wesentlichen von der Art der Hypothese, der zugrunde liegenden Verteilung und dem Skalenniveau der Merkmale ab. Die Ausgabe der Ergebnisse sieht dabei immer ähnlich aus und enthält die wesentlichen Informationen, die in einer wissenschaftlichen Arbeit berichtet werden müssen. Während wir hauptsächlich direkt in R integrierte Funktionen verwenden werden, sind in Einzelfällen zusätzliche Packages notwendig, auf die an geeigneter Stelle hingewiesen wird.
In den folgenden Kapiteln werden wir einen Datensatz namens chemo verwenden, der im remp Package enthalten ist (siehe Kapitel 2.5.2 zur Installationsanleitung). 450 PatientInnen mit einer bösartigen Krebserkrankung haben entweder nur Chemotherapie, nur eine Bestrahlung oder eine Kombination aus beidem erhalten. Neben der Beobachtungszeit (Beob_zeit) unter jeder dieser Behandlungen, wurden ebenfalls weitere Merkmale untersucht, die wir in den folgendenen Kapiteln als abhängige Variablen verwenden werden.
# A tibble: 450 × 8
Stationaer Komorb Schmerzen Lebensqualitaet Infektionen Leukos_t0 Leukos_t6 Beob_zeit
<int> <fct> <int> <int> <int> <dbl> <dbl> <dbl>
1 0 Keine 6 2 4 5.74 5.68 0.833
2 0 Lunge 2 2 3 6.71 6.64 4.80
3 0 Keine 7 2 3 6.03 5.97 1.46
4 1 Keine 2 4 1 5.6 5.55 0.922
# ℹ 446 more rowsDer stationäre Aufenthalt aufgrund einer schweren Infektion Stationaer ist binär und wird durch eine Binomialverteilung modelliert. Die Information über weitere Erkrankungen Komorb (kurz für Komorbiditäten) ist nominal skaliert aus der Multinomialverteilung. Die Variablen Schmerzen und Lebensqualitaet wurden durch einen Fragebogen auf jeweils 10 Stufen erfasst und haben daher ein ordinales Skalenniveau. Die Häufigkeit von Infektionen innerhalb von zwei Jahren hat eine negative Binomialverteilung. Schließlich wurde die Anzahl der Leukozyten zu Beginn der Therapie (Leukos_t0) und 6 Monate später (Leukos_t6) erhoben, die jeweils normalverteilt und intervallskaliert (auch metrisch genannt) sind. Leukozyten sind weiße Blutzellen, die einen integralen Bestandteil des Immunsystems darstellen.
Die beschriebenen interessierenden Größen können möglicherweise Unterschiede bzw. Zusammenhänge im Bezug auf das biologische Geschlecht, die Behandlungsart und dem Alter haben. Diese stellen unsere unabhängigen Variablen in sämtlichen Unterkapiteln dar.
# A tibble: 450 × 3
Alter Geschlecht Behandlung
<dbl> <chr> <fct>
1 49 f Radio
2 46 f Radiochemo
3 41 m Radio
4 43 f Radio
# ℹ 446 more rowsEs handelt sich beim chemo Datensatz nicht um reale Daten. Stattdessen wurden die einzelnen Variablen simuliert, was den Vorteil bietet, dass die wahren Verteilungen und Effekte bekannt sind. Für genaue Informationen über die zur Simulation verwendeten Parameter sei auf die Dokumentation verwiesen.
In R werden für inferenzstatistische Methoden zwei verschiedene Schreibweisen verwendet. Entweder müssen die interessierenden Variablen in Form von Vektoren (oder Wertereihen) aus dem Datensatz mithilfe des Dollar-Operators extrahiert werden (siehe Kapitel 4.5). Oder es wird im Kontext von (verallgemeinerten) linearen Modellen auf eine Formelschreibweise zurückgegriffen, bei der auf der linken Seite einer Tilde ~ die abhängige Variable und auf der rechten Seite eine oder mehrere unabhängige Variablen (Faktoren oder Kovariaten), getrennt von einem Pluszeichen, geschrieben werden. Interaktionsterme können wir mithilfe eines Doppelpunktes hinzufügen. In diesem Buch konzentrieren wir uns auf univariate Fragestellungen mit einer abhängigen Variable und einer oder mehreren unabhängigen Variablen.
Die innerhalb von R ausgegebenen Ergebnislisten sind zwar übersichtlich, allerdings nicht geeignet zum Weiterverarbeiten in Tabellenform beim Exportieren nach Word oder in ein PDF. Dafür gibt es ein dediziertes Package, welches abschließend in Kapitel 9.8 vorgestellt wird.
9.2 Ein- und Zweistichprobenszenarien
9.2.1 t Test und Welch Test
Zum Berechnen der Effektstärken muss das effectsize Package installiert und geladen sein.
Im Einstichproben t Test untersuchen wir ein normalverteiltes, intervallskaliertes Merkmal in Hinblick auf einen unter der Nullhypothese festgelegten Mittelwert (\(\mu\)). In Bezug auf unseren chemo Datensatz fragen wir uns, ob die mittlere Leukoyztenanzahl nach sechs Monaten Behandlung statistisch signifikant von dem Normwert von 6 E9/L abweicht. Wir würden unter der Nullhypothese demnach einen Mittelwert der Spalte Leukos_t6 in Höhe von 6 E9/L erwarten.
In R setzen wir dies mit der Funktion t.test() um, welcher wir beim Untersuchen einer Stichprobe zunächst die interessierende Variable Leukos_t6 als Vektor (bzw. Wertereihe) mithilfe des Dollar-Operators übergeben (siehe Kapitel 4.5). Anschließend müssen wir noch den Mittelwert \(\mu\) festlegen, welchen wir unter der Nullhypothese, dass es keinen Unterschied zwischen beiden Mittelwerten gibt, erwarten würden.
One Sample t-test
data: chemo$Leukos_t6
t = -1.3637, df = 449, p-value = 0.1734
alternative hypothesis: true mean is not equal to 6
95 percent confidence interval:
5.947976 6.009402
sample estimates:
mean of x
5.978689 In der Ausgabe sehen wir den t-Wert (-1.3637), die Freiheitsgrade (449), den p-Wert (0.1734), die Alternativhypothesen, das 95% Konfidenzintervall (5.95; 6.01) und den Mittelwert der Leukozyten (5.98). Bei einem in der Wissenschaft üblichen \(\alpha\)-Niveau (Fehler 1. Art) von 5% ist der auf Basis der geschätzten Teststatistik berechnete p-Wert nicht statistisch signifikant, da 0.1734 größer als 0.05 ist. Die Nullhypothese wird demnach beibehalten. Sprich, die Leukozytenanzahl nach sechs Monaten Behandlung weicht nicht signifikant vom Normwert ab.
Zweiseitige Hypothesentests (Mittelwerte sind gleich oder ungleich) sind in der wissenschaftlichen Praxis zwar Standard, weil die Richtung des Effekt vorher selten bekannt ist. Allerdings können wir in Ausnahmefällen mit dem Argument alternative allerdings auch einseitige Hypothesen testen ("less" oder "greater"). Die Standardeinstellung ist ein zweiseitiger Test, weswegen wir das Argument in unserem Beispiel nicht extra benennen müssen. Das \(\alpha\)-Niveau kann über das Argument conf.level modifiziert werden, welches mit 0.95 (also einem \(\alpha\) von 5%) voreingestellt ist.
Für den t Test nehmen wir neben einer zugrunde liegenden Normalverteilung beider Merkmale außerdem gleiche Varianzen an (siehe Kapitel 9.7). Dies müssen wir explizit mit dem Argument var.equal festlegen, da die Standardeinstellung in R von ungleichen Varianzen ausgeht. Bei ungleichen Varianzen (var.equal = FALSE) wird ein Welch Test berechnet, welcher die Freiheitsgrade entsprechend korrigiert. Die Überprüfungen der Voraussetzungen der Normalverteilung und Varianzhomogenität wird in Kapitel 9.7 erläutert.
Wir sind an dieser Stelle an der Hypothese interessiert, ob sich die mittlere Leukozytenzahl nach sechs Monaten zwischen den biologischen Geschlechtern signifikant voneinander unterscheidet. Um zwei Gruppen in R zu vergleichen, verwenden wir die Formelschreibweise. Auf die linke Seite der Tilde schreiben wir die intervallskalierte Variable Leukos_t6 und auf die rechten Seite die kategorisierende Variable Geschlecht.
Two Sample t-test
data: Leukos_t6 by Geschlecht
t = 1.456, df = 448, p-value = 0.1461
alternative hypothesis: true difference in means between group f and group m is not equal to 0
95 percent confidence interval:
-0.01591524 0.10690584
sample estimates:
mean in group f mean in group m
6.000426 5.954930 Da der p-Wert auch in diesem Beispiel mit 0.1461 größer als 0.05 ist, behalten wir auch hier bei einem \(\alpha\)-Niveau von 5% die Nullhypothese bei. In unserer Stichprobe gibt es folglich keine signifikanten Unterschiede in der Leukoyztenzahl zwischen den Gechlechtern.
Zusätzlich zum t Test sollte eine Effektstärke berechnet werden, die im Kontext von Mittelwertsvergleichen häufig Cohens d oder Hedges g darstellt. Mit der Funktion cohens_d() respektive hedges_g() aus dem effectsize Package sind diese in R implementiert. Zur Berechnung übergeben wir die Variablen fast identisch wie zuvor innerhalb der Funktion t.test(). Es ändert sich lediglich das Argument var.equal zu pooled_sd.
Cohen's d | 95% CI
-------------------------
0.14 | [-0.05, 0.32]
- Estimated using pooled SD.Beachte, dass bei einem t Test nur zwei Gruppen miteinander verglichen werden können. Für mehr als zwei Gruppen kann auf eine Varianzanalyse zurückgegriffen werden (siehe Kapitel 9.3.1).
Für einen Test auf abhängige (verbundene, gepaarte) Stichproben, müssen wir zusätzlich das paired Argument hinzufügen. Unterscheiden sich die Mittelwerte der Leukozytenanzahlen von Behandlungsbeginn im Vergleich zur Messung nach sechs Monaten? Zur Beantwortung dieser Hypothese übergeben wir der Funktion t.test() die beiden Variablen mit dem Dollar-Operator als Vektoren.
Paired t-test
data: chemo$Leukos_t0 and chemo$Leukos_t6
t = 251.83, df = 449, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
0.05926718 0.06019949
sample estimates:
mean difference
0.05973333 Die Leukoyztenanzahl ist zwischen den beiden Messzeitpunkten statistisch signifikant unterschiedlich, da der p-Wert kleiner als 2.2e-16 ist, was die wissenschaftliche Notation für eine Zahl mit 15 Nullen hinter dem Komma darstellt. Die Alternativhypothese wird in diesem Beispiel demnach angenommen. Auch hier können wir die Effektstärke in Form von Cohens d berechnen. Allerdings verwenden wir hier die Funktion repeated_measures_d() aus dem effectsize Package. An dieser Stelle ist der beobachtete Effekt deutlich größer als beim vorherigen t-Test, was wir in Hinblick auf den signifikanten p-Wert auch erwarten würden.
d (rm) | 95% CI
---------------------
0.14 | [0.14, 0.14]
- Adjusted for small sample bias.Bei einem Vergleich von Mittelwerten intervallskalierter (metrischer) Merkmale mit maximal zwei Gruppen können wir t Tests verwenden. Falls die Variablen keine ähnliche Varianz haben, wählen wir stattdessen den Welch-Test. Auch bei Verletzung der Annahme einer Normalverteilung kann der Welch-Test eine gute Wahl sein. Bei gröberen Verletzungen oder bei ordinalem Skalenniveau greifen wir hingegen auf nicht-parametrische Tests wie den Wilcoxon Test zurück, welcher im folgenden Kapitel eingeführt wird.
9.2.2 Wilcoxon Test
Zum Berechnen der Effektstärke muss das effectsize Package installiert und geladen werden.
Der Wilcoxon Test für unabhängige Stichproben ist ein nicht-parametrischer Hypothesentest für intervall- oder ordinalskalierte Variablen, welcher bei Verletzungen der Annahmen von Normalverteilung und Varianzhomogenität Anwendung findet. Dieser wird in der Literatur auch als Mann-Whitney U-Test bezeichnet.
In R ist dieser in Form der Funktion wilcox.test() enthalten, welcher äquivalent zum t Test des vorherigen Kapitels angewandt wird. Beim Vergleich einer intervallskalierten Variable zwischen zwei Gruppen verwenden wir daher auch in diesem Fall die Formelschreibweise mit kategorisierenden Variable auf der rechten Seite der Tilde. Zusätzlich lassen wir uns das 95%iger Konfidenzintervall ausgeben.
Wilcoxon rank sum test with continuity correction
data: Leukos_t6 by Geschlecht
W = 26825, p-value = 0.2571
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-0.02999439 0.10002589
sample estimates:
difference in location
0.03998897 Als Ergebnis erhalten wir die geschätzte Kenngröße W und den zugehörigen p-Wert. Der p-Wert ist mit 0.2571 größer als 0.05 und somit nicht statistisch signifikant. Es wird also die Alternativhypothese abgelehnt und die Nullhypothese beibehalten. Die Leukozytenverteilungen der beiden biologischen Geschlechter unterscheiden sich demnach in unseren Daten nicht nennenswert voneinander. Eine geeignete Effektstärke ist die biseriale Rangkorrelation, die wir mit der Funktion rank_biserial() aus dem effectsize Package berechnen.
r (rank biserial) | 95% CI
---------------------------------
0.06 | [-0.04, 0.17]Im Kontext von abhängigen (bzw. gepaarten) Stichproben findet stattdessen der Wilcoxon-Vorzeichen-Rang Test Anwendung. Dies ist häufig in Szenarien mit wiederholten Messungen für dieselbe Person der Fall. In unserem Datensatz haben wir die Leukozytenzahl zu zwei verschiedenen Zeitpunkten gemessen. Durch das Setzen des Arguments paired auf TRUE erhalten wir den entsprechenden Schätzer. Wie im vorherigen Kapitel übergeben wir die beiden Spalten als Vektoren mithilfe des Dollar-Operators (siehe Kapitel 4.5).
Wilcoxon signed rank test with continuity correction
data: chemo$Leukos_t0 and chemo$Leukos_t6
V = 101475, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
0.05996573 0.06004989
sample estimates:
(pseudo)median
0.05996195 Ähnlich wie beim t Test des vorherigen Kapitels finden wir auch hier einen statistisch signifikanten p-Wert. Der Wert ist in diesem Fall kleiner als \(2.2*10^{-16}\), was einer Zahl kleiner 0.00000000000000022 entspricht. Andere Statistikprogramm wie SPSS geben bereits bei p-Werten kleiner 0.001 keinen exakten Wert mehr aus. Auch hier können wir erneut die biseriale Rangkorrelation berechnen. Die augenscheinlich perfekte Korrelation ist der Natur der simulierten Daten geschuldet und kommt im Kontext echter Daten normalerweise nicht vor.
r (rank biserial) | 95% CI
--------------------------------
1.00 | [1.00, 1.00]Die beiden Varianten des Wilcoxon Tests finden in der Praxis häufig zur Beantwortung der gleichen Fragestellungen wie der t Test Anwendung. Die Wilcoxon Tests werden tendenziell bei kleineren Stichproben, gröberen Verletzungen der Voraussetzungen der Normalverteilung und Varianzhomogenität oder bei ordinalem Skalenniveau präferiert. Falls nur die Annahme der Varianzhomogenität verletzt ist, kann alternativ auch der Welch-Test verwendet werden (siehe Kapitel 9.2.1).
9.3 Unterschiede mehrerer Gruppen
9.3.1 Varianz- und Kovarianzanalyse
Für die Funktionen dieses Kapitels müssen die Packages car und effectsize installiert und geladen werden.
Im Vergleich zum t-Test und Wilcoxon Test der vorherigen Kapitel bietet eine Varianzanalyse drei Vorteile:
- Es können mehr als zwei Gruppen betrachtet werden.
- Mehr als eine unabhängige Variable kann in das Modell eingeschlossen werden.
- Interaktionsterme zwischen den unabhängigen Variablen können geschätzt werden.
Ansonsten setzen wir auch hier eine Normalverteilung der intervallskalierten Merkmale und Varianzhomogenität voraus (siehe Kapitel 9.7). Im Kontext der Varianzanalyse testen wir die Nullhypothese, ob alle Mittelwerte gleich sind, gegen die Alternativhypothese mindestens eines statistisch signifikant abweichenden Mittelwertes.
Wir möchten im Folgenden untersuchen, inwiefern sich die Leukozytenanzahl nach sechs Monaten Behandlung zwischen den drei Behandlungsarten und den biologischen Geschlechtern unterscheidet. Die zu überprüfende Nullhypothese des Faktors Behandlung ist demnach die Gleichheit der Mittelwerte der PatientInnen mit Radiochemotherapie, Chemotherapie oder Radiotherapie. Falls auch nur für eine dieser Gruppen ein nennenswert unterschiedlicher Mittelwert in unserer Stichprobe geschätzt wird, könnten wir bereits einen signifikanten p-Wert dieses Haupteffektes in der Varianzanalyse erhalten.
Die Berechnung einer Varianzanalyse besteht in R aus zwei Schritten. Zunächst muss das Modell mit der Funktion aov() aufgestellt werden (Akronym für analysis of variances, engl. für Varianzanalyse). Auf der linken Seite der Tilde steht dabei die abhängige Variable und auf der rechten Seite in diesem Fall unsere beiden unabhängigen Variablen, getrennt von einem Pluszeichen. Interaktionseffekte könnte man in das Modell mit einem zusätzlichen Term mit Doppelpunkt hinzufügen. Eine Interaktion (auch Moderation genannt) zwischen den Variablen Behandlung und Geschlecht würde man demnach mit Behandlung:Geschlecht erhalten. Im Rahmen von Kovarianzanalysen wird dieses Konzept später im Kapitel praktisch angewendet.
Das Ausgeben der Ergebnisse erfolgt mit der Funktion Anova() aus dem car Package. Andere Statistikprogramme wie SPSS geben standardmäßig die Typ 3 Quadratsummen aus, weswegen auch wir uns mit dem Argument type dafür entscheiden. Alternativ könnten über dieses Argument ebenfalls die Typ 2 Quadratsummen berechnet werden. Typ 1 Quadratsummen erhält man, indem man das Ergebnis von aov() der Funktion summary() übergibt (z.B. summary(result)).
result <- aov(Leukos_t6 ~ Behandlung + Geschlecht, data = chemo)
res_aov <- Anova(result, type = 3)
res_aovAnova Table (Type III tests)
Response: Leukos_t6
Sum Sq Df F value Pr(>F)
(Intercept) 4225.0 1 39004.9324 < 2e-16 ***
Behandlung 0.8 2 3.7144 0.02513 *
Geschlecht 0.2 1 2.1527 0.14302
Residuals 48.3 446
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In der Ergebnisliste erhalten wir Informationen über die Quadratsummen (Sum Sq), Freheitsgrade (df) und die Teststatistik in Form des F-Wertes mit entsprechenden p-Wert. Daraus können wir schließen, dass unsere Daten bei einem üblichen \(\alpha\)-Niveau von 5% keinen nennenswerten Unterschied zwischen den biologischen Geschlechtern in Hinblick auf die Leukozytenzahl nach sechs Monaten Therapie zeigen (da 0.14302 größer als 0.05). Bei der Behandlungsart sieht es anders aus. Weil wir hier einen p-Wert von 0.02513 geschätzt haben und dieser kleiner als 0.05 ist, weicht mindestens eine mittleren Leukozytenzahlen der Gruppen statistisch signifikant von den anderen ab.
Ob der Unterschied nun zwischen Radiochemotherapie und Chemotherapie, Radiochemotherapie und Radiotherapie oder Chemotherapie und Radiotherapie besteht, können wir aus diesem Ergebnis noch nicht ableiten. Für die genaue Untersuchung des Haupteffekts Behandlungsgruppe verwenden wir an dieser Stelle einen sogenannten Post-hoc Test nach Tukey, welcher innerhalb von R durch die Funktion TukeyHSD() implementiert ist (Akronym für Honest Significant Differences).
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Leukos_t6 ~ Behandlung + Geschlecht, data = chemo)
$Behandlung
diff lwr upr p adj
Chemo-Radiochemo -0.03151523 -0.1207342 0.05770374 0.6841799
Radio-Radiochemo -0.10123872 -0.1906078 -0.01186968 0.0217816
Radio-Chemo -0.06972349 -0.1592404 0.01979339 0.1604867Ausgeben wird eine Tabelle mit der Mittelwertsdifferenz aller Gruppenunterschiede (diff) mit 95% Konfidenzintervall (lwr, upr) und korrigierten p-Wert. Der Haupteffekt ist in der Varianzanalyse folglich aufgrund des Unterschiedes in der Leukoyztenzahl zwischen Radiochemotherapie und Radiotherapie entstanden. Die Leukozyten sind in der Gruppe mit reiner Bestrahlung (Radiotherapie) geringer als bei der Kombinationstherapie mit einem statistisch signifikanten p-Wert (da 0.022 kleiner 0.05).
Eine Alternative dazu stellen paarweise t Tests dar, welche jeweils die Mittelwerte der Behandlungsgruppen vergleichen und für jede Kombination einen p-Wert zurückgibt. Da wir an dieser Stelle mehrfach Hypothesentests anwenden und es somit zur Inflation des Fehlers 1. Art kommt, wird standardmäßig auch hier der p-Wert korrigiert. Über das p.adjust.methodArgument können unter anderem Bonferroni ("bonferroni"), Benjamini-Hochberg ("BH") oder Holm (“holm”) verwendet werden.
Pairwise comparisons using t tests with pooled SD
data: chemo$Leukos_t6 and chemo$Behandlung
Radiochemo Chemo
Chemo 0.407 -
Radio 0.024 0.102
P value adjustment method: BH Auch hier können wir mit einem p-Wert von 0.024 einen statistisch signifikanten Unterschied zwischen der mittleren Leukozytenzahl bei PatientInnen mit Radiochemotherapie und Radiotherapie erkennen.
Die übliche Effektstärke im Kontext von Varianzanalysen ist das partielle \(\eta^2\), welches Aufschluss über den korrigierten Anteil der Varianzaufklärung der jeweiligen unabhängigen Variable gibt. In R berechnen wir diese mit der Funktion eta_squared() aus dem effectsize Package und setzen das Argument partial auf TRUE.
# Effect Size for ANOVA (Type III)
Parameter | Eta2 (partial) | 95% CI
------------------------------------------
Behandlung | 0.02 | [0.00, 1.00]
Geschlecht | 4.80e-03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Wir erhalten das partielle \(\eta^2\) für beide Haupteffekte. Beachte an dieser Stelle, dass 4.80e-03 die wissenschaftliche Notation für 0.0048 ist.
Die sogenannte Kovarianzanalyse (kurz ANCOVA) korrigiert den Haupteffekt um den Einfluss einer intervallskalierten (metrischen) Kovariate. So könnten wir möglicherweise die Vermutung haben, dass der Einfluss der Behandlungsmethode auf die Leukozytenzahl maßgeblich vom Alter der PatientInnen abhängt. Zum Beispiel könnten die weißen Blutzellen bei älteren Menschen generell niedriger sein. Zur Berücksichtigung der Kovariate fügen wir sie mit einem Pluszeichen ins Modell hinzu.
Anova Table (Type III tests)
Response: Leukos_t6
Sum Sq Df F value Pr(>F)
(Intercept) 130.083 1 1208.3795 < 2e-16 ***
Behandlung 0.770 2 3.5743 0.02884 *
Alter 0.531 1 4.9323 0.02686 *
Residuals 48.012 446
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sowohl die Behandlungsart als auch das Alter haben einen statistisch signifikanten Einfluss auf die Leukoyztenzahl nach 6 Monaten Behandlung. Interaktionseffekt können wir entweder mit einem Doppelpunkt oder als Kurzform durch einen Asterisk (*) schätzen.
aov(Leukos_t6 ~ Behandlung + Geschlecht + Behandlung:Geschlecht, data = chemo)
aov(Leukos_t6 ~ Behandlung * Geschlecht, data = chemo)Die abhängige Variable muss bei der Varianzanalyse immer intervallskaliert (metrisch) und die unabhängigen Variablen nominal skaliert sein. In der Kovarianzanalyse können wir die Haupteffekte um die Einfluss einer intervallskalierten Kovariate in Form eines Interaktionsterms korrigieren. Die Varianzanalyse ist gegenüber Verletzungen der Voraussetzungen der Normalverteilung und Varianzhomogenität relativ robust. Erst bei gröberen Verletzungen müssen wir auf nicht-parametrische Tests wie dem Kruskal-Wallis Test zurückgreifen (siehe Kapitel 9.3.3).
9.3.2 Varianzanalyse mit Messwiederholung
In diesem Unterkapitel müssen die Packages afex und emmeans installiert und geladen sein.
Die Varianzanalyse für abhängige Stichproben wird häufig auch als Varianzanalyse mit Messwiederholung bezeichnet (engl. repeated-measures ANOVA). In diesem Falle haben wir, neben den bereits im vorherigen Kapitel untersuchten Zwischensubjektfaktoren (Behandlung und Geschlecht), ebenfalls Innersubjektfaktoren vorliegen. Jede Person wird also mehrfach untersucht. In unserem Beispiel wurde die Leukozytenanzahl bei Therapiebeginn und sechs Monate später erfasst. In diesem Kapitel untersuchen wir, ob sich die mittlere Leukozytenzahl zwischen diesen Zeitpunkten in Abhängigkeit der anderen unabhängigen Variablen unterscheidet.
Bevor wir zur Berechnung kommen, müssen wir den Datensatz zunächst in ein langes Format bringen (siehe Kapitel 6.6).
chemo_long <- chemo |>
select(Leukos_t0, Leukos_t6, Geschlecht, Behandlung, Pat_id, Alter) |>
pivot_longer(
cols = c(Leukos_t0, Leukos_t6),
values_to = "Leukos",
names_to = "Zeitpunkt"
)
chemo_long# A tibble: 900 × 6
Geschlecht Behandlung Pat_id Alter Zeitpunkt Leukos
<chr> <fct> <int> <dbl> <chr> <dbl>
1 f Radio 1 49 Leukos_t0 5.74
2 f Radio 1 49 Leukos_t6 5.68
3 f Radiochemo 2 46 Leukos_t0 6.71
4 f Radiochemo 2 46 Leukos_t6 6.64
# ℹ 896 more rowsAnschließend übergeben wir der Funktion aov_4() aus dem afex Package, wie gewohnt, auf der linken Seite der Tilde die abhängige Variable und auf der rechten Seite der Tilde die unabhängigen Variablen. Innerhalb runder Klammern werden die Innersubjektfaktoren (hier Zeitpunkt) definiert, gegeben (|) des Personenidentifikators (hier Pat_id). Die Ziffer 4 befindet sich im Funktionsnamen, weil die verwendete Formelschreibweise auf dem lme4 Package basiert, was uns an dieser Stelle nicht weiter interessieren muss. Zusätzlich können wir mit dem Argument anova_table() innerhalb einer Liste diverse Anpassungen vornehmen. Hier verändern wir die Effektstärke zum partiellen \(\eta^2\). Standardmäßig wird ansonsten das generalisierte \(\eta^2\) berechnet, welches die explizite Definition der beobachteten Variablen mithilfe des observed Arguments benötigt. Also alle Variablen, die wir nicht experimentell verändert haben (in unserem Fall das Geschlecht).
result3 <- aov_4(
formula = Leukos ~ Behandlung + Geschlecht + (Zeitpunkt | Pat_id),
data = chemo_long,
anova_table = list(es = "pes")
)
result3Anova Table (Type 3 tests)
Response: Leukos
Effect df MSE F pes p.value
1 Behandlung 2, 446 0.22 3.70 * .016 .026
2 Geschlecht 1, 446 0.22 2.18 .005 .141
3 Zeitpunkt 1, 446 0.00 63467.03 *** .993 <.001
4 Behandlung:Zeitpunkt 2, 446 0.00 0.43 .002 .651
5 Geschlecht:Zeitpunkt 1, 446 0.00 3.72 + .008 .054
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Die Ergebnisse der Zwischensubjektfaktoren Behandlung und Geschlecht sind in etwa wie im vorherigen Kapitel bei der Varianzanalyse für unabhängige Stichproben. Der Innersubjektfaktor Zeitpunkt zeigt einen statistisch signifikanten mit p-Wert kleiner als 001 und somit kleiner als 0.05. Insbesondere auch die Interaktion zwischen Behandlung und Zeitpunkt interessant, welche die Veränderung über die Zeit zwischen den Behandlungsarten modelliert. Ein nicht statistisch signifikanter p-Wert von 0.651 weist daraufhin, dass die Unterschiede zwischen den Behandlungen in Bezug auf die Leukozytenzahl keinen Effekt über die Zeit aufweist.
Die Werte ohne Rundung würden wir mit summary(result3) erhalten. Bei mehr als zwei Stufen des Innersubjektfaktors gibt, würde summary(result3) die Ergebnisse ohne Korrektur, mit Greenhouse-Geisser Korrektur und Hyunh-Feldt Korrektur sowie das Ergebnis des Mauchly Tests für Sphärizität ausgeben. Da wir hier nur zwei Messzeitpunkte vorliegen haben, erhalten wir nur die Schätzwerte der Varianzanalyse.
Als Post-hoc Test schätzen wir die Randmittel mithilfe der Funktion emmeans() aus dem gleichnamigen Package (Akronym für estimated marginal means, engl. für geschätzte Randmittelwerte). Neben dem Ergebnis der Varianzanalyse müssen wir der Funktion zusätzlich die zu untersuchende Gruppe nachfolgend an eine Tilde übergeben.
Behandlung emmean SE df lower.CL upper.CL
Radiochemo 6.05 0.0269 446 6.00 6.10
Chemo 6.02 0.0270 446 5.97 6.07
Radio 5.95 0.0271 446 5.90 6.00
Results are averaged over the levels of: Geschlecht, Zeitpunkt
Confidence level used: 0.95 Um die unterschiedlichen Mittelwerte durch Hypothesentests miteinander zu vergleichen, verwenden wir an dieser Stelle die Funktion pairs(). Das Ergebnis entspricht dabei dem des vorherigen Kapitel, da Behandlung als Zwischensubjektfaktor auch im Kontext einer Varianzanalyse für unabhängige Stichproben adäquat analysiert werden kann.
contrast estimate SE df t.ratio p.value
Radiochemo - Chemo 0.0315 0.0381 446 0.827 0.6868
Radiochemo - Radio 0.1015 0.0382 446 2.659 0.0221
Chemo - Radio 0.0700 0.0382 446 1.831 0.1607
Results are averaged over the levels of: Geschlecht, Zeitpunkt
P value adjustment: tukey method for comparing a family of 3 estimates Die Funktion emmeans() erlaubt es auch, die Behandlungsarten nach Zeitpunkt oder Geschlecht weiter aufzutrennen. Dabei müssen die einzelnen Faktoren mit einem Pluszeichen voneinander getrennt werden. Auch hier könnte man Hypothesentests mithilfe der pairs() Funktion durchführen.
Behandlung Zeitpunkt emmean SE df lower.CL upper.CL
Radiochemo Leukos_t0 6.08 0.0270 446 6.03 6.13
Chemo Leukos_t0 6.05 0.0271 446 6.00 6.10
Radio Leukos_t0 5.98 0.0272 446 5.93 6.03
Radiochemo Leukos_t6 6.02 0.0268 446 5.97 6.07
Chemo Leukos_t6 5.99 0.0269 446 5.94 6.04
Radio Leukos_t6 5.92 0.0270 446 5.87 5.97
Results are averaged over the levels of: Geschlecht
Confidence level used: 0.95 Alle unabhängigen Variablen werden von der Funktion aov_4() standardmäßig als nominale Merkmale angenommen und daher als Faktor umgewandelt, falls diese einen anderen Datentyp haben. Im Falle einer Kovarianzanalyse mit Messwiederholung würde somit ebenfalls die intervallskalierte Kovariate als Faktor formatiert, was keinen Sinn ergibt und in einem Fehler endet. Bei einer Kovarianzanalyse müssen folglich Faktoren im Voraus erstellt werden. In unserem Beispiel ist die Spalte Behandlung zwar bereits ein Faktor, zur besseren Übersicht wird die Umwandlung an dieser Stelle explizit gezeigt (siehe Kapitel 6.4 und 6.10). Die Kovariate (hier Alter) sollte darüber hinaus standardisiert werden.
chemo_long1 <- chemo_long |>
mutate(
Behandlung = factor(Behandlung, levels = c("Radiochemo", "Chemo", "Radio")),
Alter = as.numeric(scale(Alter))
)In der Funktion aov_4() setzen wir nun zusätzlich das Argument factorize auf FALSE, um die automatische Umwandlung der Kovariate Alter in einen Faktor zu verhindern.
aov_4(
formula = Leukos ~ Behandlung + Alter + (Zeitpunkt | Pat_id),
data = chemo_long1,
anova_table = list(es = "pes"),
factorize = FALSE
)Anova Table (Type 3 tests)
Response: Leukos
Effect df MSE F pes p.value
1 Behandlung 2, 446 0.22 3.56 * .016 .029
2 Alter 1, 446 0.22 4.93 * .011 .027
3 Zeitpunkt 1, 446 0.00 63370.93 *** .993 <.001
4 Behandlung:Zeitpunkt 2, 446 0.00 0.41 .002 .664
5 Alter:Zeitpunkt 1, 446 0.00 1.85 .004 .175
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Neben Normalverteilung und Varianzhomogenität müssen wir bei Varianzanalysen für abhängige Stichproben mit mehr als zwei Stufen des Innersubjektfaktors (z.B. drei Messzeitpunkte) ebenfalls die Sphärizität überprüfen. Das Vorliegen der Sphärizität wird automatisch von aov_4() getestet und mit entsprechenden Korrekturen zurückgegeben. Falls zu grobe Verletzungen dieser drei Annahmen vorliegen, können wir auf nicht-parametrische Verfahren wie den Friedman Test wechseln (siehe Kapitel 9.3.4).
9.3.3 Welch Approximation und Kruskal-Wallis Test
Beide in diesem Unterkapitel vorgestellte statistische Tests sind Alternativen für eine Varianzanalyse bei vermuteter grober Verletzungen der Annahmen über Normalverteilung und Varianzhomogenität, welche allerdings auf eine unabhängige Variable limitiert sind.
Falls die Voraussetzung der Normalverteilung erfüllt ist, aber keine Varianzhomogenität besteht, kann die Welch Approximation der Varianzanalyse verwendet werden. Dieser kann als Generalisierung des in Kapitel 9.2.1 kennengelernten t-Tests mit Welch Korrektur verstanden werden. Die Verallgemeinerung bezieht sich hierbei auf die Anwendung bei mehr als zwei Gruppen (hier z.B. drei Behandlungsarten). In R ist diese Methode innerhalb der Funktion oneway.test() implementiert. Wir erhalten hier ein sehr ähnliches Ergebnis zu der Varianzanalyse aus Kapitel 9.3.1, da in unserem Beispiel des chemo Datensatzes sowohl Normalverteilung als auch Varianzhomogenität vorliegt (siehe Kapitel 9.7).
One-way analysis of means (not assuming equal variances)
data: Leukos_t6 and Behandlung
F = 3.8224, num df = 2.00, denom df = 297.59, p-value = 0.02296Im Falle einer Verletzung der Annahme von Normalverteilung und Varianzhomogenität greifen wir stattdessen auf den Kruskal-Wallis Test zurück, der durch die Funktion kruskal.test() in R enthalten ist. Auch hier kriegen wir ein ähnliches Ergebnis zu der Varianzanalyse ausgegeben.
Kruskal-Wallis rank sum test
data: Leukos_t6 by Behandlung
Kruskal-Wallis chi-squared = 5.9832, df = 2, p-value = 0.05021Beide Verfahren bieten eine Alternative zur klassischen Varianzanalyse bei Verletzung der Annahmen von Normalverteilung und/oder Varianzhomogenität. Beachte hierbei, dass nur eine unabhängige, nominal skalierte Variable in Hinblick eines möglichen Einflusses auf ein intervallskaliertes Merkmal untersucht werden kann.
9.3.4 Friedman Test
Eine nicht-parametrische Alternative zur Varianzanalyse für abhängige Stichproben bei Verletzungen von Annahmen der Normalverteilung, Varianzhomogenität oder Sphärizität stellt der Friedman Test dar. Dieses Verfahren hat die Einschränkung, nur den Innersubjektfaktor einzeln untersuchen zu können.
Da auch hier häufig eine Messwiederholung die interessierende nominale unabhängige Variable darstellt, müssen wir zunächst den Datensatz in ein langes Format bringen (siehe Kapitel 6.6).
chemo_long <- chemo |>
select(Leukos_t0, Leukos_t6, Geschlecht, Behandlung, Pat_id) |>
pivot_longer(
cols = c(Leukos_t0, Leukos_t6),
values_to = "Leukos",
names_to = "Zeitpunkt"
)
chemo_long# A tibble: 900 × 5
Geschlecht Behandlung Pat_id Zeitpunkt Leukos
<chr> <fct> <int> <chr> <dbl>
1 f Radio 1 Leukos_t0 5.74
2 f Radio 1 Leukos_t6 5.68
3 f Radiochemo 2 Leukos_t0 6.71
4 f Radiochemo 2 Leukos_t6 6.64
# ℹ 896 more rowsDie Schreibweise unterscheidet sich nur insofern von der Varianzanalyse mit Messwiederholung aus Kapitel 9.3.2, als dass wir keine Klammern um den Messwiederholungsterm auf der rechten Seite der Tilde schreiben. Getrennt sind Innersubjektfaktor und Personenidentifikator durch einen vertikalen Strich.
Friedman rank sum test
data: Leukos and Zeitpunkt and Pat_id
Friedman chi-squared = 450, df = 1, p-value < 2.2e-16Die Ausgabe zeigt auch hier ein statistisch signifikantes Ergebnis. Die Leukozytenanzahl unterscheidet sich demnach zwischen dem Therapiebeginn und nach sechs Monaten Behandlung nennenswert.
9.4 Korrelationsanalysen
9.4.1 Produkt-Moment Korrelation
Bereits in Kapitel 7.3 haben wir Korrelationen in Form im Kontext der deskriptiven Statistik als Zusammenhangsmaße berechnet. Wir können allerdings auch die Nullhypothese eines nicht vorhandenen Zusammenhangs (Korrelation von 0) gegen die Alternativhypothese eines bestehenden Zusammenhangs prüfen. Die Funktion für diese Hypothesentests heißt in R cor.test().
Für die Produkt-Moment Korrelation nach Pearson können wir dem method Argument “pearson” übergeben. Wir müssen allerdings nicht, da dies die Standardeinstellung der Funktion darstellt. Mit der Produkt-Moment Korrelation vergleichen wir zwei intervallskalierte Merkmale miteinander. An dieser Stelle möchten wir herausfinden, ob es einen Zusammenhang zwischen der Leukoszytenanzahl nach sechs Monaten Behandlung und dem Alter gibt. Die Merkmale übergeben wir dabei als Vektoren (bzw. Wertereihen) mithilfe des Dollar-Operators (siehe Kapitel 4.5).
Pearson's product-moment correlation
data: chemo$Leukos_t6 and chemo$Alter
t = 2.2785, df = 448, p-value = 0.02317
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.01473737 0.19751435
sample estimates:
cor
0.10703 Die Ausgabe enthält die Teststatistik in Form eines t-Wertes mit Freiheitsgraden und entsprechendem statistisch signifikanten p-Wert, welcher mit 0.02317 kleiner als das in der Wissenschaft übliche \(\alpha\)-Niveau von 0.05 ist. Außerdem wird die Korrelation in Höhe von 0.10703 mit 95% Konfidenzintervall ausgegeben. Je älter die PatientInnen sind, desto mehr Leukozyten haben wir nach sechs Monaten in unserer Stichprobe gemessen (positive Korrelation). Vielleicht haben die Behandlungen eine einschneidenere Auswirkung auf das Immunsystem jüngerer PatientInnen. Das ist allerdings reine Mutmaßung, da ein vorhandener Zusammenhang nicht zwingend auch eine Kausalität bedeutet.
Ein Sonderfall der Pearson Korrelation ist die Punkt-biseriale Korrelation, welche den Zusammenhang zwischen einem intervallskaliertes Merkmal und einem nominalen Merkmal mit zwei Ausprägungen quantifiziert. Hierfür müssen wir zunächst die Inhalte der Spalte Geschlecht in Zahlen umwandeln (siehe Kapitel 6.4.3). In diesem Fall steht eine 1 für weiblich und eine 0 für männlich.
Anschließend übergeben wir die beiden Spalten auf dieselbe Art und Weise wie zuvor der Funktion cor.test().
Pearson's product-moment correlation
data: chemo1$Leukos_t6 and chemo1$Geschlecht
t = 1.456, df = 448, p-value = 0.1461
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.02396547 0.16004821
sample estimates:
cor
0.06862502 Als Ergebnis erhalten wir eine Korrelation in Höhe von 0.07, die nicht statistisch signifikant von 0 abweicht, da der p-Wert mit 0.1461 größer als 0.05 ist. Dies stimmt mit den Ergebnissen des t Tests aus Kapitel 9.2.1 überein, indem wir keinen Unterschied zwischen den Geschlechtern in Bezug auf die Leukozytenanzahl finden konnten.
Die Produkt-Moment Korrelation quantifiziert den Zusammenhang zwischen zwei intervallskalierten Merkmalen. Falls ordinales Skalenniveau vorliegt, muss stattdessen auf Rangkorrelation zurückgegriffen werden.
9.4.2 Rangkorrelationen
Bei Vorliegen zweier ordinaler Merkmale können wir nicht-parametrischen Rangkorrelationen verwenden. Eine weit verbreitete Methode dafür ist die Rangkorrelation nach Spearman. Im Vergleich zum vorherigen Kapitel verändert sich lediglich das method Argument. Die Merkmale müssen auch hier in Form von Vektoren (bzw. Wertereihen) mithilfe des Dollar-Operators aus dem Datensatz extrahiert werden (siehe Kapitel 4.5).
Wir möchten an dieser Stelle untersuchen, ob das subjektive Schmerzausmaß der PatientInnen mit der empfundenen Lebensqualität zusammenhängt. Beide Merkmale sind durch einen Fragebogen mit zehn Stufen erhoben worden. Auf der Schmerzskala entspricht eine 1 keinen und eine 10 unerträglichen Schmerzen. Bei der Lebensqualität steht eine 1 für unzumutbare Einschränkungen und eine 10 für uneingeschränktes Glück. Die Nullhypothese ist auch bei der Rangkorrelation ein fehlender Zusammenhang, während die Alternativhypothese eine Korrelation ungleich 0 beschreibt.
Warning in cor.test.default(chemo$Schmerzen, chemo$Lebensqualitaet, method = "spearman"): Cannot
compute exact p-value with ties
Spearman's rank correlation rho
data: chemo$Schmerzen and chemo$Lebensqualitaet
S = 23595117, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.5535956 Die Rangkorrelation nach Spearman wird als \(\rho\) (rho) bezeichnet. In der Ausgabe sehen wir eine statistisch signifikante Abweichung der Korrelation von 0, da ein p-Wert von 2.2e-16 die wissenschaftliche Notation für 0.00000000000000022 darstellt und somit deutlich kleiner als das übliche \(\alpha\)-Niveau von 0.05 ist. Wir finden eine negative Korrelation von -0.55 vor. Je mehr Schmerzen die PatientInnen haben, desto geringer ist ihre Lebensqualität in unserer Stichprobe. Umgekehrt, je höher die Lebensqualität ist, desto niedriger die Schmerzen. Auch hier muss bei kausalen Schlussfolgerungen Vorsicht geboten werden, weil Korrelationskoeffizienten lediglich Zusammenhänge quantifizieren.
Zusätzlich erhalten wir eine Warnmeldung, dass sogenannte Ties vorliegen. Ties beschreiben im Kontext von Rangreihungen mehrere Beobachtungen mit demselben Rang. In diesem Fall ist die Berechnung der Kendall-Korrelation angemessener, wofür wir lediglich das method Argument entsprechend anpassen müssen.
Kendall's rank correlation tau
data: chemo$Schmerzen and chemo$Lebensqualitaet
z = -12.197, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
tau
-0.4351842 Die Höhe des negativen Zusammenhangs ist hier mit -0.44 etwas kleiner, aber ebenfalls statistisch signifikant von 0 abweichend.
9.5 Regressionsmodelle
9.5.1 Lineare Regression
Lineare Regressionmodelle werden zur Untersuchung des Zusammenhangs zwischen einer intervallskalierten (metrischen) abhängigen Variable und einer bis zu mehreren unabhängigen Variablen verwendet. Der Hauptunterschied zur einfachen Korrelation des vorherigen Kapitels besteht in der Möglichkeit, mehr als einen möglichen Einfluss auf das interessierende Merkmal zu überprüfen. Die unabhängigen Variablen müssen dabei intervallskaliert (z.B. Alter) oder nominal mit zwei Ausprägungsgraden (z.B. Geschlecht) sein. Bei mehr als zwei Ausprägungsgraden müssen wir Faktoren erstellen (siehe Kapitel 4.3.2 und 6.10). Dabei sollte man sich immer darüber bewusst sein, welche Faktorstufe die Referenzkategorie darstellt, mit denen man die anderen Faktorstufen vergleicht. Die Referenzkategorie ist immer die erste Faktorstufe.
In R erstellt die Funktion lm() ein lineares Regressionsmodell (Akronym für lineares Modell). Auf der linken Seite der Tilde (~) wird die abhängige Variable und auf der rechten Seite eine oder mehrere unabhängige Variablen der Funktion übergeben. Wir stellen uns an dieser Stelle die Frage eines möglichen Zusammenhangs zwischen der Leukozytenanazahl nach sechs Monaten Behandlung und den unabhängigen Variablen Alter und Geschlecht. Eine Übersicht über die Ergebnisse erhalten wir mit der Funktion summary().
Call:
lm(formula = Leukos_t6 ~ Alter + Geschlecht, data = chemo)
Residuals:
Min 1Q Median 3Q Max
-0.92152 -0.23057 0.01753 0.21903 0.87584
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.614790 0.162217 34.613 <2e-16 ***
Alter 0.008676 0.003618 2.398 0.0169 *
Geschlechtm -0.051062 0.031170 -1.638 0.1021
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3294 on 447 degrees of freedom
Multiple R-squared: 0.01735, Adjusted R-squared: 0.01296
F-statistic: 3.947 on 2 and 447 DF, p-value: 0.01998In der Ausgabe sehen wir zunächst die Verteilung der Residuen. Darunter sind die \(\beta\)-Koeffizienten als estimates für jede Kovariate aufgelistet. Der Intercept (\(\beta_0\)) ist entsprechend darüber angegeben. In den nächsten Spalten ist der Standardfehler und der t-Wert mit entsprechendem p-Wert angegeben. Die Teststatistik als t-Wert basiert in diesem Falle auf dem Wald Test, der die Nullhypothese prüft, ob der \(\beta\)-Koeffizient der jeweiligen unabhängigen Variable gleich 0 ist (kein Zusammenhang). Die Alternativhypothese ist hier, ein \(\beta\)-Koeffizient ungleich 0. Beachte außerdem, dass es sich hierbei um partielle Korrelationen handelt. Der Einfluss der jeweils anderen unabhängigen Variablen auf die abhängige Variable wird dabei herausgerechnet.
Der \(\beta\)-Koeffizient der Kovariate Alter weicht statistisch signifikant von 0 ab, da der p-Wert mit 0.0169 kleiner als das übliche \(\alpha\)-Niveau von 0.05 ist. Da der \(\beta\)-Koeffizient positiv ist, haben wir hier einen positiven Zusammenhang vorliegen: je älter die PatientInnen sind, desto höher ist die Leukozytenanzahl und umgekehrt. Beim Geschlecht finden wir hingegen keinen nennenswerten Einfluss auf die Leukozytenanzahl in unserer Stichprobe (0.1021 ist größer als 0.05). Der \(\beta\)-Koeffizient ist negativ, aber welches der Geschlechter hat eine geringer Leukozytenanzahl? Hinter dem Geschlecht ist ein kleines m angehängt, womit R uns sagt, dass das weibliche Geschlecht innerhalb der Funktion mit 0 und das männliche Geschlecht mit 1 kodiert wird. Wenn die Kovariate Geschlecht einen Wert von 1 hat, reduziert sich die Leukozytenanzahl. Folglich ist die Anzahl weißer Blutkörperchen bei Männern nach sechs Monaten Behandlung tendenziell geringer. Allerdings ist dieser Zusammenhang in unserem Fall nicht statistisch signifikant.
Weiter unten in der Ausgabe finden wir noch die Effektstärke der Regressionsanalyse: die Varianzaufklärung in Form des (korrigierten) \(R^2\). Darunter ist das Ergebnis des globalen F-Tests ausgegeben, welcher die Alternativhypothese prüft, ob mindestens einer der \(\beta\)-Koeffizienten ungleich 0 ist. Da eine der Kovariaten einen nennenswerten Zusammenhang zeigt mit der Leukozytenanzahl zeigt, finden wir auch hier einen statistisch signifikanten p-Wert vor (0.01998).
Die Voraussetzungen der linearen Regression sind Normalverteilung und Varianzhomogenität der Residuen. Diese lassen sich am einfachsten und besten graphisch überprüfen (siehe Kapitel 8.13.2 und 9.7).
Auch im Kontext von Regression können wir mithilfe eines Doppelpunktes Interaktionsterme zum Modell hinzufügen. Hier schauen wir uns zusätzlich zu den bisherigen Kovariaten ebenfalls die Interaktion zwischen Alter und Geschlecht an. Hier untersuchen wir also, ob der Alterseffekt möglicherweise je nach Geschlecht anders ausgeprägt ist.
Call:
lm(formula = Leukos_t6 ~ Alter + Geschlecht + Alter:Geschlecht,
data = chemo)
Residuals:
Min 1Q Median 3Q Max
-0.9323 -0.2266 0.0190 0.2186 0.8763
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.871429 0.229189 25.618 <2e-16 ***
Alter 0.002902 0.005134 0.565 0.5721
Geschlechtm -0.562745 0.324841 -1.732 0.0839 .
Alter:Geschlechtm 0.011431 0.007223 1.582 0.1143
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3288 on 446 degrees of freedom
Multiple R-squared: 0.02284, Adjusted R-squared: 0.01627
F-statistic: 3.475 on 3 and 446 DF, p-value: 0.01606Wir können beobachten, dass der Alterseffekt nach Berücksichtigung einer möglichen Interaktion mit dem biologischen Geschlecht nun nicht mehr statistisch signifikant ist. Auch wenn der Interaktionsterm an sich nicht signifikant ist, könnte dies ein Hinweis auf einen gewissen Zusammenhang zwischen Alter und Geschlecht hindeuten.
Falls die nominale unabhängige Variable mehr als zwei Stufen hat, solltest du vor der Analyse einen Faktor daraus machen (siehe Kapitel 6.10). Definiere dabei explizit, welche Gruppe die Referenz darstellen soll (auch Baseline genannt). Die \(\beta\)-Koeffizienten der jeweils anderen Gruppen beziehen sich dann auf Veränderungen im Vergleich zur Referenzgruppe (Stufe 1 des Faktors). An dieser Stelle untersuchen wir mögliche Zusammenhänge zwischen den Behandlungsarten und der Leukozytenanzahl nach sechs Monaten Therapie. Die Referenzgruppe ist in diesem Fall die Behandlung mit Radiochemotherapie.
Call:
lm(formula = Leukos_t6 ~ Behandlung, data = chemo)
Residuals:
Min 1Q Median 3Q Max
-0.9727 -0.2224 0.0088 0.2387 0.8588
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.02272 0.02682 224.58 < 2e-16 ***
BehandlungChemo -0.03152 0.03799 -0.83 0.40722
BehandlungRadio -0.10124 0.03805 -2.66 0.00808 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3295 on 447 degrees of freedom
Multiple R-squared: 0.01629, Adjusted R-squared: 0.01189
F-statistic: 3.701 on 2 and 447 DF, p-value: 0.02545Während beide \(\beta\)-Koeffizienten negativ sind, ergibt nur der Wald Tests der Gruppe der Radiotherapie einen statistisch signifikanten p-Wert (0.00808 ist kleiner als 0.05). Die Leukoyztenzahl ist folglich sowohl bei alleiniger Chemotherapie als auch alleiniger Radiotherapie kleiner als bei der Kombinationstherapie. Allerdings scheint der Zusammenhang beim Vergleich zwischen Radiochemotherapie und Chemotherapie kleiner ausgeprägt zu sein.
Im linearen Regressionsmodell untersuchen wir den Zusammenhang zwischen einer intervallskalierten abhängigen Variable und einer oder mehrerer intervallskalierter oder nominaler unabhängiger Variable. Falls die abhängige Variable nicht intervallskaliert ist, müssen wir auf ein anderes Regressionsmodell zurückgreifen.
9.5.2 Logistische Regression
Logistische Regression verwenden wir zur Untersuchung einer abhängigen Variablen mit nominalem Skalenniveau mit zwei Ausprägungsgraden (0 und 1), die aus einer Binomialverteilung stammen (Synonym: Binomiales Logit Modell).
In R trägt die Funktion zur Modellierung einer logistischen Regression den Namen glm() (Akronym für generalisiertes lineares Modell). Wir möchten an dieser Stelle die Frage beantworten, ob das biologische Geschlecht und die Infektionshäufigkeit innerhalb von zwei Jahren einen Zusammenhang mit dem stationären Aufenthalt aufweist. Die abhängige Variable des stationären Aufenthalts (Stationaer) beinhaltet die Information, ob jemand im Beobachtungszeitraum einen Krankenhausaufenthalt hatte (1) oder nicht (0).
Das Aufstellen der Regressionsgleichung funktioniert hier äquivalent zu den bisher kennengelernten linearen Modellen. Auf der linken Seite der Tilde schreiben wir die abhängige Variable und auf der rechten Seite die unabhängigen Variablen, getrennt von einem Pluszeichen. Zusätzlich definieren wir über das family Argument die Verteilung als binomial (ohne Anführungszeichen). Beachte, dass wir hier die Funktion glm() anstelle von lm() verwenden. Die Ergebnisse werden von der Funktion summary() zusammengefasst.
res_reg4 <- glm(Stationaer ~ Geschlecht + Infektionen, data = chemo, family = binomial)
summary(res_reg4)
Call:
glm(formula = Stationaer ~ Geschlecht + Infektionen, family = binomial,
data = chemo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.99282 0.32229 -9.286 < 2e-16 ***
Geschlechtm 0.19232 0.33572 0.573 0.56674
Infektionen 0.21365 0.06936 3.080 0.00207 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 269.96 on 449 degrees of freedom
Residual deviance: 260.76 on 447 degrees of freedom
AIC: 266.76
Number of Fisher Scoring iterations: 5Während der \(\beta\)-Koeffizient für das Geschlecht nicht statistisch signifikant von 0 abweicht, erhalten wir für die Infektionshäufigkeit einen nennenswerten Zusammenhang mit der Notwendigkeit eines Krankenhausaufenthalts (0.00207 ist kleiner als 0.05). Der positive \(\beta\)-Koeffizient besagt, dass die Wahrscheinlichkeit eines Krankenhausaufenthalts mit steigender Infektionshäufigkeit ansteigt. Das Odds Ratios erhalten wir durch das Exponieren der \(\beta\)-Koeffizienten, die wir mit dem Dollar-Operator aus der Ergebnisliste herausziehen (siehe Kapitel 4.5).
(Intercept) Geschlechtm Infektionen
0.05014566 1.21206341 1.23818556 Im Unterschied zu den linearen Modellen wird uns an dieser Stelle zusätzlich ein Informationskriterium (AIC) ausgegeben, mit welchem wir mehrere ähnliche Modelle miteinander vergleichen könnten. Für eine genaue Beschreibung der Ausgabe von Regressionsmodellen in R und deren Voraussetzungen sei auf Kapitel 9.5.1 verwiesen.
9.5.3 Multinomiales Logit Modell
Für dieses Unterkapitel muss das Package VGAM installiert und geladen sein. Beachte hier, dass der Name des Packages vollständig in Großbuchstaben geschrieben ist.
Möchten wir eine abhängige Variable mit nominalem Skalenniveau mit mehr als zwei Ausprägungen untersuchen, greifen wir auf ein Logit Modell unter multinomialer Verteilungsannahme zurück. Da wir ein nominales Skalenniveau und somit ungeordnete Kategorien vorliegen haben, müssen wir uns auch hier für eine Referenzkategorie entscheiden. Falls die erste Kategorie als Referenz gewählt wird, ist das Modell in der Literatur auch als sogenanntes Baseline Logit Modell beschrieben.
Wir möchten untersuchen, inwiefern die unabhängigen Variablen Alter und Geschlecht einen Einfluss auf drei mögliche Komorbiditäten der PatientInnen hat. Hierbei handelt es sich um weitere Erkrankungen, die neben der Haupterkrankung auftreten können. In unserem Beispiel haben 271 keine weitere Erkrankung, während bei insgesamt 179 zusätzlich entweder eine Erkrankungen der Lungen, des Herzens oder Gehirns auftritt.
Keine Lunge Herz Gehirn
271 67 63 49 Umgesetzt ist das multinomiale Logit Modell in R in Form der vglm() Funktion aus dem VGAM Package. Wir übergeben der Funktion die abhängige Variable auf der linken Seite der Tilde und die unabhängigen Variablen auf der rechten Seite. Entscheiden ist hier das family Argument. Innerhalb der Helferfunktion multinomial() legen wir die Referenzkategorie fest. Hier entscheiden wir uns für die erste Faktorstufe (keine Komorbidität). Die drei weiteren Faktorstufen Lunge, Herz und Gehirn werden im Vergleich zur Referenz betrachtet. Wir erhalten folglich einen Intercept und zwei \(\beta\)-Koeffizienten pro Kategorie. Welche Faktorstufe die erste ist, muss explizit im Zuge der Erstellung des Faktors definiert werden (siehe Kapitel 6.10). Die Zusammenfassung der Ergebnisse wird mit der Funktion summary() erstellt.
reg_res5 <- vglm(
formula = Komorb ~ Geschlecht + Alter,
data = chemo,
family = multinomial(refLevel = 1)
)
summary(reg_res5)
Call:
vglm(formula = Komorb ~ Geschlecht + Alter, family = multinomial(refLevel = 1),
data = chemo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 0.545334 1.425047 0.383 0.702
(Intercept):2 -0.696183 1.466621 -0.475 0.635
(Intercept):3 0.008255 1.619572 0.005 0.996
Geschlechtm:1 -0.071655 0.274680 -0.261 0.794
Geschlechtm:2 -0.097160 0.281295 -0.345 0.730
Geschlechtm:3 -0.047474 0.312092 -0.152 0.879
Alter:1 -0.042804 0.032045 -1.336 0.182
Alter:2 -0.015986 0.032713 -0.489 0.625
Alter:3 -0.037991 0.036388 -1.044 0.296
Names of linear predictors: log(mu[,2]/mu[,1]), log(mu[,3]/mu[,1]), log(mu[,4]/mu[,1])
Residual deviance: 992.3763 on 1341 degrees of freedom
Log-likelihood: -496.1882 on 1341 degrees of freedom
Number of Fisher scoring iterations: 5
No Hauck-Donner effect found in any of the estimates
Reference group is level 1 of the responseDie Referenzkategorie wurde zuvor als keine Komorbiditäten definiert. Wir erhalten drei Intercepts, drei \(\beta\)-Koeffizienten für Geschlecht und drei für das Alter. Dabei wird die Reihenfolge der Faktorenstufen eingehalten. (Intercept):1, Geschlechtm:1 und Alter:1 beschreiben folglich mögliche Geschlechts- oder Alterseffekte auf die Wahrscheinlichkeit, eine Lungenerkrankungen zusätzlich zu der bereits bestehenden Krebserkrankungen zu haben (im Vergleich zu keiner Komorbidität). Äquivalent steht die Zahl 2 für Herzerkrankungen und 3 für Erkrankungen des Gehirns jeweils im Vergleich zu unserer Baseline (keine Komorbidität). Da keiner der \(\beta\)-Koeffizienten signifikant von 0 abweicht, liegt kein Zusammenhang zwischen den unabhängigen Variablen und dem Auftreten einer Komorbidität vor.
Wenn die abhängige Variable nominal skaliert mit mehr als zwei Ausprägungen ist, können wir das multinomiale Logit Modell verwenden. Falls nur zwei Ausprägungen vorliegen, ist die logistische Regression angemessener. Bei ordinalem Skalenniveau greifen wir auf das kumulative Logit Modell zurück.
9.5.4 Kumulatives Logit Modell
In diesem Kapitel muss das VGAM Package installiert und geladen werden. Der Packagename wird in Großbuchstaben geschrieben.
Bei einem interessierenden Merkmal, welches geordnete Kategorien beinhaltet, haben wir ein ordinales Skalenniveau vorliegen und verwenden daher ein kumulatives Logit Modell. Ein anderer etwas unscharfer Begriff dafür ist Ordinales Regressionsmodell. Grundsätzlich gibt es verschiedene kumulative Logit Modelle wie das Adjacent Category Model oder das Proportional Odds Model. Wir werden uns im weiteren Verlauf auf letzteres konzentrieren.
Wir möchten im weiteren Verlauf einen möglichen Zusammenhang zwischen dem subjektiv empfundenen Schmerz während der Behandlung und den unabhängigen Variablen Alter und Geschlecht untersuchen. Die Schmerzen wurden mithilfe eines Fragebogens erfasst, bei dem 1 für keine Schmerzen und 10 für unerträgliche Schmerzen steht.
1 2 3 4 5 6 7 8 9 10
36 46 37 40 86 62 61 56 17 9 Wie im vorherigen Kapitel greifen wir auch hier auf die Funktion vglm() aus dem VGAM Package zurück. Die abhängige Variable wird auf die linke Seite der Tilde und die unabhängigen, getrennt von einem Pluszeichen, auf die rechte Seite geschrieben. Das family Argument wird beim Proportional Odds Model durch die Helferfunktion propodds() definiert. Die Zusammenfassung der Ergebnisse erhalten wir durch die Funktion summary().
reg_parallel <- vglm(
formula = Schmerzen ~ Alter + Geschlecht,
data = chemo,
family = propodds
)
summary(reg_parallel)
Call:
vglm(formula = Schmerzen ~ Alter + Geschlecht, family = propodds,
data = chemo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 2.90034 0.87618 3.310 0.000932 ***
(Intercept):2 1.95733 0.86692 2.258 0.023958 *
(Intercept):3 1.47788 0.86452 1.709 0.087361 .
(Intercept):4 1.05810 0.86311 1.226 0.220232
(Intercept):5 0.27075 0.86174 0.314 0.753379
(Intercept):6 -0.31889 0.86184 -0.370 0.711376
(Intercept):7 -1.05919 0.86392 -1.226 0.220186
(Intercept):8 -2.34958 0.87822 -2.675 0.007465 **
(Intercept):9 -3.44873 0.91845 -3.755 0.000173 ***
Alter -0.01246 0.01920 -0.649 0.516203
Geschlechtm 0.22298 0.16561 1.346 0.178152
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of linear predictors: 9
Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3]), logitlink(P[Y>=4]),
logitlink(P[Y>=5]), logitlink(P[Y>=6]), logitlink(P[Y>=7]), logitlink(P[Y>=8]), logitlink(P[Y>=9]),
logitlink(P[Y>=10])
Residual deviance: 1957.469 on 4039 degrees of freedom
Log-likelihood: -978.7347 on 4039 degrees of freedom
Number of Fisher scoring iterations: 4
No Hauck-Donner effect found in any of the estimates
Exponentiated coefficients:
Alter Geschlechtm
0.987615 1.249799 Für jeden Logit wird ein Intercept ausgegeben. Wir erhalten über alle Logits einen \(\beta\)-Koeffizienten für Alter und einen für Geschlecht. Beide zeigen in diesem Fall keinen signifikanten Zusammenhang mit den einzelnen Logits.
In der Praxis sollte immer die sogenannte Proportional Odds Assumption überprüft werden. Schließlich wissen wir im vornherein nicht, ob die Annahme der gleichen Effekte der Kovariate über alle Kategorien zutrifft. Hierfür stellen wir ein Partial Proportional Odds Model auf, in dem wir für jede Kategorie einen eigenen Effekt schätzen. Hierfür benötigen wir für das family Argument die cumulative() Funktion, in der wir das Argument auf FALSE setzen. Getrennt von einer Tilde wird das partielle Modell definiert. In diesem Fall exkludieren wir den Intercept (-1) und schätzen die Effekte für beide Kovariaten (Alter und Geschlecht).
reg_notparallel <- vglm(
formula = Schmerzen ~ Alter + Geschlecht,
data = chemo,
family = cumulative(parallel = FALSE ~ -1 + Alter + Geschlecht, reverse = TRUE)
)Anschließend vergleichen wir beide Modelle mithilfe des Likelihood Ratio Tests, welcher durch die Funktion lrtest() ebenfalls im VGAM Package enthalten ist. Die Argumente sind hierbei lediglich beide Modelle.
Likelihood ratio test
Model 1: Schmerzen ~ Alter + Geschlecht
Model 2: Schmerzen ~ Alter + Geschlecht
#Df LogLik Df Chisq Pr(>Chisq)
1 4039 -978.73
2 4031 -978.90 -8 0.334 1Da wir einen p-Wert von 1 erhalten, welcher größer als 0.05 ist, können wir von zutreffender Proportional Odds Annahme ausgehen.
Eine andere häufig verwendete Funktion ist polr() aus dem MASS Package, welche Ergebnisse mit umgekehrten Vorzeichen im Vergleich zu der hier kennengelernten Funktion ausgibt. Dies liegt daran, dass in VGAM die Logits als \(logit[P(Y \leq j)] = \alpha_j + \beta x\) definiert sind und im MASS Package als \(logit[P(Y \leq j)] = \alpha_j - \beta x\) mit \(j =1, \dots ,J −1\).
Falls die geordneten Kategorien des interessierenden Merkmals viele Unterstufen haben, ist möglicherweise das Schätzen eines linearen Modell möglich, welches einfacher zu interpretieren ist. In der wissenschaftlichen Praxis werden vor allem im Kontext von Fragebogenstudien häufig Verfahren für intervallskalierte Merkmale angewandt, obwohl Messungen durch einen Fragebogen im Regelfall nur ordinales Skalenniveau vorweisen.
9.5.5 Poisson Regression
Für dieses Kapitel sind die Installation und das Laden des MASS Packages notwendig.
Von sogenannten Zähldaten spricht man normalerweise im Kontext von Häufigkeiten. Die klassische Verteilung für Häufigkeiten ist die Poisson Verteilung, welche den Mittelwert gleich der Varianz postuliert. Dies ist in der Realität dennoch selten der Fall, weil die Varianz häufig bei höherer Häufigkeit steigt. Wir werden im weiteren Verlauf dieses Kapitels zwei mögliche Lösungen für das Problem kennenlernen. Ein weiterer Anwendungsfall für die Possion Regression ist die Modellierung seltener Ereignisse.
Um den Einfluss der Behandlungsarten auf das Immunsystem zu beurteilen, wurde im chemo Datensatz ebenfalls die Infektionshäufigkeit innerhalb von zwei Jahren erhoben. Da die Behandlungsarten mehr als zwei nominale Ausprägungen haben, müssen wir diese als Faktor übergeben (wurde bereits erstellt) (siehe Kapitel 6.10). Die Referenz ist die Radiochemotherapie. Zusätzlich möchten wir noch etwaige Alterseffekte auf die Infektionshäufigkeit ausschließen. Als family Argument verwenden wir poisson (ohne Anführungszeichen). Mit der Funktion summary() werden die Ergebnisse übersichtlich zusammengefasst.
Call:
glm(formula = Infektionen ~ Behandlung + Alter, family = poisson,
data = chemo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.671640 0.331482 2.026 0.0427 *
BehandlungChemo 0.161160 0.078867 2.043 0.0410 *
BehandlungRadio 0.229135 0.077726 2.948 0.0032 **
Alter 0.000258 0.007307 0.035 0.9718
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 874.75 on 449 degrees of freedom
Residual deviance: 865.62 on 446 degrees of freedom
AIC: 1862.6
Number of Fisher Scoring iterations: 5Im Vergleich zur Radiochemotherapie haben PatientInnen unter alleiniger Chemotherapie oder Radiotherapie mehr Infektionen innerhalb von 2 Jahren. Dies erkennen wir an den positiven \(\beta\)-Koeffizienten, die statistisch signifikant von 0 abweichen, weil die p-Werte mit 0.0410 und 0.0032 kleiner als das übliche \(\alpha\)-Niveau von 0.05 sind. Das Alter scheint hingegen keinen relevanten Zusammenhang mit der Infektionshäufigkeit zu haben. Als Effektstärke berechnen wir in diesem Kontext die sogenannte Incidence Rate Ratio, welche wir durch Exponieren der \(\beta\)-Koeffizienten erhalten. Die \(\beta\)-Koeffizienten extrahieren wir anhand des Dollar-Operators aus dem Modell.
(Intercept) BehandlungChemo BehandlungRadio Alter
1.957445 1.174873 1.257512 1.000258 Der Dispersionsparameter ist hier als 1 festgelegt, da in der Poisson Verteilung der Mittelwert gleich der Varianz postuliert wird. Wie bereits angesprochen, ist die aus den Daten geschätzte Varianz häufig höher als durch die Poisson Verteilung angenommen. Man redet in diesem Kontext auch von einer sogenannten Überdispersion (engl. over-dispersion). Eine Möglichkeit zur Korrektur ist die Verwendung der Quasipoisson Verteilung. Ein Nachteil hierbei ist, dass Quasi-Likelihood Methoden keine Möglichkeit zur Berechnung von Informationskriterien bieten.
reg_res8 <- glm(Infektionen ~ Behandlung + Alter, data = chemo, family = quasipoisson)
summary(reg_res8)
Call:
glm(formula = Infektionen ~ Behandlung + Alter, family = quasipoisson,
data = chemo)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.671640 0.453395 1.481 0.1392
BehandlungChemo 0.161160 0.107873 1.494 0.1359
BehandlungRadio 0.229135 0.106312 2.155 0.0317 *
Alter 0.000258 0.009995 0.026 0.9794
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 1.870828)
Null deviance: 874.75 on 449 degrees of freedom
Residual deviance: 865.62 on 446 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5Im Vergleich zur vorherigen Poisson Regression stellen wir bei den Ergebnissen hier einen entscheidenden Unterschied fest. Nur die Radiotherapie hat im Vergleich zur Radiochemotherapie noch einen nennenswerten positiven Zusammenhang mit der Infektionshäufigkeit, die Chemotherapie hingegen nicht mehr. In den Daten scheint folglich eine Überdispersion vorzulegen, worauf auch der auf 1.87 geschätzte Dispersionsparameter hindeutet.
Eine weitere Alternative ohne die Einschränkungen der Quasipoisson Verteilung bei bestehender Überdispersion ist die sogenannte negative Binomial Verteilung. Diese ist in Form der glm.nb() im MASS Package implementiert. Bei der Verwendung bleibt abgesehen vom Wegfall des family Arguments alles gleich.
Call:
glm.nb(formula = Infektionen ~ Behandlung + Alter, data = chemo,
init.theta = 2.592130972, link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6602195 0.4533734 1.456 0.1453
BehandlungChemo 0.1611161 0.1065635 1.512 0.1306
BehandlungRadio 0.2292369 0.1057676 2.167 0.0302 *
Alter 0.0005127 0.0100032 0.051 0.9591
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(2.5921) family taken to be 1)
Null deviance: 516.19 on 449 degrees of freedom
Residual deviance: 511.28 on 446 degrees of freedom
AIC: 1767.4
Number of Fisher Scoring iterations: 1
Theta: 2.592
Std. Err.: 0.408
2 x log-likelihood: -1757.383 Das Ergebnis ist hier sehr ähnlich des Regressionsmodell unter Annahme einer Quasipoisson Verteilung. Hier wird die Überdispersion mit 2.59 sogar noch höher geschätzt. Für eine genaue Beschreibungen der Ausgabe von Regressionsmodellen in R und deren Voraussetzungen sei auch an dieser Stelle auf Kapitel 9.5.1 verwiesen.
Beim Laden des MASS Packages für die negative Binomial Regression wird select() aus dem tidyverse maskiert, weswegen die Reihenfolge beim Laden der beiden Packages relevant ist. Andernfalls hat die kennengelernte Funktion select() aus dem tidyverse nicht mehr die gewohnte Funktionalität. Lade also zunächst MASS und anschließend das tidyverse. Für mehr Informationen über den dafür verantwortlichen Namespace, schaue dir Kapitel 2.5.3 erneut an. Führe im Zweifel einen Neustart von RStudio durch und lade die beiden Packages nun in der für uns richtigen Reihenfolge.
Die klassische Verteilung für Zähldaten ist die Poisson Verteilung. Man sollte diese aber wegen häufig auftretender Überdispersion immer mit den Ergebnissen der Regression unter Annahme einer Quasipoisson oder negativen Binomial Verteilung vergleichen. Falls Häufigkeit nicht als Spalte sondern in Form von Kontingenztafeln vorliegen, muss man stattdessen auf den Exakten Fisher Test oder den \(\chi^2\) Test zurückgreifen (siehe Kapitel 9.6).
9.5.6 Cox Regression
In diesem Kapitel muss das survival Package installiert und geladen werden.
Wir haben in Kapitel 8.13.1 bereits gelernt, wie Überlebenszeiten in Form von Kaplan-Meier Kurven graphisch dargestellt werden können. In diesem Kapitel beschäftigen wir uns mit der statistischen Modellierung durch die Cox Regression (engl. Cox Proportional-Hazards Model) und den Log Rank Test. Uns interessiert die Frage, ob die drei Behandlungsmethoden und das Geschlecht einen Einfluss auf die Überlebenszeit haben.
Im chemo Datensatz beinhaltet die Spalte Beob_zeit die Beobachtungszeit innerhalb eines Zeitraums von 15 Jahren und Status die Information über das Versterben oder sonstige Ausscheiden aus dieser Studie. Eine 0 steht in dieser Spalte für ein sogenanntes zensiertes Ereignis. Der Patient oder die Patientin ist also aus irgendeinem Grund aus der Studie ausgeschieden (z.B. durch Umzug in ein anderes Land). Mit 1 ist der jeweilige Todesfall dokumentiert. Grundsätzlich wäre es außerdem möglich zwei Ereignisse miteinander zu vergleichen, wobei das zweite Ereignis dann entsprechend als 2 kodiert werden muss. Es ist wichtig, dass die Spalte mit den Informationen über die Ereignisse immer genau diese Kodierung hat. Funktionen zur Umkodierung haben wir bereits in Kapitel 6.4.3 kennengelernt.
# A tibble: 450 × 4
Beob_zeit Status Behandlung Geschlecht
<dbl> <dbl> <fct> <chr>
1 0.833 1 Radio f
2 4.80 1 Radiochemo f
3 1.46 1 Radio m
4 0.922 1 Radio f
# ℹ 446 more rowsIm Kontext von Überlebenszeitanalysen interessiert uns häufig zunächst das Überleben nach einer bestimmten Zeit. Hierfür übergeben wir die abhängigen Variablen (Beob_zeit und Status) der Helferfunktion Surv() auf der linken Seite der Tilde und eine 1 auf der rechten Seite. Schließlich sind wir zunächst am gesamten Überleben ohne Unterteilung interessiert.
Call: survfit(formula = Surv(time = Beob_zeit, event = Status) ~ 1,
data = chemo)
n events median 0.95LCL 0.95UCL
[1,] 450 407 1.99 1.82 2.21Das mediane Überleben liegt bei dieser Erkrankung bei knapp 2 Jahren. Für die Überlebenswahrscheinlichkeit nach 5 Jahren übergeben wir die Variable zeit an die Funktion summary(). Durch das times Argument kann die interessierende Überlebenszeit festgelegt werden (hier 5 Jahre).
Call: survfit(formula = Surv(time = Beob_zeit, event = Status) ~ 1,
data = chemo)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
5 47 361 0.129 0.017 0.0993 0.167Nur knapp 13% überleben die ersten fünf Jahre nach Diagnosestellung. Als nächsten Schritt beantworten wir die Frage nach möglichem Einfluss der unabhängigen Variablen Behandlungsart und Geschlecht, welche wir ebenfalls auf die rechte Seite der Tilde schreiben. Die Cox Regression ist mit der Funktion coxph() aus dem survival Package in R verfügbar. Falls eine Modellierung in Abhängigkeit der Start- und Endzeit gewünscht ist, kann die Startzeit dem time und die Endzeit dem time2 Argument übergeben werden. Die Zusammenfassung der Ergebnisse erfolgt durch summary().
reg_res10 <- coxph(
formula = Surv(time = Beob_zeit, event = Status) ~ Behandlung + Geschlecht,
data = chemo
)
summary(reg_res10)Call:
coxph(formula = Surv(time = Beob_zeit, event = Status) ~ Behandlung +
Geschlecht, data = chemo)
n= 450, number of events= 407
coef exp(coef) se(coef) z Pr(>|z|)
BehandlungChemo 1.1703 3.2230 0.1364 8.577 <2e-16 ***
BehandlungRadio 2.2387 9.3816 0.1558 14.372 <2e-16 ***
Geschlechtm 0.8596 2.3622 0.1041 8.260 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
BehandlungChemo 3.223 0.3103 2.467 4.211
BehandlungRadio 9.382 0.1066 6.913 12.731
Geschlechtm 2.362 0.4233 1.926 2.897
Concordance= 0.732 (se = 0.012 )
Likelihood ratio test= 261.1 on 3 df, p=<2e-16
Wald test = 241.1 on 3 df, p=<2e-16
Score (logrank) test = 275.7 on 3 df, p=<2e-16Oben in der Ausgabe erfahren wir, dass nach 15 Jahren nur 43 der 450 Behandelten, die bis zum Ende Teil der Studie waren, noch Leben. Anschließend werden die \(\beta\)-Koeffizienten (coef), Hazards Ratios (exp(coef)), Standardfehler (se), z-Wert und p-Wert.
Die Kovariate Behandlung ist bereits von vornherein als Faktor definiert worden, da hier drei Gruppen vorliegen. Dabei ist die Referenzgruppe die Radiochemotherapie. Wir erhalten also die Informationen, dass unter der Behandlung mit alleiniger Chemotherapie dreimal mehr PatientInnen sterben, mit alleiniger Radiotherapie sogar neunmal mehr. In der Kovariate Geschlecht sehen wir am Ende des Wortes den Buchstaben m. Die Referenzkategorie ist hier folglich das weibliche Geschlecht (f). Männer haben in unserer Stichprobe mehr als die doppelte Gefahr, an der Erkrankung zu sterben. Alle drei \(\beta\)-Koeffizienten weichen statistisch signifikant von 0 ab, da die p-Werte mit 2.2e-16 (\(2.2*10^{-16}\)) deutlich kleiner als das übliche Signifikanzniveau von 5% ist.
Anschließend werden die Hazard Ratios diesmal mit 95% Konfidenzintervall gezeigt. Zuletzt werden die Ergebnisse der drei üblichsten Likelihood Tests für die globale Nullhypothese, dass alle \(\beta\)-Koeffizienten gleich 0 sind, angegeben. Die Alternativhypothese stellt hier das Abweichen mindestens eines \(\beta\)-Koeffizienten von 0 dar.
Eine zentrale Voraussetzung steckt bereits im Namen des Cox Proportional Hazards Modells. Man geht davon aus, dass sich die Hazard Ratio über den zeitlichen Verlauf zwischen den Gruppen proportional zueinander verhält. Es sollte also nicht sein, dass das Risiko in einer Gruppe zu Beginn der Zeitmessung deutlich höher ist und sich dieser Trend zum Ende der Zeitmessung umdreht. Überprüft wird diese Annahme mit dem Score Test für die Schoenfeld Residuen. Im survival Package ist dieser in Form der Funktion cox.zph() umgesetzt.
chisq df p
Behandlung 1.503 2 0.47
Geschlecht 0.846 1 0.36
GLOBAL 2.395 3 0.49Da der p-Wert hier nicht statistisch signifikant ist, behalten wir die Nullhypothese bei, dass die Annahme der Proportional Hazards nicht verletzt ist. Graphisch würde man Verletzungen dieser Annahme an einem Kreuzen der Kaplan-Meier Kurven erkennen.
Ergänzend sei an dieser Stelle noch der Log-Rank Test vorgestellt, der ohne Berücksichtigung von Kovariaten zwei oder mehr Gruppen in Hinblick Unterschiede zwischen den Kaplan-Meier Kurven vergleicht. So könnten wir beispielsweise die Kaplan-Meier Kurven der beiden erhobenen Geschlechter miteinander vergleichen. Der p-Wert des Log-Rank Tests wird wahlweise direkt in der Abbildung eingeblendet (siehe Kapitel 8.13.1). Generell kann man in Abbildung 9.1 schon mit bloßem Auge eine Assoziation zwischen dem Überleben und dem biologischen Geschlecht vermuten, da die Prognose für Männer schlechter zu sein scheint.
Abbildung 9.1: Kaplan-Meier-Kurve unterteilt nach biologischem Geschlecht.
Den dazugehörigen Log-Rank Test berechnen wir mit der Funktion survdiff() aus dem survival Package. Ansonsten ändert sich nichts im Bezug auf die Schreibweise im Vergleich zur zuvor kennengelenrten Cox Regression.
Call:
survdiff(formula = Surv(time = Beob_zeit, event = Status) ~ Geschlecht,
data = chemo)
N Observed Expected (O-E)^2/E (O-E)^2/V
Geschlecht=f 235 205 265 13.5 40.5
Geschlecht=m 215 202 142 25.3 40.5
Chisq= 40.5 on 1 degrees of freedom, p= 2e-10 Auch hier finden wir einen statistisch signifikanten Zusammenhang zwischen dem Geschlecht und der Überlebenszeit. Zusätzlich werden die beobachteten und erwarteten Häufigkeiten ausgegeben, die zur Berechnung miteinander verglichen wurden.
9.6 Kontingenztafeln
9.6.1 Exakter Fisher Test
Zur Berechnung der Effektstärken muss das effectsize Package installiert und geladen werden.
Die Erstellung von Kontingenztafeln wurde bereits in Kapitel 7.2 eingeführt (Synonym: Kreuztabellen). Durch Erstellung dieser Kontingenztafeln können wir die Unabhängigkeit zweier kategoriale Merkmale systematisch überprüfen. Wir möchten wissen, ob ein stationärer Aufenthalt aufgrund irgendeiner Erkrankung abhängig vom biologischen Geschlecht ist. Zunächst müssen wir dafür die Kontingenztafel erstellen.
0 1
f 216 19
m 194 21Im Kontext von \(2 \times 2\) Tabellen und relativ kleinen Stichproben verwenden wir den exakten Test nach Fisher. Es gibt eine umstrittene Daumenregel, die besagt, dass dieser Test nur bei einer Stichprobengröße kleiner als 1000 mit mindestens einer Zelle mit einer Häufigkeit von weniger als fünf verwendet wird. Ansonsten wird der im folgenden Kapitel erklärte \(\chi^2\) Test nach Pearson empfohlen. Heutige Computer können den exakten Fisher Test zwar auch problemlos bei größeren Stichproben berechnen, allerdings wird im Regelfall schon bei moderater Stichprobengröße eine \(\chi^2\) Verteilung approximiert.
Die Funktion innerhalb von R heißt fisher.test(), welcher wir lediglich die Kontingenztafel übergeben müssen. Ein großer Vorteil des exakten Fisher Tests ist die Möglichkeit des gerichteten Hypothesentestens. Anstelle der Abhängigkeit beider Merkmale könnten wir uns also auch fragen, ob der stationäre Aufenthalt häufiger bei Männern oder weniger häufig ist. Dafür müssten wir das Argument alternative von "two.sided" auf "less" oder "greater" ändern.
Fisher's Exact Test for Count Data
data: tbl1
p-value = 0.6195
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.6091001 2.4985032
sample estimates:
odds ratio
1.23003 Als Ausgabe erhalten wir das Odds Ratio mit 95% Konfidenzintervall und einen p-Wert. In unserem Beispiel ist dieser mit 0.6195 größer als das übliche \(\alpha\)-Niveau von 0.05 und somit statistisch nicht signifikant. Es findet sich folglich keine Abhängigkeit zwischen dem Krankenhausaufenthalt und dem biologischen Geschlecht wieder. Effektstärken erhalten wir mit Funktionen aus dem effectsize Package. Im Falle einer \(2 \times 2\) Tafel ist \(\phi\) als Zusammenhangsmaß empfohlen, welches durch die Funktion phi() implementiert ist.
Phi | 95% CI
-------------------
0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Der marginale Effekt nahe Null passt zu unserem nicht-signifikanten Hypothesentest. Alternativ könnten wir bspw. auch Pearson’s C mit der Funktion pearsons_c() berechnen. Das relative Risiko erhalten wir durch riskratio() ebenfalls aus dem effectsize Package.
Risk ratio | 95% CI
-------------------------
1.11 | [0.79, 1.56]9.6.2 Chi Quadrat Tests
Für dieses Kapitel muss das effectsize Package installiert und geladen sein.
Bei größeren Stichproben verwenden wir zum Vergleich zweier kategorialer Merkmale den \(\chi^2\) Test nach Pearson. Auch hier wird eine mögliche Abhängigkeit bzw. ein Zusammenhang zwischen den Merkmalen durch die Alternativhypothese überprüft. Der Vorteil des \(\chi^2\) Tests ist die Anwendung auch auf größeren Tafeln als \(2 \times 2\), der Nachteil ist die fehlende Möglichkeit gerichtete Hypothesen zu überprüfen.
Der \(\chi^2\) Test nach Pearson ist in R in Form der Funktion chisq.test() enthalten. Zunächst möchten wir einen möglichen Zusammenhang zwischen dem stationären Aufenthalt im Krankenhaus und dem biologischen Geschlecht untersuchen. Dafür nehmen wir die im vorherigen Kapitel erstelle Kontingenztafel namens tbl1.
Pearson's Chi-squared test with Yates' continuity correction
data: tbl1
X-squared = 0.21214, df = 1, p-value = 0.6451Das Ergebnis ist in diesem Fall fast identisch mit dem des exakten Fisher Tests (siehe Kapitel 9.6.1). Da die Funktion hier nicht direkt das Odds Ratio mit ausgibt, können wir diese mithilfe der Funktion oddsratio() aus dem effectsize Package berechnen. Als Alternativ wäre an dieser Stelle je nach Studiendesign auch das relative Risiko als Effektstärke denkbar (siehe Kapitel 9.6.1).
Odds ratio | 95% CI
-------------------------
1.23 | [0.64, 2.36]Wir können mit dem \(\chi^2\) Test auch zwei kategoriale Merkmale mit mehr als zwei Ausprägungen untersuchen. Uns interessiert an dieser Stelle, ob ein stationäre Aufenthalt von der Behandlungsart abhängt. Dafür erstellen wir uns zuerst die zugehörige Kontingenztafel (siehe Kapitel 7.2).
0 1
Radiochemo 144 7
Chemo 139 11
Radio 127 22Anschließend verwenden wir, wie zuvor, die Funktion chisq.test() mit der Kreuztabelle als alleiniges Argument.
Pearson's Chi-squared test
data: tbl2
X-squared = 10.174, df = 2, p-value = 0.006178Zurückgegeben wird die Teststatistik in Form des \(\chi^2\)-Wertes, die Freiheitsgrade (df) und der p-Wert. In diesem Fall liegt ein Zusammenhang zwischen dem stationären Aufenthalt und der Behandlungsform vor, da der p-Wert mit 0.006178 kleiner als 0.05 ist. Welche Behandlung häufiger zu einem stationären Aufenthalt führt, vermögen wir aufgrund des \(\chi^2\) Tests nicht zu sagen. Eine geeignete Effektstärke für Tafeln, die größer als \(2 \times 2\) sind, bietet Cramers V (cramers_v) aus dem effectsize Package.
Cramer's V (adj.) | 95% CI
--------------------------------
0.13 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Bei abhängigen (gepaarten) Merkmalen ist der \(\chi^2\) Test nach McNemar die Methode der Wahl. Ein häufiger Anwendungsfall ist daher der Kontext einer Messung zu zwei verschiedenen Zeitpunkten oder mittels zwei verschiedenen Verfahren. Hier schauen wir uns an, ob die Qualität zwischen PCR-Test und Antigen-Schnelltest bei der Diagnostik des SARS-CoV-2 Virus nennenswert unterschiedlich ist. Der zugehörige Datensatz aus dem remp heißt covid.
Schnelltest
PCR Positiv Negativ
Positiv 35 14
Negativ 8 2421Anschließend übergeben wir die Kontingenztafel der Funktion mcnemar.test().
McNemar's Chi-squared test with continuity correction
data: covid
McNemar's chi-squared = 1.1364, df = 1, p-value = 0.2864Da der p-Wert mit 0.2864 in diesem Fall größer als das in der Wissenschaft übliche \(\alpha\)-Niveau von 0.05 ist, behalten wir die Nullhypothese bei. Es besteht also kein nennenswerter Unterschied zwischen beiden diagnostischen Verfahren, auch wenn der Antigen-Schnelltest offensichtlich 8 falsch positive und 14 falsch negative Ergebnisse produziert hat. Bei beiden Verfahren können wir mithilfe des correct Arguments die Korrekturen ausschalten.
Eine adäquate Effektstärke wäre in diesem Fall Cohens g, welche ebenfalls im effectsize Package implementiert ist.
Cohen's g | 95% CI
-------------------------
0.14 | [-0.07, 0.30]Durch die Quadrierung geht bei \(\chi^2\) Tests jede Richtung verloren. Als Konsequenz können ausschließlich zweiseitige Hypothesen getestet werden. Bei kleinen bis moderaten Stichproben in Form von \(2 \times 2\) Tafeln bietet der exakte Fisher Test eine Alternative für gerichtete Hypothesen.
9.7 Voraussetzungen überprüfen
9.7.1 Normalverteilung
Zur graphischen Überprüfung der Verteilungsannahmen muss das Package ggfortify installiert und geladen werden.
Eine der gängigsten statistischen Methoden zur Überprüfung des Vorliegens einer Normalverteilung ist der Shapiro-Wilk Test. Dieser ist in Form der Funktion shapiro.test() direkt in R integriert. Als Argument übergeben wir die entsprechende Variable als Vektor (bzw. Wertereihe) mithilfe des Dollar-Operators (siehe Kapitel 4.5). Die Leukozytenanzahl nach sechs Monaten Behandlung wurde in diversen Kapitels dieses Buches als normalverteilt angenommen. Nun können wir überprüfen, ob diese Annahme zutreffen war.
Shapiro-Wilk normality test
data: chemo$Leukos_t6
W = 0.99747, p-value = 0.7284Der p-Wert ist mit 0.7284 größer als das gängige \(\alpha\)-Niveau von 0.05, weswegen wir die Nullhypothese beibehalten und von normalverteilten Daten ausgehen können. Ein weiterer Anwendungsfall ist die Überprüfung der Normalverteilung der Residuen eines Regressionsmodells. Die Berechnung der Residuen erfolgt durch die in R enthaltene Funktion residuals().
Anschließend übergeben wir die in resi gespeicherten Residuen der Funktion shapiro.test(). Auch hier können wir von einer Normalverteilung ausgehen, weil der p-Wert in Höhe von 0.6904 größer als 0.05 ist.
Shapiro-Wilk normality test
data: resi
W = 0.99734, p-value = 0.6904Eine Alternative der Hypothesentests zur Überprüfung der Normalverteilung stellen graphische Vergleiche der Verteilungen mithilfe eines Quantil-Quantil Plots dar. In Kapitel 8.8 haben wir bereits die Darstellung eines Q-Q Plots für intervallskalierte Variablen gesehen. Auch Residuen wurden bereits mit der Funktion autoplot() aus dem Package ggfortify unter anderem auf Normalverteilung überprüft (siehe Kapitel 8.13.2). In Abbildung 9.2 werden die Quantile der standardisierten Residuen (y-Achse) mit den theoretischen Quantilen der Normalverteilung durch autoplot() verglichen. Solange die Werte auf der Diagonalen liegen, kann von einer Normalverteilung ausgegangen werden.
Abbildung 9.2: Graphische Überprüfungen der Normalverteilung der Residuen
Die meisten kennengelernten Verfahren sind relativ robust gegenüber leichteren Verletzungen der Annahme der Normalverteilung. Hypothesentests wie der Shapiro-Wilk Test sind in solchen Szenarien sehr sensibel gegenüber Verletzungen und verwerfen daher unter Umständen zu häufig die Nullhypothese des Vorliegens einer Normalverteilung. Daher sollte in jedem Fall eine graphische Überprüfung mittels Q-Q Plots erfolgen.
9.7.2 Varianzhomogenität und Homoskedastizität
Für dieses Kapitel muss das car Package installiert und geladen werden.
Viele parametrische statistische Verfahren können nur bei gegebener Varianzhomogenität fehlerfrei geschätzt werden. Zur Überprüfung der Varianzgleichheit zweier Gruppen verwenden wir den F Test, welcher durch die Funktion var.test() in R implementiert ist.
An dieser Stelle möchten wir die Varianz der Leukozytenanzahl nach sechs Monaten Therapie zwischen zwei biologischen Geschlechtern vergleichen. Dafür wird die intervallskalierte Variable auf die linke Seite der Tilde und die nominale auf die rechte Seite geschrieben. Beachte, dass nur zwei Gruppen verglichen werden können. Schließlich werden die beiden Varianzen direkt miteinander ins Verhältnis gesetzt. Ein F-Wert von 1 entspricht demnach einer Varianzhomogenität (Nullhypothese).
F test to compare two variances
data: Leukos_t6 by Geschlecht
F = 0.85053, num df = 234, denom df = 214, p-value = 0.2256
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.6534313 1.1053574
sample estimates:
ratio of variances
0.8505318 Als Ergebnis erhalten wir den F Wert, die Freiheitsgrade, p-Wert und das 95% Konfidenzintervall. Da in diesem Fall die kleiner Varianz der Frauen (0.101) durch die der Männer (0.119) geteilt wird, erhalten wir einen F-Wert kleiner 1. Welche Gruppe im Zähler und welche im Nenner steht, kann explizit durch Umwandlung in einen Faktor kontrolliert werden (siehe Kapitel 6.10). Der F-Wert ist zwar kleiner als 1, allerdings nicht statistisch signifikant, da der p-Wert in Höhe von 0.2256 größer als das häufig gewählte \(\alpha\)-Niveau von 0.05 ist. Die Nullhypothese der Varianzhomogenität wird demnach beibehalten.
Beim Vergleich der Varianz von drei oder mehr Gruppen können wir zwischen dem Levene Test und Bartlett Test wählen. Für den Levene Test verwenden wir die Funktion leveneTest() aus dem car Package. Beachte an dieser Stelle, die abweichende Schreibweise der Funktion mit einem Großbuchstaben beim Wort Test. Wir möchten nun überprüfen, ob die Varianz der Leukozytenanzahl zwischen den Behandlungsarten gleich ist. Dafür übergeben wir auf der linken Seite der Tilde erneut die intervallskalierte Variable und auf der rechten Seite das nominalskalierte Merkmal.
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 0.3454 0.7081
447 Wir erhalten auch hier die Teststatistik eines F-Wertes, die Freiheitsgrade und den p-Wert. Da der p-Wert mit 0.7081 bei eine \(\alpha\)-Niveau von 5% nicht statistisch signifikant ist, behalten wir auch hier die Nullhypothese der Varianzhomogenität bei.
Für die Durchführung des Bartlett Test ändert sich im Vergleich zum Levene Test nur der Funktionsaufruf zu bartlett.test().
Bartlett test of homogeneity of variances
data: Leukos_t6 by Behandlung
Bartlett's K-squared = 1.785, df = 2, p-value = 0.4096Hier kommen wir ebenfalls zur Schlussfolgerung vorliegender Varianzhomogenität, da der p-Wert von 0.4096 ebenfalls nicht statistisch signifikant ist.
Eine Voraussetzung für die Regressionsanalysen ist die Homoskedastizität, welche die Gleichheit der Varianz der Residuen über den gesamten Parameterraum beschreibt. Die beste Überprüfung erfolgt graphisch durch das Auftragen der Residuen gegen die vorhergesagten Werte. Homoskedastizität ist gegeben, solange die Werte gleichmäßig um die horizontale Nulllinie der y-Achse streuen. Die dafür verwendete Funktion autoplot() wurde bereits im Kapitel 8.13.2 vorgestellt.
Abbildung 9.3: Graphische Überprüfungen der Homoskedastizität der Residuen.
Bei Überprüfung der Varianzhomogenität von zwei Gruppen greifen wir auf den klassischen F Test zurück. Bei nominalen Merkmalen mit drei oder mehr Ausprägungen sind sowohl der Levene Test als auch der Bartlett Test Methoden der ersten Wahl. Im Kontext von Regressionsanalysen kann die Varianzhomogenität der Residuen (Homoskedastizität) am besten graphisch überprüft werden (siehe Abbildung 9.3).
9.8 Ergebnisse formatieren
Für das Umformatieren in Tabellenform muss das broom Package installiert und geladen sein.
Als erstes Beispiel verwenden wir das lineare Regressionsmodell aus Kapitel 9.5.1 und nennen die Ergebnisliste model1.
Die wichtigsten Informationen werden mit der Funktion tidy() ausgegeben, welche als einziges Argument die Variable mit dem Modell benötigt. Die \(\beta\)-Koeffizienten sind in der Spalte namens estimate. Die Spalte statistic bezieht sich auf die t-Werte basierend auf den z-Werten der Wald Tests. Die p-Werte werden leider in tibbles immer in wissenschaftlicher Notation dargestellt. Bei der Kovariate Alter ist der p-Wert 1.69e-2, also \(1.69 * 10^{-2}\), und entspricht somit der Zahl 0.0169. Äquivalent dazu ist der p-Wert für die Kovariate Geschlecht 0.102, weil dies eine andere Schreibweise für \(1.02 * 10^{-1}\) darstellt, welches in R als 1.02e-1 angezeigt wird.
# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 5.61 0.162 34.6 1.49e-128
2 Alter 0.00868 0.00362 2.40 1.69e- 2
3 Geschlechtm -0.0511 0.0312 -1.64 1.02e- 1Zusätzliche Informationen wie die Varianzaufklärung (\(R^2\)) oder Informationskriterien in Form des AIC und BIC erhalten wir durch die Funktion glance(). Außerdem ist hier die Teststatistik (F-Wert) des globalen F-Tests als statistic mit zugehörigem p-Wert und Freiheitsgrad (df) ausgegeben. Auch hier muss nur das Modell als Argument übergeben werden.
# A tibble: 1 × 12
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 0.0174 0.0130 0.329 3.95 0.0200 2 -137. 283. 299. 48.5 447
# ℹ 1 more variable: nobs <int>Die Funktion augment() gibt unter anderem die geschätzten (.fitted) Werte, Cooks Distanz (.cooksd) und die Residuen (.resid) zurück.
# A tibble: 450 × 9
Leukos_t6 Alter Geschlecht .fitted .resid .hat .sigma .cooksd .std.resid
<dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 5.68 49 f 6.04 -0.360 0.00676 0.329 0.00273 -1.10
2 6.64 46 f 6.01 0.626 0.00455 0.328 0.00553 1.91
3 5.97 41 m 5.92 0.0505 0.00667 0.330 0.0000530 0.154
4 5.55 43 f 5.99 -0.438 0.00451 0.329 0.00268 -1.33
# ℹ 446 more rowsWir betrachten an dieser Stelle noch ein weiteres Beispiel anhand der Ergebnisse des t-Tests aus Kapitel 9.2.1. Hier entspricht estimate der Mittelwertsdifferenz, estimate1 die mittlere Leukozytenzahl der Frauen und estimate2 die der Männer. In der Spalte statistic wird der t-Wert abgebildet und als parameter die Freiheitsgrade. Leider sind die Bezeichnungen nicht immer eindeutig, weswegen ein Vergleich mit der originalen Ausgabe empfohlen ist.
# A tibble: 1 × 10
estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 0.0455 6.00 5.95 1.46 0.146 448 -0.0159 0.107 Two Sample t-test
# ℹ 1 more variable: alternative <chr>Die Funktion glance() gibt hier dasselbe Ergebnis zurück wie tidy() und augment() kann in diesem Fall nicht verwendet werden.
Während tidy() mit Ausnahme der Varianzanalyse mit Messwiederholung für alle hier vorgestellten statistischen Modellen funktioniert, können glance() und augment() nicht immer angewendet werden. Die aktuelle Liste der unterstützten Funktionen kann in der Dokumentation des broom Packages eingesehen werden.
Beachte, dass die genaue Ausgabe von tidy(), glance() und augment() je nach statistischem Modell unterschiedlich ist. So kann die Spalte estimate bei einem Modell für einen t-Wert stehen, während damit beim Betrachten einer Regression \(\beta\)-Koeffizienten gemeint sind.
Da wir die Ergebnisse nun als tibble in Tabellenform vorliegen haben, kannst du diese im nächsten Schritt in das Textverarbeitungsprogramm deiner Wahl exportieren.