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[48]. 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.