23 Lineare Kontraste
Bei der Analyse eines CRDs haben wir bisher nur untersucht, ob es irgendwo zwischen den verschiedenen Faktorstufen einen relevanten Unterschied gibt. Im nächsten Schritt wollen wir daher bestimmen zwischen welchen Faktorstufen tatsächlich ein Unterschied besteht. Im einfachen Fall aus unserem Laufbeispiel mit drei Faktorstufen (Koffein, Placebo und Control, siehe Abbildung 23.1) können drei verschiedene, sogenannte paarweise Vergleiche durchgeführt werden.
Wir können Koffein mit Placebo, Koffein mit Control und Placebo mit Control vergleichen. In allen drei Fällen vergleichen wir die Mittelwerte \(\mu_i\) der jeweiligen Faktorstufen miteinander. Allerdings haben wir hier das Problem, dass jede Faktorstufe in mehreren (im Beispiel zwei) Vergleichen beteiligt ist. Wenn wir mehr Faktorstufen haben entsprechend noch öfter. Allgemein gilt bei \(K\)-Faktorstufen:
\[\begin{equation} \text{Anzahl der Paarvergleiche} = K(K-1) \end{equation}\]
Dadurch, das jeder Mittelwert \(\mu_i\) in mehreren Vergleichen beteiligt ist, führt dies zu dem Problem der Mehrfachvergleiche.
23.1 Das multiple-comparison Problem
Allgemein gilt, wenn wir mehrere statistische Tests durchführen, und uns den \(\alpha\)-Fehler betrachten, dass über die Gesamtheit der Test der \(\alpha\)-Fehler nicht mehr gleich dem nominalen Level ist. Die Überlegung dazu ist die Folgende. Wir testen \(m\) Hypothesen, jede mit einer Irrtumswahrscheinlichkeit von \(\alpha\) und in allen Fällen sei die \(H_0\) zutreffend. Dann haben wir für jeden einzelnen Test eine Wahrscheinlichkeit von \(1 - \alpha\) die korrekte Entscheidung zu treffen. Wenn wir nun \(m\) unabhängige Hypothesentestungen durchführen, wird die Wahrscheinlichkeit sich korrekt zu entscheiden insgesamt \(m\)-Mal miteinander multipliziert.
\[\begin{equation} P(\text{korrekte Entscheidung}) = (1-\alpha)\cdots(1-\alpha) = (1-\alpha)^m \end{equation}\]
Daraus folgt, dass die Wahrscheinlichkeit mindestens einen \(\alpha\)-Fehler zu machen \(1 - P(\text{korrekte Entscheidung})\) ist.
\[\begin{equation} P(\text{min. }1\text{ Type-I Fehler}) = 1 - (1-\alpha)^m \end{equation}\]
Wird diese Wahrscheinlichkeit gegen die Anzahl der Test abgetragen, ergibt sich der folgende Zusammenhang (siehe Abbildung 23.2).
In Abbildung 23.2 ist zu erkennen, dass die Wahrscheinlichkeit für mindestens einen \(\alpha\) sehr schnell mit der Anzahl der Tests ansteigt. D.h. ab etwa 15 Tests liegt die Wahrscheinlichkeit für ein statistisch signifikantes Ergebnis bei etwas über \(50\%\). Daher, umso mehr Tests in einer Untersuchung durchgeführt werden, umso höher die Wahrscheinlichkeit, dass ein Test anspringt. Um dieses Problem einzufangen, sollte daher für die Anzahl der Tests kontrolliert werden (siehe Rothman (1990) für eine Gegenposition).
Ein Umstand macht die Kontrolle der \(\alpha\)-Fehlerrate allerdings etwas undurchsichtig, da nicht eindeutig ist auf welchem Level die Fehlerkontrolle stattfinden soll. In Abbildung 23.3 ist die typische Struktur eines Experiments dargestellt.
In dem Beispiel sind zwei Endpunkte in dem Experiment geplant gewesen (Primary und Secondary Endpoint). Dies könnte zum Beispiel die Squad und die Hip thrust Leistung in der Stichprobe sein. Beide Endpunkte wurden unter jeweils zwei verschiedenen Faktoren untersucht (z.B. Ernährungsergänzung und Kontrolle). Faktor 1 besitzt drei Faktorstufen. D.h. wir können wiederum drei paarweise Vergleiche zwischen den Mittelwerten (MW) der drei Stufen durchführen. Ähnliches trifft aber auch auf Faktor 2 zu. Bei zwei Endpunkten verdoppelt sich noch einmal die Anzahl der Test bzw. wenn wir noch weitere Endpunkte haben entsprechend mehr. Dies führt dazu, dass wir mehrere Ebenen haben auf denen ein Fehlerkontrollmechanismus ansetzen kann. Die Kontrolle kann entweder auf der Ebene der Endpunkte, der Ebene der Faktoren innerhalb eines Endpunktes oder auf der Ebene der einzelnen Faktorstufen innerhalb eines Faktors innerhalb eines Endpunktes kontrolliert werden. Von oben nach unten werden diese Ebenen als Experiment , family und test wise error rates bezeichnet. Die Kontrolle ist dabei von unten nach oben immer strikter da die Anzahl der zu berücksichtigenden Test immer größer wird. Dazu kommt, dass es leider in der Literatur dazu keinen klaren Konsens gibt auf welcher Ebene die Mehrfachtestung berücksichtigt werden sollte. In den meisten Fällen wird dher auf der test-Ebene kontrolliert.
23.2 Kontraste \(\psi\)
Im folgenden wird ein etwas allgemeineres Rahmenwerk aufgebaut um Vergleiche zwischen Faktorstufen durchzuführen. Dazu wird zunächst das Konzept eines Kontrasts für Mittelwertsvergleiche benötigt.Ein Kontrast \(\psi\) ist dabei ein Spezialfall einer sogenannten Linearkombination . Allgemein hat eine Linearkombination die folgende Form.
\[\begin{equation*} \sum_{i=1}^K c_i \cdot a_i \end{equation*}\]
Es wird eine Summe Elementen \(a_i\) multipliziert mit einem Faktor \(c_i\). Die \(c_i\)s werden als Gewichte bezeichnet. Sei zum Beispiel ein Vektor \(\boldsymbol{a}\) gegeben mit den Elementen \(\boldsymbol{a} = (1,2,3)\) und ein Gewichtsvektor mit den Elementen \(\boldsymbol{c} = (1/2,1/2,1/3)\), dann hat die entsprechende Linearkombination die folgende Form:
\[\begin{equation} \sum_{i=1}^K c_i \cdot a_i = \frac{1}{2}\cdot 1 + \frac{1}{2}\cdot 2 + \frac{1}{3}\cdot 3 \label{eq-ed-lc-ex-01} \end{equation}\]
Ein Kontrast \(\psi\) ist dann ein Spezialfall einer Linearkombination bei dem die Gewichte \(a_i\) sich zu Null addieren. D.h. es gilt:
\[\begin{equation*} \sum_{i=1}^K c_i \cdot a_i \quad \text{mit } \sum_{i=1}^K c_i = 0 \end{equation*}\]
Das Beispiel der Linearkombination in Formel \(\eqref{eq-ed-lc-ex-01}\) stellt somit keinen Kontrast dar, da \(1/2 + 1/2 + 1/3 \neq 0\) gilt. Dagegen wäre die folgende Linearkombination mit \(\boldsymbol{c} = (1/2,1/2,-1)\) ein Kontrast.
\[\begin{equation*} \sum_{i=1}^K c_i \cdot a_i = \frac{1}{2}\cdot 1 + \frac{1}{2}\cdot 2 - \frac{1}{3}\cdot 3 \end{equation*}\]
Tipp
Linearkombination können in R einfach über das Skalarprodukt berechnet werden.
c_i <- c(1/2,1/2,1/3)
a_i <- 1:3
t(a_i) %*% c_i [,1]
[1,] 2.5Was einen Spezialfall des Matrizenprodukts darstellt.
mat <- matrix(1:6, nr=2, byrow=T)
mat [,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6mat %*% c_i [,1]
[1,] 2.5
[2,] 6.523.2.1 Struktur von Kontrasten
Seien nun die Elemente von \(\boldsymbol{a}\) die Modellparamter eines statistischen Models dann erhalten ergibt sich die folgende Definition:
Definition 23.1 (Kontrast) Ein Kontrast \(\psi\) ist eine Linearkombination von Modellparametern deren Koeffizienten sich zu Null addieren.
Nehmen wir zum Beispiel den Fall des CRD Modells mit einem Faktor und \(K\) Faktorstufen.
\[\begin{equation*} Y_{ij} = \mu + \tau_i + \epsilon_{ij} \end{equation*}\]
Die jeweiligen Effekte der Faktorstufen werden durch die Parameter \(\tau_i,i=1,\ldots,K\) repräsentiert. Die allgemeine Form eines Kontrasts \(\psi\) für dieses Modell hat dann die Form:
\[\begin{equation} \psi = \sum_{i=1}^k c_i \tau_i, \quad \text{mit } \sum_{i=1}^k c_i=0 \label{eq-ed-lc-crd-con} \end{equation}\]
In den meisten Fällen werden mittels eines Kontrasts die Mittelwerte \(\mu_i\) der Faktorstufen miteinander verglichen, da die Mittelwert \(\bar{y}_i\) als Schätzer für die Faktorstufeneffekte \(\hat{\tau}_i = \bar{y}_i\) dienen. Der Effekt der \(i\)-ten Stufe setzt sich nach dem Modell zusammen aus \(\mu + \tau_i\) und dieser ist bei einem CRD gleich dem Gruppenmittelwert \(\bar{y}_i\). D.h. ein Schätzer für einen Kontrast \(\hat{\psi}\) hat die Form.
\[\begin{equation} \hat{\psi} = \sum_{i=1}^K c_i (\hat{\mu} + \hat{\tau_i}) = \sum_{i=1}^K c_i \bar{y}_{i.} \end{equation}\]
Unter Verwendung der Eigenschaft, dass ein Kontraste \(\psi\) zu \(0\) aufsummiert gilt:
\[\begin{equation*} \psi=\sum_{i=1}^Kc_i(\mu + \tau_i) = \sum_{i=1}^Kc_i \mu + \sum_{i=1}^K c_i\tau_i=\mu\underbrace{\sum_{i=1}^K c_i}_{=0} + \sum_{i=1}^K c_i\tau_i = \sum_{i=1}^K c_i\tau_i \end{equation*}\]
D.h. über den Vergleich der Mittelwerte \(\bar{y}_i\) erhalten wir den Vergleich der Modellparameter \(\tau_i\).
Durch eine geschickte Kombination von \(c_i\) können nun gewünschte Gruppenvergleiche durchgeführt werden. Im Beispiel der Laufdaten, können mittels des folgenden Schemas von \(c_i\)s die Gruppenvergleich als Kontraste definiert werden (siehe Tabelle 23.1).
| Kontrast | Koffein \((c_1)\) | Placebo \((c_2)\) | Control \((c_3)\) | \(\sum_{i=1}^K c_i\) |
|---|---|---|---|---|
| \(\Delta_{\text{Koffein-Placebo}}\) | \(1\) | \((-1)\) | 0 | 0 |
| \(\Delta_{\text{Koffein-Control}}\) | \(1\) | \(0\) | \((-1)\) | 0 |
| \(\Delta_{\text{Placebo-Control}}\) | \(0\) | \(1\) | \((-1)\) | 0 |
Schauen wir uns konkret anhand des Beispiels an wie die Kontraste berechnet werden. In Tabelle 23.2 sind die Gruppenmittelwerte aus dem Koffeinbeispiel abgebildet.
| Gruppe | \(\bar{y}_{i.}\) |
|---|---|
| Koffein | 1914.2 |
| Placebo | 1924.7 |
| Control | 1927.5 |
Wollen wir nun einen Vergleich zwischen den Gruppen Koffein und Placebo bestimmen, dann wählen wir aus Tabelle 23.1 den ersten Kontrast \(c_{\text{K-C}} = (1,-1,0)\) aus. Angewendet auf die Mittelwert der Gruppen folgt daraus.
\[\begin{equation*} \hat{\psi}_{\text{K-C}} = \sum_{i=1}^k c_i \bar{y}_{i.} = 1 \cdot 1914.2 + (-1) \cdot 1924.7 + 0 \cdot 1927.5 = -10.5 \end{equation*}\]
D.h. wir beobachten einen Unterschied von \(-10.5s\) zwischen den beiden Gruppen. Die beiden verbleibenden paarweisen Vergleichen können parallel durchgeführt werden.
Allgemein wird zwischen zwei Arten von Kontrasten unterschieden: Einfache und komplexe Kontraste. Einfache Kontraste sind dabei die paarweisen Vergleiche und alle anderen Arten von Kontrasten als komplexe Kontraste bezeichnet werden. Im Beispiel ist zu erkennen, dass bei paarweisen Kontrasten immer ein Kontrastgewicht \(c_i = 1\) gesetzt ist und ein weiterer auf \(c_j = -1\) gesetzt wird, während alle anderen Gewichte \(c_l = 0, \forall l \neq i,j\) sind.
Erstellen wir nun einen komplexen Kontrast. Zum Beispiel möchten wir untersuchen ob es einen Unterschied zwischen der Kontrollgruppe und den beiden Gruppen die eine Pille bekommen gibt. D.h. wir wollen Untersuchen ob die Gabe einer Pille unabhängig davon ob in der Pille ein Wirkstoff (Koffein) oder nicht (Placebo) gewesen ist. Mittels eines Kontrasts können wir diesen Vergleich wie folgt durchführen.
\[\begin{equation*} c_1 = \frac{1}{2}, c_2 = \frac{1}{2}, c_3 = -1 \end{equation*}\]
D.h. wir Mitteln die Mittelwerte von Koffein und Placebo und vergleichen diesen Mittelwert mit dem Mittelwert der Kontrollgruppe. Wenn wir die Summe der Gewichte bilden erhalten wir \(\frac{1}{2} + \frac{1}{2} - 1 = 0\). Es handelt sich also immer noch um eine Kontrast. Angewendet auf das Beispiel erhalten wir:
\[\begin{equation*} \hat{\psi} = \sum_{i=1}^k c_i \bar{y}_{i.} = \frac{1}{2} \cdot 1914.2 + \frac{1}{2} \cdot 1924.7 + (-1) \cdot 1927.5 = -8.05 \end{equation*}\]
D.h. der Effekt der Pillengabe unabhängig von Wirkstoff ist \(-8.05s\).
Zusammengefasst ermöglichen Kontraste beliebige Vergleiche zwischen mehreren Faktorstufen durchzuführen um Modellparameter miteinander zu vergleichen. Die üblichen paarweisen Vergleich bilden hierbei einen Spezialfall (einfache Kontraste) der allgemeineren komplexen Kontrasten. Um nun zu bewerten ob ein beobachteter Kontrast \(\psi_i\) auch wirklich bedeutsam ist, wird nun noch Information über die Varianz \(Var(\psi_i)\) benötigt.
23.2.2 Varianz von Kontrasten
Um die Varianz von Kontrasten herzuleiten hier noch einmal kurz zur Wiederholung der Standardfehler des Mittelwerts
\[\begin{equation} \sigma_{\bar{x}} = \frac{\sigma}{\sqrt{n}} \end{equation}\]
Die Varianz des Mittelswerts ist dementsprechend das Quadrat des Standardfehlers \(\sigma_{\bar{x}}^2\). Sei beispielsweise der Mittelwert in der Stichprobe unter der \(i\)-ten Faktorstufe \(\bar{Y}_i\) mit der Stichprobengröße \(n_i\), dann gilt:
\[\begin{equation*} Var(\bar{Y}_i) = \frac{\sigma^2}{n_i} \end{equation*}\]
Als zweites Wiederholungsitem die Varianz einer Linearkombination von unabhängigen Zufallsvariablen berechnet sich nach (siehe Formel \(\eqref{eq-stats-hypo-sum-var}\))
\[\begin{equation*} Var\left(\sum_{i=1}^n a_i X_i\right) = \sum_{i=1}^n a_i^2 Var(X_i) \end{equation*}\]
Wichtig für die weiteren Herleitungen ist insbesondere das Beispiel Beispiel 12.2. Angewendet auf einen Kontrast \(\psi\) folgt für die Varianz des Kontrast unter der Annahme das die eingehenden Mittelwerte \(\bar{Y}_{i.}\) alle die gleiche Varianz \(\sigma^2\) haben (unsere übliche Annahme der Homoskedastizität), also \(Var(\bar{Y}_{i.}) = \sigma^2\).
\[\begin{equation} \begin{aligned} Var(\psi) &= Var\left(\sum c_i \bar{Y}_{i.}\right) = \sum c_i^2 Var(\bar{Y}_{i.}) \\ &= \sum c_i^2(\sigma^2/n_i) = \sigma^2\sum(c_i^2/n_i) \end{aligned} \label{eq-ed-lc-var-con-pop} \end{equation}\]
D.h. die Kontrastgewichte \(c_i\) werden quadriert und durch die jeweilige Stichprobengröße \(n_i\) der \(i\)-ten Faktorstufe geteilt. Wie immer schätzen wir die Varianz \(\sigma^2\) anhand der Stichprobe mittels des \(MSE\) bzw. in der Terminologie des CRD \(MS_{\text{within}}\). Zur Erinnerung der \(MSE\) ist die Quadratsumme der Residuen \(e_i\) geteilt durch die Stichprobengröße \(N\) minus der Anzahl der Parameter \(p\). Daraus folgt für die den Schätzer der Varianz des Kontrasts \(\widehat{Var}(\psi)\).
\[\begin{equation} \widehat{Var}(\psi) = \widehat{Var}\left(\sum c_i \bar{y}_{i.}\right) = MS_w\sum (c_i^2/n_i) \label{eq-ed-lc-var-con} \end{equation}\]
Letztendlich ist dies die Formel die wir in der Einführung zu Experimentellen Design verwendet haben, um die Verteilungen der Stichprobengrößen in die beiden Gruppen zu verteilen. Um nun von der Varianz auf den Standardfehler des Kontrasts \(s_{\psi}\) zu kommen, ziehen wir wie immer die Wurzel aus der Varianz:
\[\begin{equation} s_{\psi} = \sqrt{MS_w\sum (c_i^2/n_i)} \label{eq-ed-lc-se-contrast} \end{equation}\]
\(n_i\) = Stichprobengröße in Faktorstufe \(i\)
\(\hat{\sigma}^2 = MS_w = MSE\)
\(\hat{\sigma}^2 = MS_w = MSE\)
Angewendet auf den Kontrast bei Paarvergleichen folgen für den Vergleich der Faktorstufen \(i\) und \(j\) dargestellt mittels der Kontrastgewichte \(c_i\) die folgenden Terme \(c_i^2 = (1)^2, c_j^2 = (-1)^2\) während alle anderen \(c_l = 0\) sind.
\[\begin{equation} Var(\psi) = Var\left(\tau_i - \tau_j\right) = \sigma^2\left(\frac{1}{n_i} + \frac{1}{n_j}\right) = \end{equation}\]
Wenn wir ein Design haben, bei dem die Stichprobengrößen unter allen Faktorstufen gleich groß sind, also \(n_i = n_j = n\), gilt, dann folgt mit der Festsetzung \(2n = N\) für den Standardfehler des Kontrasts \(\sigma_{\psi}\) in der Population:
\[\begin{equation*} \sigma_{\psi} = s_{\Delta} = \sqrt{\sigma^2\left(\frac{1}{N/2} +\frac{1}{N/2}\right)}=\sqrt{\sigma^2\frac{2}{N/2}} = \sigma\sqrt\frac{4}{N}=\sigma \frac{2}{\sqrt{N}} \end{equation*}\]
Beziehungsweise in der Stichprobe dann entsprechend.
\[\begin{equation} s_{\psi} = s_{\Delta} = \sqrt{MSE}\frac{2}{\sqrt{N}} \end{equation}\]
D.h. der Kontrast für den Vergleich zwischen zwei Mittelwerten hat im Vergleich zu Standardfehler des Mittelwerts die doppelte Größe. Insgesamt erhalten wir aber wieder den gleichen Wurzelzusammenhang der dazu führt, dass der Standardfehler für den Kontrast mit zunehmender Stichprobengröße stetig abnimmt, aber die Abnahme mit größer werdenden Stichproben immer kleiner wird. Anders herum bei kleinen Stichprobengrößen führt jeder Hinzunahme von Stichproben zu einer relativ stärkeren Verkleinerung des Standardfehlers und damit zu einer Erhöhung der Power.
23.2.3 Hypothesentests für Kontraste
Ein Hypothesentest für den Kontrast kann nun wieder nach dem üblichen Muster konstruiert werden indem wir den Kontrast \(\hat{\psi}\) durch seinen Standardfehler teilen \(s_{\psi}\). Die entsprechenden Hypothesen sind:
\[\begin{equation} \begin{split} H_{0,\psi}: \sum_{i=1}^k c_i \tau_i = 0 \\ H_{1,\psi}: \sum_{i=1}^k c_i \tau_i \neq 0 \end{split} \end{equation}\]
Die Teststatistik \(t\) ist wiederum nichts anderes als ein Verhältnis von Effektstärke normiert durch deren Standardfehler. Wenn dieses Verhältnis dann wieder größer als ein kritischer Wert \(w\) ist, kann die entsprechende \(H_0\) abgelehnt werden. Es folgt somit:
\[\begin{equation} \text{Lehne }H_0\text{ ab, wenn} \left|\frac{\sum_{i=1}^K c_k \tau_i}{\sqrt{MSE\sum_{i=1}^K c_i^2/n_i}}\right| > w \label{eq-ed-lc-w} \end{equation}\]
Der kritische Wert \(w\) ist somit letztendlich wieder nichts anders als die Quantile \(q_{\alpha}\) einer entsprechenden theoretischen Verteilung. Welche theoretische Verteilung gewählt werden muss, hängt nun davon welche Art der \(\alpha\)-Kontrolle angewendet werden soll. Im folgenden werden verschiedene Methoden vorgestellt \(\alpha\)-Kontrolle durchzuführen und entsprechend sollten die Überlegungen zur Auswahl so durchgeführt werden, dass eine möglichst hohe Power für das gewählte Design erreicht werden kann.
23.2.4 Konfidenzintervalle von Kontrasten
Nach dem uns bekannten Muster können wir auch Konfidenzintervalle für einen Kontrast \(\psi\) erstellen. Das Berechnungsmuster ist dabei immer das Gleiche:
\[\begin{gather*} \text{CI}(\hat{\psi}) = \sum_i c_i \hat{\tau}_i \pm w \sqrt{\widehat{Var}\left(\sum c_i \hat{\tau}_i\right)} = \hat{\psi} + w \times s_{\psi}\\ \text{estimate } \pm (\text{kritischer Wert}) \times (\text{Standardfehler}) \end{gather*}\]
Der Kontrastschätzer \(\hat{\psi}\) bleibt dabei gleich und es ändert sich je nach gewählter Methode der kritische Wert \(w\).
23.3 Einteilung der Mehrfachvergleiche
Mehrfachvergleiche können zunächst in zwei Arten unterschieden werden. Vergleiche die vor dem Experiment geplant wurden und Vergleiche die nach einsehen der Daten auf Basis der Ergebnisse durchgeführt werden. Vergleiche die vor der Durchführung des Experiments festgelegt werden, werden als pre-planned bezeichnet. Vergleiche die dagegen nach der Durchführung des Experiments nach der Einsicht der Daten durchgeführt werden, werden als post-hoc Vergleiche bezeichnet. Zwischen diesen beiden Vergleichen wird unterschieden, da bei den post-hoc dadurch, das die Daten schon bekannt sind, die Vergleiche durch die Daten beeinflusst werden. Daher muss in diesem Fall die Kontrolle strikter sein, als wenn die Vergleiche im Voraus geplant werden.
Es gibt mittlerweile zahllose Arten die Kontrolle bei Mehrfachvergleiche durchzuführen. Daher besprechen wir hier nur eine kleine Auswahl der üblichsten Methoden. In Tabelle 23.3 sind die Bonferroni, Tukey, Scheffé, Dunnet und die FisherLSD Korrekturmethoden inklusive der Art der Kontraste tabellarisch dargestellt.
| Name | Zeitpunkt | Kontraste | Kontrolliert |
|---|---|---|---|
| Bonferroni | pre-planned | einfache und komplexe | Ja |
| Tukey | pre-planned | alle paarweisen | Ja |
| Scheffé | post-hoc | einfache und komplexe | Ja |
| Dunnet | pre-planned | paarweise TRT gegen CON | Ja |
| FisherLSD | post-hoc | einfache und komplexe | Nein1 |
Die Bonferroni-Korrektur ist in der Literatur weit verbreitet, da sie relativ einfach durchzuführen ist. Die \(\alpha\)-Korrektur wird erreicht, indem ein korrigiertes Signifikanzlevel \(\alpha^*\) nach der folgenden Formel berechnet wird.
\[\begin{equation*} \alpha^* = \alpha / m \end{equation*}\]
\(m\) ist die Anzahl der Kontraste. Die Bonferroni-Korrektur kontrolliert dadurch das Gesamt-\(\alpha\) und ist anwendbar bei einfachen und/oder komplexen Kontraste. Allerdings, müssen diese Kontraste pre-planned sein, d.h. die Vergleiche müssen vor Einsicht in die Daten bestimmt worden sein. Die Bonferroni-Methode ist üblichweise eher auf der konservativen Seite, vor alle da mit der Anzahl der Tests \(\alpha^*\) relativ schnell sehr klein wird und dementsprechend die Power abnimmt. Der kritische Wert in Formel \(\eqref{eq-ed-lc-w}\) wird unter der Bonferroni-Korrektur mittels der \(t\)-Verteilung mit \(N-K\) Freiheitsgraden und \(1-\alpha/(2m)\)-Quantile bestimmt.
\[\begin{equation} w_{\text{Bonferroni}} = t_{N-K,1-\alpha/(2m)} \end{equation}\]
Eine weitere weit verbreitete Methode ist die Tukey-Pairwise difference Methode . Diese Methode ist optimiert darauf alle paarweisen Vergleiche durchzuführen und daher nur für einfache Kontraste anwendbar. Durch die Optimierung sind die Konfidenzintervalle unter der Tukey-Methode überlicherweise schmaler als dies unter der Bonferroni-Methode der Fall wäre. Die Tukey-Methode ist ebenfalls pre-planned und im Falle von ungleichen Stichproben in den Faktorstufen ist die Methode nur asymptotisch exact. Unter der Tukey-Methode berechnet sich der kritische Wert \(w\) anhand eine speziellen Verteilung, der Studentized range statistics.
\[\begin{equation} w_{\text{Tukey}} = q_{K,N-K,\alpha/\sqrt{2}} \label{eq-ed-lc-tukey-w} \end{equation}\]
In R kann die Quantile der Studentized range distribution aus Formel \(\eqref{eq-ed-lc-tukey-w}\) mittels der Funktion qtukey() bestimmen werden.
Warnung
Die Tukeykorrektur ist nur für balancierte Vergleiche, d.h. gleiche Stichprobengrößen in allen Gruppen, hergeleitet und für ungleiche Stichprobengrößen nur annährend exakt .
Die Scheffé-Methode ist ebenfalls sehr oft in der Literatur anzutreffen, da sie post-hoc, also nach Einsicht der Daten, angewendet werden kann. Mit der Scheffé-Methode können auch alle Arten von Vergleichen, einfache wie auch komplexe Vergleiche, durchgeführt werden. Der kritische Wert \(w\) wird bei Scheffé-Methode mittels der \(F\)-Verteilung mit \(K-1\) und \(N-k\) Freiheitsgraden bei \(1-\alpha\) ermittelt. Wobei noch ein Korrekturfaktor \(K-1\) multipliziert werden muss. Dabei ist die Anzahl der Vergleiche irrelevant, da der kritische Wert alle möglichen Vergleiche einbezieht.
\[\begin{equation} w_{\text{Scheffé}} = \sqrt{(K-1)F_{K-1,N-K,1-\alpha}} \end{equation}\]
Dementsprechend kann die Scheffé-Korrektur eigentlich immer angewendet werden, egal ob pre-planned oder post-hoc und unabhängig von der Art der Vergleichen ohne das das Risiko einer fehlerhaften Korrektur eingegangen wird.
Um noch einmal die Unterschiede der verschiedenen Methoden stärker zu betonen, sind in Tabelle 23.4 die verschiedene kritische Werte \(w\) der drei Methoden für verschiedenen Anzahl \(K\) von Faktorstufen bzw. der resultierenden paarweisen Vergleichen aufgelistet.
| Stufen | \(\alpha_{PC}\) | Tukey | Bonferroni | Scheffé |
|---|---|---|---|---|
| 2 | 4.75 | 4.75 | 4.75 | 4.75 |
| 3 | 4.75 | 7.12 | 7.73 | 7.77 |
| 4 | 4.75 | 8.81 | 9.94 | 10.47 |
| 5 | 4.75 | 10.16 | 11.75 | 13.04 |
| 6 | 4.75 | 11.28 | 13.31 | 15.53 |
| 7 | 4.75 | 12.25 | 14.69 | 17.98 |
Tabelle 23.4 zeigt, dass der Tukey durchgehend bei den drei Methoden die kleinsten kritischen Werte \(w\) verwendet. Dadurch wird entsprechend das Konfidenzintervall auch schmaler. Für diesen Fall liegt die Bonferroni-Methode immer zwischen Tukey und Scheffé die immer den größten kritischen Wert \(w\) hat. Die kritischen Werte unter der Scheffé werden relativ schnell größer. Da die Größe des kritischen Werts \(w\) die Power beeinflusst, können wir anhand der Werte ableiten, dass die Scheffé-Methode entsprechend die geringste Power hat.
Zwei weitere Methoden sind noch die Methode nach Dunnet und Fisher least significant difference (FisherLSD) . Die Dunnet Methode ist ebenfalls pre-planned und ist darauf optimiert ein Kontrollkondition mit den anderen Konditionen zu vergleichen. Die FisherLSD-Methode führt keine Anpassung durch sondern ist ein Spezialfall wenn der Faktor \(K = 3\) Stufen hat und im ersten Schritt ein statistisch Signifikanter \(F\)-Wert gefunden wurde. Die drei individuellen Tests werden dann als normale t-Tests durchgeführt mit dem jeweils unkorrigierten \(\alpha\). In diesem Fall, wenn der paarweise Vergleich nach einen signifikanten \(F\)-Test bei \(K=3\) durchgeführt wird, ist eine gesonderte \(\alpha\)-Kontrolle nicht mehr notwendig.
Zusammenfassend sollte vor der Durchführung des Experiment schon eine klare Idee formuliert werden, welche Vergleiche von primären Interesse sind und entsprechend in einer Prä-Registrierung festgehalten werden. Dadurch lässt sich die Power maximieren.
23.4 Mehrfachvergleiche in R
23.4.1 Das package emmeans
In R können Mehrfachvergleiche über verschiedene Pakete durchgeführt werden. Wir beschränken uns hier hauptsächlich auf das Paket emmeans, da es sehr flexibel ist und mit einer Vielzahl von verschiedenen Modellen über das allgemeine lineare Modell hinaus zusammenarbeitet. Alternativen zu emmeans sind multcomp und gmodels.
Bei der Verwendung von emmeans ist immer zweistufiger Prozess durchzuführen. Zunächst müssen die Zellmittelwerte mit der Funktion emmeans(<MODEL>, ~<FAKTOR>) berechnet werden. Anschließen erfolgt die Berechnung der Vergleiche entweder mittels mittels der Funktion pairs() bei paarweisen Vergleichen oder mit contrast() für beliebige Vergleiche.
Wollen wir zum Beispiel mittels Bonferroni-Korrektur einen einfachen Vergleich zwischen Koffein und Placebo sowie einen komplexen Vergleich der Mittelwerte von Koffein und Placebo gegen Kontrolle durchführen, dann verwenden wir die folgenden Befehle nachdem wir das package emmeans geladen haben:
mod_em <- emmeans(mod_aov, specs=~Gruppe)
contrast(mod_em,
adjust='bonferroni',
method = list(
"Koffein vs. Placebo" = c(1, -1, 0),
"Tablette vs. CON" = c(1/2, 1/2, -1)
),
infer=TRUE) contrast estimate SE df lower.CL upper.CL t.ratio p.value
Koffein vs. Placebo -10.5 2.09 57 -15.3 -5.67 -5.014 <.0001
Tablette vs. CON -8.1 1.81 57 -12.3 -3.93 -4.476 0.0001
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 2 estimates
P value adjustment: bonferroni method for 2 tests Im ersten Schritt haben wir die Mittelwerte mittels emmeans() berechnet. Dazu haben wir das gefitte Modell mod_aov übergeben und als zweiten Parameter specs übergeben wir eine Formel für welchen Faktor wir die Mittelwerte berechnen wollen. Wie wir später bei mehrfaktoriellen Modellen sehen gibt es hier verschiedene Möglichkeiten. Bei diesem Beispiel sind hier aber nur die Mittelwerte für die Faktorstufen von Gruppe möglich. Schauen wir uns kurz an, welchen Output `emmeans()´ produziert.
mod_em Gruppe emmean SE df lower.CL upper.CL
Koffein 1914 1.478 57 1911 1917
Placebo 1925 1.478 57 1922 1928
Control 1928 1.478 57 1925 1931
Confidence level used: 0.95 Letztendlich haben wir hier die Mittelwerte der drei verschiedenen Gruppen die anhand des Modells berechnet werden.
Warnung
In diesem Fall sind die Mittelwerte \(\bar{y}_i\) der Rohdaten gleich denjenigen des Modell des Modells, da es sich um ein balancierten Datensatz handelt, d.h. in allen Gruppen sind die Stichproben gleich. Wenn dies nicht der Fall ist, dann müssen die geschätz ten Mittelwerte \(\hat{y}_i\) anhand des Modells nicht unbedingt gleich zu den Mittelwerte \(\bar{y}_i\) der Rohdaten sein. Dies trifft vor allem zu, wenn die Modell komplexer werden. Die Werte von emmeans() werden anhand des Modells geschätzt und werden daher auch als expected marginal means bezeichnet (daher auch der Name des packages e(expected)m(arginal)means).
Nachdem wir die Mittelwerte mit emmeans() berechnet haben, können wir nun die Kontraste mit der Funktion contrast() aufsetzen. Dazu wird das emmeans-Objekt als erster Parameter übergeben. Über den Parameter adjust wird die Korrekturmethode angegeben. In diesem Fall bonferroni. Anschließend übergeben wir eine Liste der Kontraste die berechnet werden sollen. Hierzu muss die Abfolge der Mittelwerte laut der Reihenfolge in der Ausgabe in emmeans() verwendet werden, die üblicherweise die Gleiche wie auch bei lm() ist. Die Kontraste werden als Liste mit <NAME> = <KONTRAST> übergeben. Der Name ist willkürlich und sollte so gewählt werden, dass der eben Wert möglichst einfach interpretiert werden kann. Der Parameter infer kontrolliert ob direkt ein Konfidenzintervall berechnet werden soll. Alternativ könnte auch wie folgt vorgegangen werden:
contrast_list <-list(
"Koffein vs. Placebo" = c(1, -1, 0),
"Tablette vs. CON" = c(1/2, 1/2, -1)
)
mod_contrast <- contrast(mod_em,
adjust = 'bonferroni',
method = contrast_list)
mod_contrast contrast estimate SE df t.ratio p.value
Koffein vs. Placebo -10.5 2.09 57 -5.014 <.0001
Tablette vs. CON -8.1 1.81 57 -4.476 0.0001
P value adjustment: bonferroni method for 2 tests Mit nachfolgender Berechnung der Konfidenzintervalle.
confint(mod_contrast) contrast estimate SE df lower.CL upper.CL
Koffein vs. Placebo -10.5 2.09 57 -15.3 -5.67
Tablette vs. CON -8.1 1.81 57 -12.3 -3.93
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 2 estimates Die Ergebnisse sind natürlich in beiden Folgen gleich und zeigen einen statistisch signifikanten Unterschied zwischen den Koffein und der Placebo Kondition, sowie einen statistisch signifikanten Unterschied zwischen der Tablettengabe unabhängig vom Wirkstoff und der Kontrollkondition.
Möchten wir dagegen alle paarweisen Vergleiche beispielsweise mittels der Tukey-Methode berechnen können wir die pairs() Methode auf die mittels emmeans() gefitten Faktormittelwerte anwenden.
pairs(mod_em, adjust="tukey") contrast estimate SE df t.ratio p.value
Koffein - Placebo -10.48 2.09 57 -5.014 <.0001
Koffein - Control -13.34 2.09 57 -6.383 <.0001
Placebo - Control -2.86 2.09 57 -1.369 0.3639
P value adjustment: tukey method for comparing a family of 3 estimates Die Tukey Methode wird standardmäßig von pairs() verwendet, wenn keine Argument für adjust angegeben wird. Alternativ geht es auch über die contrast() Funktion wenn für method die Zeichenkette pairwise übergeben wird.
contrast(mod_em, method='pairwise') contrast estimate SE df t.ratio p.value
Koffein - Placebo -10.48 2.09 57 -5.014 <.0001
Koffein - Control -13.34 2.09 57 -6.383 <.0001
Placebo - Control -2.86 2.09 57 -1.369 0.3639
P value adjustment: tukey method for comparing a family of 3 estimates Soll die Korrektur mittels der Scheffé-Methode durchgeführt werden, kann dies über adjust=scheffe erreicht werden.
contrast(mod_em,
adjust='scheffe',
method = contrast_list) contrast estimate SE df t.ratio p.value
Koffein vs. Placebo -10.5 2.09 57 -5.014 <.0001
Tablette vs. CON -8.1 1.81 57 -4.476 0.0002
P value adjustment: scheffe method with rank 2 Um die Dunnettkorrektur zu verwenden, entsprechend contrast() mit dem Parameter method='trt.vs.ctrl'. Hier kann zusätzlich über den Parameter ref explizit angeben werden, welche Faktorstufe die Kontrollbedingung darstellt.
contrast(mod_em, method='trt.vs.ctrl', ref='Control') contrast estimate SE df t.ratio p.value
Koffein - Control -13.34 2.09 57 -6.383 <.0001
Placebo - Control -2.86 2.09 57 -1.369 0.3029
P value adjustment: dunnettx method for 2 tests Wenn wir die FisherLSD-Methode verwenden wollen, dann kann beispielweise wieder pairs() mit dem Parameter adjust='none' verwendet werden.
pairs(mod_em, adjust='none', infer=T) contrast estimate SE df lower.CL upper.CL t.ratio p.value
Koffein - Placebo -10.48 2.09 57 -14.66 -6.29 -5.014 <.0001
Koffein - Control -13.34 2.09 57 -17.53 -9.16 -6.383 <.0001
Placebo - Control -2.86 2.09 57 -7.05 1.32 -1.369 0.1765
Confidence level used: 0.95 Wie schon erwähnt ist das package emmeans extrem flexibel und bietet unzählige Möglichkeiten verschiedene Mehrfachvergleiche zu spezifizieren und entsprechend zu berechnen. Daher lohnt sich die intensive Auseinandersetzung mit der Dokumentation in jedem Fall. Wahrscheinlich geht ein Teil des Erfolges von R in den letzten Jahren auf die Funktionalität dieses packages zurück.
23.4.2 Cohen’s d für Mehrfachvergleiche
Über das Paket emmeans können wir auch eine Cohen’s D Typ Effektstärke mittels der Funktion eff_size() berechnen. Allerdings müssen wir als Parameter die zu verwendende Standardabweichung sowie die Freiheitsgrade selbst angeben. Dies können aber über die Funktionen sigma() und df.residual() aus dem gefitteten lm-Objekt extrahiert werden.
eff_size(mod_em,
sigma=sigma(mod_aov),
edf=df.residual(mod_aov)) contrast effect.size SE df lower.CL upper.CL
Koffein - Placebo -1.586 0.349 57 -2.29 -0.886
Koffein - Control -2.019 0.368 57 -2.76 -1.281
Placebo - Control -0.433 0.319 57 -1.07 0.206
sigma used for effect sizes: 6.609
Confidence level used: 0.95 Standardmäßig werden paarweise Vergleiche berechnet, aber über den Parameter method können die Effektstärken auch für beliebige Kontraste berechnet werden.
eff_size(mod_em,
sigma=sigma(mod_aov),
edf=df.residual(mod_aov),
method = contrast_list) contrast effect.size SE df lower.CL upper.CL
Koffein vs. Placebo -1.59 0.349 57 -2.29 -0.886
Tablette vs. CON -1.23 0.297 57 -1.82 -0.631
sigma used for effect sizes: 6.609
Confidence level used: 0.95 23.5 Bestimmung des Replikationsfaktors anhand des Konfidenzintervalls
Im Kapitel zum CRD haben wir den Replikationsfaktor für die experimentellen Einheiten anhand der Effektstärke bestimmt. Ein weitere Möglichkeit den Replikationsfaktor zu bestimmen besteht mittels der Kontraste bzw. deren Konfidenzintervallen. Die angesteuerte Größe ist hierbei die Breite des Konfidenzintervalls. Der direkteste Weg ist hierbei den minimal signifikanten Unterschied (minimal significant difference: msd) . Der Ansatz basiert im Grunde genommen auf der Dualität von Hypothesentests und Konfidenzintervallen.
In Abbildung 23.4 ist noch einmal der Zusammenhang zwischen dem Hypothesentest unter der \(H_0\) und der Breite des Konfidenzintervalls dargestellt. Wenn das Konfidenzintervall die \(H_0\) nicht enthält, dann gehen wir von einem statistisch signifikantem Ergebnis aus. Daher wenn wir von einem gegebenem Effekt, im vorliegenden Fall einem Kontrast \(\psi\) ausgehen und die Residualvarianz \(\sigma\) kennen, können wir über den Standardfehler \(s_{\psi}\) die notwendige Stichprobengröße bestimmen. Der minimale signifikante Unterschied kann über die folgende Formel berechnet werden.
\[\begin{equation} msd = w \sqrt{\widehat{Var}\left(\sum c_i \hat{\tau}_i\right)}=w \times s_{\psi} \end{equation}\]
Wir wollen zum Beispiel den Effekt eines Nahrungsergänzungsmittels auf die Laufleistung bei einem Mitteldistanzwettkampf ähnlich wie bei unserem Koffeinexperiment durchführen. Wir wollen das gleiche Design mit Interventions-, Placebo- und Kontrollgruppe durchführen (\(K = 3\)). Wir gehen davon aus, dass ein Unterschied zwischen der Placebogruppe und der Interventionsgruppe von \(5\)s mindestens gefunden werden sollte um als praktisch relevant zu gelten. Damit gilt
\[\begin{equation*} msd = 5 \end{equation*}\]
Im Koffeinbeispiel hatten wir für \(\sigma^2\) einen Wert von \(\hat{\sigma^2} = MSE \approx 30\)s bestimmt. Daraus folgt für einen paarweisen Vergleich mit gleicher Stichprobengröße \(n\) in beiden Gruppen:
\[\begin{equation*} msd = w \sqrt{30\left(\frac{1}{n} + \frac{1}{n}\right)} = w\sqrt{30\frac{2}{n}} \leq 5 \end{equation*}\]
Wir wissen a-priori das wir nur paarweise Vergleiche durchführen wollen, daher wissen wir, dass wir die Quartile mittels der Tukey-Methode bestimmen werden. D.h. der kritische Wert \(w\) berechnet sich nach Formel \(\eqref{eq-ed-lc-tukey-w}\).
\[\begin{equation*} w = q_{K,N-K,\alpha}/\sqrt{2} \end{equation*}\]
Damit erhalten wir:
\[\begin{align*} \frac{q_{K,N-K,\alpha}}{\sqrt{2}} \sqrt{30\frac{2}{n}} &= q_{K,N-K,\alpha} \sqrt{\frac{30}{n}} \leq 5 = msd \\ \Leftrightarrow q_{K,N-K,\alpha}^2 \frac{30}{n} &\leq 25 \\ \Leftrightarrow q_{K,N-K,\alpha}^2 \leq \frac{25}{30}n &= \frac{5}{6}n \\ \end{align*}\]
D.h. das Quadrat der studentized range quantile muss kleiner als \(\frac{5}{6}n\) sein, damit das Konfidenzintervall schmal genug ist. Da sich der Wert von \(q_{K,N-K,\alpha}\) mit der Stichprobengröße ändert führen wir eine try-and-error Methode durch und berechnen die Breite des Konfindenzintervalls für verschiedene Wert von \(n\) mittels qtukey()-Funktion.
K <- 3
n <- 10:20
q_t <- qtukey(0.975,K,K*n-K)**2
n_w <- 5/6*nWenn wir die Werte in einer Tabelle abtragen erhalten wir die folgende Liste (siehe Tabelle 23.5).
| n | \(q_{n,Kn-K,\alpha}^2\) | \(5/6n\) |
|---|---|---|
| 10 | 15.58 | 8.33 |
| 11 | 15.36 | 9.17 |
| 12 | 15.18 | 10.00 |
| 13 | 15.04 | 10.83 |
| 14 | 14.92 | 11.67 |
| 15 | 14.82 | 12.50 |
| 16 | 14.73 | 13.33 |
| 17 | 14.65 | 14.17 |
| 18 | 14.58 | 15.00 |
| 19 | 14.52 | 15.83 |
| 20 | 14.47 | 16.67 |
In Tabelle 23.5 sehen wir, dass wir eine Stichprobengröße von \(n = 18\) also insgesamt \(N = 54\) Athleten brauchen um die gewünschte Präzision für das Konfidenzintervall zu erreichen.
23.6 Vollständige Dokumentation eines CRDs
Hier einmal beispielhaft Dokumentation eines CRD inklusive verwendeter Mehrfachvergleiche.
Eine einfaktorielle ANOVA mit dem Faktor Gruppe ergabe einen statistisch signifikanten Haupteffekt für Gruppe \(F(2, 57) = 22,6, p < 0,001\). Überprüfung auf Varianzgleichheit zwischen den Gruppen mittels eines Levene-Tests deutete auf keine Verletzung der Voraussetzungen hin, \(F(2, 57) = 0,38, p = 0,69\). Daher wird die \(H_0\), das kein Unterschied zwischen den Gruppen besteht, abgelehnt. Die beobachtete Effektstärke \(\omega^2 = 0,42\), CI\(95\%[0,22, 0,56]\) ist als großer Effekt zu interpretieren. Pre-planned Paarweisetestung mittels Tukey-Korrektur deutete auf statistisch signifikate Unterschiede zwischen den Gruppen Koffein und Placebo \(z = -10.4\), CI95%\([-15,5, -5,5], p < 0,001\), und Koffein und Kontrolle, \(z = -13,3\), CI95\(\%[-18,4, -8,3], p < 0,001\), hin. Insgesamt deuten die Ergebnisse daher darauf hin, dass die Gabe von Koffein zu einer bedeutsamen Leistungssteigerung \((>5-10s)\) in der beobachteten Untersuchungsgruppe geführt hat.
23.7 Zum Nacharbeiten
Weiter Informationen zu Mehrfachvergleiche und den möglichen Problemen findet ihr in beispielsweise in Feise (2002), Rothman (1990).
Außer wenn \(K=3\)↩︎