El problema de las comparaciones múltiples

En los contrastes de hipótesis se controlan habitualmente dos posibles tipos de error:

En la realización de un test de hipótesis se define siempre un estadístico de contraste, que es una función que mide la discrepancia entre lo que muestran los datos y lo que se esperaría observar si la hipótesis nula fuera cierta. Si, una vez tomados los datos, y supuesto que la hipótesis nula es cierta, resulta poco probable observar una discrepancia tan grande o mayor que la observada, entonces la hipótesis nula se rechaza; en caso contrario se acepta. Se considera poco probable si dicha probabilidad es menor que el nivel de significación \(\alpha\) del test, en muchos casos fijado en el valor 0.05.

Ejemplo.

Se desea contrastar si una moneda está trucada. Desde la perspectiva del contraste de hipótesis, se debe elegir entre las dos siguientes opciones:

\[\begin{cases} H_{0}: & \textrm{La moneda está equilibrada}\\ H_{1}: & \textrm{La moneda está trucada} \end{cases}\]

Para ello se realiza el experimento de lanzar la moneda 10 veces. El estadístico de contraste en este caso consiste simplemente en contar el número de caras. Si este número es menor que 2 o mayor que 8 la moneda se declara trucada. El número de caras que ocurren al lanzar 10 veces una moneda equilibrada es una variable aleatoria con distribución binomial de parámetros \(p=0.5\) y \(n=10\). La probabilidad de que una moneda equilibrada sea declarada como trucada es entonces (siendo \(X\) el número de caras):

\[P\left(X<2\right)+P\left(X>8\right)=P\left(X=0\right)+P\left(X=1\right)+P\left(X=9\right)+P\left(X=10\right)\] que en R puede calcularse mediante:

sum(dbinom(c(0,1,9,10),10,0.6))
## [1] 0.04803512

Por tanto, de acuerdo con el criterio elegido, la probabilidad de decidir que la moneda está trucada (hipótesis alternativa) cuando en realidad es equilibrada es 0.0480351, algo menos del 5%.

Pero ¿que ocurre si se nos presentan 10 monedas y tenemos que decidir si hay alguna trucada?. Aplicando el criterio anterior, si las 10 monedas estuviesen equilibradas, la probabilidad de que las 10 fuesen declaradas no trucadas seria:

\[\left(1-0.0480351\right)^{10}=0.9519649^{10}= 0.6112365\]

Por tanto, la probabilidad de que al menos una sea declarada trucada es

\[1-0.6112365=0.3887635\]

Por tanto, ¡aunque las 10 monedas sean correctas, la probabilidad de que al menos una sea declarada como trucada es casi del 40%! Por tanto corremos un riesgo de casi un 40% de cometer un error tipo I (rechazar la hipótesis nula cuando es cierta), muy lejos del 5% planificado inicialmente.

 

 

Procedimientos para comparaciones múltiples

Los tests stepwise son más potentes que los de paso único (rechazan más hipótesis nulas falsas; toda hipótesis rechazada por un stepwise es también rechazada por un single-step, pero no a la inversa). Sin embargo los intervalos de confianza para los parámetros de interés son más difíciles de construir en los tests stepwise.

 

 

Tests de Bonferroni

Si se realizan \(m\) contrastes, el procedimiento de Bonferroni consiste simplemente en comparar el p-valor (sin ajustar) de cada test con el valor \(\alpha/m\), siendo \(\alpha\) el nivel de significación global para todos los tests, esto es, la probabilidad de rechazar al menos una hipótesis nula bajo la hipótesis de que todas las hipótesis nulas son ciertas. Una hipótesis nula se rechaza si el p-valor del test asociado es menor que \(\alpha/m\).

De manera equivalente, se puede definir el p-valor ajustado del test \(i\) como \(m\cdot p_i\) siendo \(p_i\) el p-valor sin ajustar; la hipótesis nula \(i\) se rechaza entonces si su p-valor ajustado es menor que \(\alpha\). En R, si para varios contrastes se han obtenido los siguientes p-valores:

p <- c( 0.01, 0.015, 0.02, 0.04)

Los p-valores ajustados por Bonferroni se obtienen como:

p.adjust(p,"bonferroni")
## [1] 0.04 0.06 0.08 0.16

Como vemos, aunque con los p-valores “crudos”, se rechazarían las cuatro hipótesis nulas, tras el ajuste por Bonferroni, la única que se rechazaría es la primera.

En cualquier caso, el test de Bonferroni es conservador: tiene poca potencia.

 

Procedimiento de Holm

Partiendo del método de Bonferroni, el método de Holm es un método paso a paso que mejora la potencia de aquél. El procedimiento es el siguiente:

  1. Ordenar los p-valores de menor a mayor: \(p_1 \le p_2 \le \dots \le p_m\)

  2. Si \(p_1 > \alpha/m\) aceptar todas las hipótesis nulas (pues si el menor p-valor es mayor que \(\alpha/m\) también lo serán todos los demás)

  3. En caso contrario rechazar la hipótesis nula asociada al p-valor \(p_1\). Quedan ahora \(m-1\) hipótesis por contrastar. Si \(p_2 > \alpha/(m-1)\) aceptar las hipótesis nulas restantes. Si no, rechazar la hipótesis nula asociada al p-valor \(p_2\).

  4. Repetir, disminuyendo el valor del divisor una unidad en cada paso.

p.adjust(p,"holm")
## [1] 0.040 0.045 0.045 0.045

 

Procedimientos para comparaciones múltiples en modelos lineales.

En el caso de los contrastes asociados a los parámetros en modelos lineales, se puede conseguir mayor potencia en los contrastes si se tiene en cuenta la estructura de correlación entre las estimaciones de dichos parámetros.

Análisis de la varianza

# ANOVA DATA
sd=1
dd <- data.frame(grupo = gl(4,5)) # balanced 2-way
M=model.matrix(~ grupo , dd)
beta=c(8,0,0,0)
dd$x=as.vector(as.matrix(M)%*%beta+rnorm(nrow(M),0,sd))
dd
##    grupo        x
## 1      1 9.037792
## 2      1 7.861754
## 3      1 8.026735
## 4      1 8.610903
## 5      1 7.370009
## 6      2 9.558112
## 7      2 7.899853
## 8      2 9.575086
## 9      2 7.550465
## 10     2 7.656223
## 11     3 8.906318
## 12     3 7.368612
## 13     3 7.389259
## 14     3 8.868736
## 15     3 7.175251
## 16     4 7.487441
## 17     4 8.277980
## 18     4 8.234794
## 19     4 8.240724
## 20     4 7.537432
aggregate(x~grupo,dd,function(x) c(media=mean(x),sd=sd(x)))
##   grupo   x.media      x.sd
## 1     1 8.1814386 0.6523669
## 2     2 8.4479477 1.0290309
## 3     3 7.9416354 0.8676041
## 4     4 7.9556742 0.4053426
with(dd,boxplot(x~grupo,col="lightblue3"))

library(car)
leveneTest(x~grupo,data=dd)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.5543 0.6526
##       16
adeva=aov(x~grupo,data=dd)
shapiro.test(residuals(adeva))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(adeva)
## W = 0.89035, p-value = 0.02729
qqPlot(residuals(adeva))

summary(adeva)
##             Df Sum Sq Mean Sq F value Pr(>F)
## grupo        3  0.848  0.2827   0.471  0.707
## Residuals   16  9.606  0.6004
TukeyHSD(adeva,"grupo")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = x ~ grupo, data = dd)
## 
## $grupo
##            diff       lwr       upr     p adj
## 2-1  0.26650909 -1.135544 1.6685625 0.9468721
## 3-1 -0.23980326 -1.641857 1.1622502 0.9603532
## 4-1 -0.22576440 -1.627818 1.1762891 0.9665269
## 3-2 -0.50631234 -1.908366 0.8957411 0.7329864
## 4-2 -0.49227348 -1.894327 0.9097800 0.7490527
## 4-3  0.01403886 -1.388015 1.4160923 0.9999912
# Con la librería multcomp
library(multcomp)
summary(glht(adeva,linfct=mcp(grupo="Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = x ~ grupo, data = dd)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0  0.26651    0.49005   0.544    0.947
## 3 - 1 == 0 -0.23980    0.49005  -0.489    0.960
## 4 - 1 == 0 -0.22576    0.49005  -0.461    0.967
## 3 - 2 == 0 -0.50631    0.49005  -1.033    0.733
## 4 - 2 == 0 -0.49227    0.49005  -1.005    0.749
## 4 - 3 == 0  0.01404    0.49005   0.029    1.000
## (Adjusted p values reported -- single-step method)
confint(glht(adeva,linfct=mcp(grupo="Tukey")))
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = x ~ grupo, data = dd)
## 
## Quantile = 2.8593
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##            Estimate lwr      upr     
## 2 - 1 == 0  0.26651 -1.13471  1.66772
## 3 - 1 == 0 -0.23980 -1.64102  1.16141
## 4 - 1 == 0 -0.22576 -1.62698  1.17545
## 3 - 2 == 0 -0.50631 -1.90753  0.89490
## 4 - 2 == 0 -0.49227 -1.89349  0.90894
## 4 - 3 == 0  0.01404 -1.38718  1.41525
groupSizes=table(dd$grupo)
contrastes=contrMat(groupSizes,type="Dunnet",base=3)
summary(glht(adeva,linfct=mcp(grupo=contrastes)))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = x ~ grupo, data = dd)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)
## 1 - 3 == 0  0.23980    0.49005   0.489    0.928
## 2 - 3 == 0  0.50631    0.49005   1.033    0.613
## 4 - 3 == 0  0.01404    0.49005   0.029    1.000
## (Adjusted p values reported -- single-step method)
confint(glht(adeva,linfct=mcp(grupo=contrastes)))
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = x ~ grupo, data = dd)
## 
## Quantile = 2.5928
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##            Estimate lwr      upr     
## 1 - 3 == 0  0.23980 -1.03083  1.51043
## 2 - 3 == 0  0.50631 -0.76432  1.77694
## 4 - 3 == 0  0.01404 -1.25659  1.28467

Modelo de regresión lineal

library(ISwR)
## 
## Attaching package: 'ISwR'
## The following object is masked from 'package:survival':
## 
##     lung
data("thuesen")
plot(thuesen)
thuesen.lm=lm(short.velocity~blood.glucose,thuesen)
summary(thuesen.lm)
## 
## Call:
## lm(formula = short.velocity ~ blood.glucose, data = thuesen)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.40141 -0.14760 -0.02202  0.03001  0.43490 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.09781    0.11748   9.345 6.26e-09 ***
## blood.glucose  0.02196    0.01045   2.101   0.0479 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2167 on 21 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1737, Adjusted R-squared:  0.1343 
## F-statistic: 4.414 on 1 and 21 DF,  p-value: 0.0479
abline(thuesen.lm)

Ajuste de Bonferroni (Procedimiento más malo)

library("multcomp")
thuesen.mc <- glht(thuesen.lm, linfct = diag(2))
summary(thuesen.mc, test = adjusted(type = "bonferroni"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = short.velocity ~ blood.glucose, data = thuesen)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  1.09781    0.11748   9.345 1.25e-08 ***
## 2 == 0  0.02196    0.01045   2.101   0.0958 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Método de paso único:

Se utiliza la distribución T de Hotelling (equivalente a la t de Student unidimensional); tiene en cuenta que los coeficientes del modelo pueden estar correlados

library("multcomp")
thuesen.mc <- glht(thuesen.lm, linfct = diag(2))
summary(thuesen.mc)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = short.velocity ~ blood.glucose, data = thuesen)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  1.09781    0.11748   9.345    1e-08 ***
## 2 == 0  0.02196    0.01045   2.101   0.0645 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(thuesen.mc)
## 
##   Simultaneous Confidence Intervals
## 
## Fit: lm(formula = short.velocity ~ blood.glucose, data = thuesen)
## 
## Quantile = 2.23
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##        Estimate  lwr       upr      
## 1 == 0  1.097815  0.835837  1.359793
## 2 == 0  0.021963 -0.001348  0.045274

Método paso a paso (método de Westfall)

Este método es más potente

summary(thuesen.mc, test = adjusted(type = "Westfall"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = short.velocity ~ blood.glucose, data = thuesen)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  1.09781    0.11748   9.345    1e-08 ***
## 2 == 0  0.02196    0.01045   2.101   0.0479 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- Westfall method)