Family-wise error rate

Wanneer \(m>1\) toetsen worden aangewend om 1 beslissing te vormen, is het noodzakelijk te corrigeren voor het risico op vals positieve resultaten. Meeste procedures voor meervoudig toetsen gaan ervan uit dat alle \(m\) nulhypotheses waar zijn. Er wordt dan geprobeerd om het risico op minstens 1 vals positief resultaat te controleren op experimentgewijs significantieniveau \(\alpha_E\), typisch \(\alpha_E=0.05\). In de engelstalige literatuur wordt het experimentgewijs significantieniveau family-wise error rate (FWER) genoemd.

Bonferroni correctie

Bij het uitvoeren van \(m\) onafhankelijke toetsen met elk significantieniveau \(\alpha\), is

\[\begin{eqnarray*} \alpha_E&=&\text{P}[\text{minstens 1 Type I fout}]\\ &=&1-(1-\alpha)^m \leq m\alpha \end{eqnarray*}\]

Als we 5 toetsen uitvoeren op het 5% significantieniveau is FWER \(\approx 25\%\).
Door ze op het 1% significantieniveau uit te voeren, bekomen we FWER \(\approx 5\%\).

De Bonferroni correctie houdt de FWER begrensd op \(\alpha_E\) door \(\alpha=\alpha_E/m\) te kiezen voor het uitvoeren van de \(m\) paarsgewijze vergelijkingen. Als alternatieve methode kunnen we ook

aangepaste p-waarden rapporteren zodat we deze met het experimentgewijze \(\alpha_E\) niveau kunnen vergelijken:

\[\tilde{p}=min(m\times p,1)\]

\((1-\alpha_E/m)100\%\) betrouwbaarheidsintervallen rapporteren.

Het gebruik van aangepaste p-waarden heeft als voordeel dat de lezer deze zelf kan interpreteren en hij/zij een maat kan geven voor de significantie op het experimentsgewijs significantie niveau. Merk op dat we de aangepaste p-waarden begrenzen op \(\tilde{p}=1\) omdat p-waarden kansen zijn en steeds tussen 0 en 1 dienen te liggen.

Onderstaande R code geeft de resultaten (gecorrigeerde \(p\)-waarden) na correctie met de methode van Bonferroni.

with(
  prostacyclin,
  pairwise.t.test(
    prostac,
    dose,
    p.adjust.method="bonferroni")
  )

## 
## 	Pairwise comparisons using t tests with pooled SD 
## 
## data:  prostac and dose 
## 
##    10      25     
## 25 1.00000 -      
## 50 6e-05   0.00094
## 
## P value adjustment method: bonferroni

De conclusies blijven hetzelfde behalve dat de FWER nu gecontroleerd is \(\alpha=5\%\) en de \(\tilde{p}\)-waarden een factor 3 groter zijn.

Dezelfde analyse kan uitgevoerd worden met het multcomp R package dat speciaal werd ontwikkeld voor multipliciteit in lineaire modellen.

library(multcomp)
model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey"))
summary(model1.mcp, test = adjusted("bonferroni"))

## 
## 	 Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = prostac ~ dose, data = prostacyclin)
## 
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(>|t|)    
## 25 - 10 == 0    8.258      8.698   0.949 1.000000    
## 50 - 10 == 0   43.258      8.698   4.974 5.98e-05 ***
## 50 - 25 == 0   35.000      8.698   4.024 0.000943 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Om Bonferroni aangepaste betrouwbaarheidsintervallen te verkrijgen moeten we eerst zelf functie definiëren in R om bonferroni kritische waarde te bepalen. We noemen deze functie calpha_bon_t.

calpha_bon_t <- function(object, level)
{
    abs(
      qt(
        (1-level)/2/nrow(object$linfct),
        object$df)
      )
}
confint(model1.mcp, calpha = calpha_bon_t)

## 
## 	 Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = prostac ~ dose, data = prostacyclin)
## 
## Quantile = 2.5222
## 95% confidence level
##  
## 
## Linear Hypotheses:
##              Estimate lwr      upr     
## 25 - 10 == 0   8.2583 -13.6790  30.1957
## 50 - 10 == 0  43.2583  21.3210  65.1957
## 50 - 25 == 0  35.0000  13.0626  56.9374

We zullen nu het effect van de Bonferroni correctie opnieuw nagaan via simulaties.

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, "bonferroni")
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.0457

We vinden dus een FWER van 4.6% (een beetje conservatief). Wanneer we de simulaties doen voor \(g=5\) groepen, vinden we een FWER van \(4.1\%\) (conservatiever). Door de Bonferroni correctie is de kans op minstens één vals positief resultaat \(< \alpha_E\). Hoewel de FWER wordt gecontroleerd door de Bonferroni methode, kan een verlies aan power worden verwacht aangezien het werkelijke niveau lager is dan het vooropgestelde 5% experimentsgewijs significantieniveau.