Methode van Tukey

De methode van Tukey is een minder conservatieve methode voor het uitvoeren van post hoc testen. De implementatie benadert de nuldistributie van de posthoc test d.m.v. simulaties. De resultaten kunnen daarom lichtjes verschillen wanneer je de posthoc analyse opnieuw uitvoert.
De details van de methode vallen buiten het bestek van deze cursus. Via de implementatie in het multcomp package kunnen we opnieuw aangepaste p-waarden verkrijgen en aangepaste betrouwbaarheidsintervallen voor alle \(m\) paarsgewijze testen. We hoeven zelf geen functies te definiëren voor het verkrijgen van Tukey gecorrigeerde BIs.

model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey"))
summary(model1.mcp)

## 
##   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 0.613390    
## 50 - 10 == 0   43.258      8.698   4.974  < 1e-04 ***
## 50 - 25 == 0   35.000      8.698   4.024 0.000835 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

confint(model1.mcp)

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = prostac ~ dose, data = prostacyclin)
## 
## Quantile = 2.4539
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##              Estimate lwr      upr     
## 25 - 10 == 0   8.2583 -13.0849  29.6016
## 50 - 10 == 0  43.2583  21.9151  64.6016
## 50 - 25 == 0  35.0000  13.6567  56.3433

Merk op dat de Tukey methode smallere BIs en kleinere aangepaste p-waarden teruggeeft dan Bonferroni en dus minder conservatief is. De betrouwbaarheidsintervallen kunnen ook grafisch worden weergegeven wat handig is als er veel vergelijkingen worden uitgevoerd (zie Figuur 48).

model1.mcp %>%
  confint %>%
  plot

Figuur 48: $$95\%$$ experimentsgewijze betrouwbaarheidsintervallen voor de paarsgewijze verschillen in gemiddeld prostacycline niveau tussen alle arachidonzuur dosisgroepen. De BIs zijn gecorrigeerd voor multipliciteit via de Tukey methode.

Hierop zien we onmiddellijk dat het effect van de hoogste dosisgroep verschillend is van de laagste en middelste dosisgroep en dat er geen significant verschil is tussen de laagste en de middelste dosisgroep op het 5% experimentsgewijze significantieniveau.

Tenslotte gaan we ook via simulatie na of de Tukey methode de FWER correct kan controleren.

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)
m <- lm(y ~ trt)
m.mcp <- glht(m, linfct = mcp(trt = "Tukey"))
tests <- summary(m.mcp)$test
verwerp <- min(
  as.numeric(tests$pvalues),
  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.0503

We vinden dus een FWER van \(5.03\%\) wat heel dicht hij het nominale FWER\(=5\%\) ligt. Voor \(g=5\) groepen, vinden we een FWER van \(5.2\%\), wat ook vrij goed is.[49]