Post hoc analyse: Meervoudig Vergelijken van Gemiddelden

Naïeve methode

In het eerste deel van dit hoofdstuk hebben we een \(F\)-test besproken die gebruikt kan worden voor het testen van

\[H_0: \mu_1=\cdots = \mu_g \text{ versus } H_1: \text{niet } H_0.\]

Dus als de nulhypothese verworpen wordt, dan wordt besloten dat er minstens twee gemiddelden verschillen van elkaar. De methode stelt ons echter niet in staat om te identificeren welke gemiddelden van elkaar verschillen.

Een eerste, maar naïeve benadering van het probleem bestaat erin om de nulhypothese op te splitsen in partiële hypotheses

\[H_{0jk}: \mu_j=\mu_k \text{ versus } H_{1jk}: \mu_j \neq \mu_k\]

en deze partiële hypotheses te testen met two-sample \(t\)-testen. Voor het vergelijken van groep \(j\) met groep \(k\) wordt de klassieke two-sample \(t\)-test onder de veronderstelling van homoscedasticiteit gegeven door

\[T_{jk} = \frac{\bar{Y}_j-\bar{Y}_k}{S_p\sqrt{\frac{1}{n_j}+\frac{1}{n_k}}} \sim t_{n-2}\]

waarin \(S_p^2\) de gepoolde variantieschatter is,

\[S_p^2 = \frac{(n_j-1)S_j^2 + (n_k-1)S_k^2}{n_j+n_k-2}\]

met \(S_j^2\) en \(S_k^2\) de steekproefvarianties van respectievelijk de uitkomsten uit groep \(j\) en \(k\).

In een ANOVA context wordt echter verondersteld dat in alle \(g\) groepen de variantie van de uitkomsten dezelfde is (de residuele variantie \(\sigma^2\)). Indien we dus \(S_p^2\) gebruiken, dan is dit niet de meest efficiënte schatter omdat deze niet van alle data gebruik maakt[47]. We kunnen dus efficiëntie winnen door MSE te gebruiken. Ter herinnering, MSE kan geschreven worden als

\[\text{MSE}= \sum_{j=1}^g \frac{(n_j-1)S_j^2}{n-g}.\]

De \(t\)-testen voor het twee-aan-twee vergelijken van alle gemiddelden worden dus best gebaseerd op

\[T_{jk} = \frac{\bar{Y}_j-\bar{Y}_k}{\text{MSE}\sqrt{\frac{1}{n_j}+\frac{1}{n_k}}} \sim t_{n-g}.\]

We zullen hier eerst demonstreren dat het werken met \(m\)-testen op het \(\alpha\) significantieniveau een foute aanpak is die de kans op een type I fout niet onder controle kan houden. Dit zal aanleiding geven tot een meer algemene definitie van de type I fout.

Alvorens de denkfout in de naïeve aanpak te demonsteren via simulaties, tonen we hoe de naïeve benadering in zijn werk zou gaan voor het prostacycline voorbeeld.

with(
  prostacyclin,
  pairwise.t.test(prostac, dose, "none")
  )

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  prostac and dose 
## 
##    10      25     
## 25 0.34927 -      
## 50 2e-05   0.00031
## 
## P value adjustment method: none

Deze output toont de tweezijdige \(p\)-waarden voor het testen van alle partiële hypotheses. We zouden hier kunnen besluiten dat het gemiddelde prostacycline niveau extreem significant verschillend is tussen de hoge en de lage dosis groep en tussen de hoge en de matige dosis groep (beide \(p<<0.001\)). Verder is het gemiddelde prostacycline niveau niet significant verschillend is tussen de matige en de lage dosis groep.

In onderstaande R code wordt een simulatiestudie opgezet (herhaalde steekproefname).

We simuleren uit een ANOVA model met \(g=3\) groepen.
De gemiddelden in het ANOVA model zijn gelijk aan elkaar, zodat de nulhypothese

\[H_0: \mu_1=\mu_2=\mu_3\]

opgaat.

Voor iedere gesimuleerde dataset zijn er \(m=3\) paarsgewijze two-sample \(t\)-testen
Zodra minstens één van de \(p\)-waarden kleiner is dan het significantieniveau \(\alpha=5\%\), wordt de nulhypothese \(H_0: \mu_1=\mu_2=\mu_3\) verworpen omdat er minstens twee gemiddelden verschillend zijn volgens de \(t\)-testen.
We rapporteren de relatieve frequentie van het verwerpen van de globale nulhypothese, meer bepaald de kans op een type I fout van de test voor \(H_0: \mu_1=\mu_2=\mu_3\).

g <- 3 # aantal behandelingen (g=3)
ni <- 12 # aantal herhalingen in iedere groep
n <- g*ni # totaal aantal observaties
alpha <- 0.05 # significantieniveau van een individuele test
N <- 10000 #aantal simulaties
set.seed(302) #seed zodat resultaten exact geproduceerd kunnen worden
trt <- factor(rep(1:g, ni)) #factor
cnt <- 0 #teller voor aantal foutieve verwerpingen

for(i in 1:N) {
if (i %% 1000 == 0) cat(i, "/", N, "\n")
y <- rnorm(n)
tests <- pairwise.t.test(y, trt, "none")
verwerp <- min(tests$p.value, na.rm = T) < alpha
if(verwerp) cnt <- cnt+1
}

## 1000 / 10000 
## 2000 / 10000 
## 3000 / 10000 
## 4000 / 10000 
## 5000 / 10000 
## 6000 / 10000 
## 7000 / 10000 
## 8000 / 10000 
## 9000 / 10000 
## 10000 / 10000

cnt/N

## [1] 0.1209

De simulatiestudie toont aan dat de kans op een type I fout gelijk is aan 12.1%, wat meer dan dubbel zo groot is dan de vooropgestelde \(\alpha=5\%\). Als we de simulatiestudie herhalen met \(g=5\) groepen (i.e. m=g(g-1)/2=10 paarsgewijze \(t\)-testen) dan vinden we \(28.0\%\) in plaats van de gewenste \(5\%\). Deze simulaties illustreren het probleem van multipliciteit (Engels: multiplicity): de klassieke \(p\)-waarden mogen enkel met het significantieniveau \(\alpha\) vergeleken worden, indien het besluit op exact één \(p\)-waarde gebaseerd is. Hier wordt het finale besluit (aldanniet verwerpen van \(H_0: \mu_1=\cdots =\mu_g\)) gebaseerd op \(m=g\times(g-1)/2\) \(p\)-waarden, met \(g\) het aantal groepen.

In de volgende sectie breiden we het begrip van type I fout uit en introduceren we enkele oplossingen om met multipliciteit om te gaan.