Wanneer
Bij het uitvoeren van
De Bonferroni correctie houdt de FWER begrensd op
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
Onderstaande R code geeft de resultaten (gecorrigeerde
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
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