14.6: Comparaciones múltiples y pruebas post hoc

Última actualización
Guardar como PDF

Page ID: 151304

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

Cada vez que ejecutas un ANOVA con más de dos grupos, y terminas con un efecto significativo, lo primero que probablemente querrás preguntar es qué grupos son realmente diferentes entre sí. En nuestro ejemplo de drogas, nuestra hipótesis nula fue que los tres medicamentos (placebo, Anxifree y Joyzepam) tienen exactamente el mismo efecto en el estado de ánimo. Pero si lo piensas, la hipótesis nula en realidad es reclamar tres cosas diferentes todas a la vez aquí. En concreto, afirma que:

El medicamento de su competidor (Anxifree) no es mejor que un placebo (es decir, μ _A = μ _P)
Su medicamento (Joyzepam) no es mejor que un placebo (es decir, μ _J = μ _P)
Anxifree y Joyzepam son igualmente efectivos (es decir, μ _J = μ _A)

Si alguna de esas tres afirmaciones es falsa, entonces la hipótesis nula también es falsa. Entonces, ahora que hemos rechazado nuestra hipótesis nula, estamos pensando que al menos una de esas cosas no es verdad. Pero, ¿cuáles? Estas tres proposiciones son de interés: ciertamente quieres saber si tu nuevo medicamento Joyzepam es mejor que un placebo, y sería bueno saber qué tan bien se compara con una alternativa comercial existente (es decir, Anxifree). Incluso sería útil verificar el desempeño de Anxifree frente al placebo: incluso si Anxifree ya ha sido ampliamente probado contra placebos por otros investigadores, aún puede ser muy útil verificar que su estudio esté produciendo resultados similares a trabajos anteriores.

Cuando caracterizamos la hipótesis nula en términos de estas tres proposiciones distintas, queda claro que hay ocho posibles “estados del mundo” que necesitamos distinguir entre:

posibilidad:	es μ _P =μ _A?	es μ _P =μ _J?	es μ _A =μ _J?	¿cuál hipótesis?
1	✓	✓	✓	null
2	✓	✓		alternativa
3	✓		✓	alternativa
4	✓			alternativa
5		✓	✓	alternativa
6		✓		alternativa
7			✓	alternativa
8				alternativa

Al rechazar la hipótesis nula, hemos decidido que no creemos que #1 sea el verdadero estado del mundo. La siguiente pregunta es, ¿cuál de las otras siete posibilidades creemos que es la correcta? Ante esta situación, suele ayudar a mirar los datos. Por ejemplo, si miramos las parcelas de la Figura 14.1, es tentador concluir que Joyzepam es mejor que el placebo y mejor que Anxifree, pero no hay diferencia real entre Anxifree y el placebo. No obstante, si queremos obtener una respuesta más clara al respecto, podría ser útil ejecutar algunas pruebas.

Ejecución de pruebas t “por pares”

¿Cómo podríamos resolver nuestro problema? Dado que tenemos tres pares separados de medias (placebo versus Anxifree, placebo versus Joyzepam y Anxifree versus Joyzepam) para comparar, lo que podríamos hacer es ejecutar tres pruebas t separadas y ver qué pasa. Hay un par de formas en las que podríamos hacer esto. Un método sería construir nuevas variables correspondientes a los grupos que desea comparar (por ejemplo, anxifree, placebo y joyzepam), y luego ejecutar una prueba t sobre estas nuevas variables:

t.test( anxifree, placebo, var.equal = TRUE )   # Student t-test

anxifree <- with(clin.trial, mood.gain[drug == "anxifree"])  # mood change due to anxifree
placebo <- with(clin.trial, mood.gain[drug == "placebo"])    # mood change due to placebo

o bien, podrías usar el argumento subconjunto en la función t.test () para seleccionar solo aquellas observaciones correspondientes a uno de los dos grupos que nos interesan:

t.test( formula = mood.gain ~ drug, 
       data = clin.trial, 
       subset = drug %in% c("placebo","anxifree"), 
       var.equal = TRUE 
)

Consulta el Capítulo 7 si has olvidado cómo funciona el operador %en%. Independientemente de la versión que hagamos, R imprimirá los resultados de la prueba t, aunque no he incluido esa salida aquí. Si continuamos haciendo esto para todos los pares de variables posibles, podemos mirar para ver cuáles (si los hay) pares de grupos son significativamente diferentes entre sí. Esta “idea de muchas pruebas t” no es una mala estrategia, aunque como veremos más adelante hay algunos problemas con ella. Sin embargo, por el momento nuestro mayor problema es que es un dolor tener que escribir un comando tan largo una y otra vez: por ejemplo, si tu experimento tiene 10 grupos, entonces tienes que ejecutar 45 pruebas t. Eso es demasiado mecanografiar.

Para ayudar a mantener la escritura al mínimo, R proporciona una función llamada pairwise.t.test () que ejecuta automáticamente todas las pruebas t por usted. Hay tres argumentos que necesitas especificar, la variable resultado x, la variable de grupo g y el argumento p.adjust.method, que “ajusta” el valor p de una forma u otra. Voy a explicar el ajuste del valor p en un momento, pero por ahora solo podemos establecer p.adjust.method = “none” ya que no estamos haciendo ningún ajuste. Para nuestro ejemplo, esto es lo que hacemos:

pairwise.t.test( x = clin.trial$mood.gain,   # outcome variable
                  g = clin.trial$drug,        # grouping variable
                  p.adjust.method = "none"    # which correction to use?
 )

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  clin.trial$mood.gain and clin.trial$drug 
## 
##          placebo anxifree
## anxifree 0.15021 -       
## joyzepam 3e-05   0.00056 
## 
## P value adjustment method: none

Una cosa que me fastidia un poco sobre la función pairwise.t.test () es que no puedes simplemente darle un objeto aov, y hacer que produzca esta salida. Después de todo, me metí en todos esos problemas antes de conseguir que R creara la variable my.anova y —como vimos en la Sección 14.3.2— R en realidad ha almacenado suficiente información dentro de ella como para que solo pueda hacer que ejecute todas las pruebas por pares usando my.anova como entrada. Para ello, he incluido una función posthocpairWiset () en el paquete lsr que te permite hacer esto. La idea detrás de esta función es que solo puedes ingresar el propio objeto aov, ²⁰⁸ y luego obtener las pruebas por pares como salida. A partir de la escritura actual, posthocPairWiset () es en realidad solo una forma sencilla de llamar a la función pairwise.t.test (), pero debes tener en cuenta que tengo la intención de hacer algunos cambios más adelante. Aquí hay un ejemplo:

posthocPairwiseT( x = my.anova, p.adjust.method = "none" )

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  mood.gain and drug 
## 
##          placebo anxifree
## anxifree 0.15021 -       
## joyzepam 3e-05   0.00056 
## 
## P value adjustment method: none

En versiones posteriores, planeo agregar más funcionalidad (por ejemplo, intervalos de confianza ajustados), pero por ahora creo que es al menos algo útil. Para ver por qué, supongamos que has ejecutado tu ANOVA y has almacenado los resultados en my.anova, y estás contento usando la corrección Holm (el método predeterminado en pairwise.t.test (), que voy a explicar esto en un momento). En ese caso, todo lo que tienes que hacer es escribir esto:

posthocPairwiseT( my.anova )

y R dará salida a los resultados de la prueba. Mucho más conveniente, creo.

Correcciones para pruebas múltiples

En la sección anterior insinué que hay un problema con solo correr lotes y muchas pruebas t. La preocupación es que al ejecutar estos análisis, lo que estamos haciendo es ir a una “expedición pesquera”: estamos haciendo montones y montones de pruebas sin mucha orientación teórica, con la esperanza de que algunas de ellas salgan significativas. Este tipo de búsqueda libre de teoría de diferencias de grupo se conoce como análisis post hoc (“post hoc” es latino para “después de esto”). ²⁰⁹

Está bien realizar análisis post hoc, pero se requiere mucho cuidado. Por ejemplo, el análisis que ejecuté en la sección anterior es en realidad bastante peligroso: cada prueba t individual está diseñada para tener una tasa de error de tipo I del 5% (es decir, α=.05), y ejecuté tres de estas pruebas. Imagínese lo que hubiera pasado si mi ANOVA involucrara a 10 grupos diferentes, y hubiera decidido ejecutar 45 pruebas t “post hoc” para tratar de averiguar cuáles eran significativamente diferentes entre sí, se esperaría que 2 o 3 de ellas surgieran significativas solo por casualidad. Como vimos en el Capítulo 11, el principio central de organización detrás de las pruebas de hipótesis nulas es que buscamos controlar nuestra tasa de errores Tipo I, pero ahora que estoy ejecutando muchas pruebas t a la vez, para determinar la fuente de mis resultados ANOVA, mi tasa de error Tipo I real en todo este conjunto familia de pruebas se ha salido completamente de control.

La solución habitual a este problema es introducir un ajuste al valor p, que tiene como objetivo controlar la tasa de error total en toda la familia de pruebas (ver Shaffer 1995). Un ajuste de esta forma, que suele aplicarse (pero no siempre) porque se está haciendo análisis post hoc, a menudo se le denomina corrección para comparaciones múltiples, aunque a veces se le denomina “inferencia simultánea”. En cualquier caso, hay bastantes formas diferentes de hacer este ajuste. Voy a discutir algunos de ellos en esta sección y en la Sección 16.8, pero debe tener en cuenta que hay muchos otros métodos por ahí (ver, por ejemplo, Hsu 1996).

Correcciones Bonferroni

El más simple de estos ajustes se llama la corrección Bonferroni (Dunn 1961), y en verdad es muy, muy simple. Supongamos que mi análisis post hoc consiste en m pruebas separadas, y quiero asegurar que la probabilidad total de cometer algún error de Tipo I es como máximo α. ²¹⁰ Si es así, entonces la corrección de Bonferroni solo dice “multiplica todos tus valores p brutos por m”. Si dejamos que p denote el valor p original, y dejamos que p′ _j sea el valor corregido, entonces la corrección de Bonferroni dice que:

p′=m×p

Y por lo tanto, si estás usando la corrección Bonferroni, rechazarías la hipótesis nula si p′<α. La lógica detrás de esta corrección es muy sencilla. Estamos haciendo m pruebas diferentes; así que si lo organizamos para que cada prueba tenga una tasa de error Tipo I de como máximo α/m, entonces la tasa de error total de Tipo I en estas pruebas no puede ser mayor que α. Eso es bastante simple, tanto es así que en el artículo original, el autor escribe:

El método dado aquí es tan sencillo y tan general que estoy seguro de que debió haber sido usado antes de esto. No lo encuentro, sin embargo, así que sólo puedo concluir que quizás su misma sencillez ha impedido a los estadísticos darse cuenta de que es un método muy bueno en algunas situaciones (pp 52-53 Dunn 1961)

Para usar la corrección Bonferroni en R, puede usar la función parwise.t.test (), ²¹¹ asegurándose de que establece p.adjust.method = “bonferroni”. Alternativamente, ya que toda la razón por la que estamos haciendo estas pruebas por pares en primer lugar es porque tenemos un ANOVA que estamos tratando de entender, probablemente sea más conveniente usar la función posthocpairWiset () en el paquete lsr, ya que podemos usar my.anova como la entrada:

posthocPairwiseT( my.anova, p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  mood.gain and drug 
## 
##          placebo anxifree
## anxifree 0.4506  -       
## joyzepam 9.1e-05 0.0017  
## 
## P value adjustment method: bonferroni

Si comparamos estos tres valores p con los que vimos en la sección anterior cuando no hicimos ningún ajuste en absoluto, es claro que lo único que ha hecho R es multiplicarlos por 3.

Correcciones Holm

Si bien la corrección de Bonferroni es el ajuste más simple que existe, no suele ser el mejor para usar. Un método que a menudo se usa en su lugar es la corrección Holm (Holm 1979). La idea detrás de la corrección Holm es fingir que estás haciendo las pruebas secuencialmente; comenzando con el valor p más pequeño (crudo) y pasando al más grande. Para el j-ésimo mayor de los valores p, el ajuste es

p′ _j = j×p _j

(es decir, el valor p más grande permanece sin cambios, el segundo valor p más grande se duplica, el tercer valor p más grande se triplica, y así sucesivamente), o

p′ _j =p′ _j+1

el que sea más grande. Esto puede sonar un poco confuso, así que vamos a repasarlo un poco más despacio. Esto es lo que hace la corrección Holm. Primero, ordenas todos tus valores p en orden, desde el más pequeño hasta el más grande. Para el valor p más pequeño todo lo que haces es multiplicarlo por m, y ya terminaste. Sin embargo, para todos los demás es un proceso de dos etapas. Por ejemplo, cuando se mueve al segundo valor de p más pequeño, primero lo multiplica por m−1. Si esto produce un número que es mayor que el valor p ajustado que obtuviste la última vez, entonces lo mantienes. Pero si es más pequeño que el último, entonces se copia el último valor p. Para ilustrar cómo funciona esto, considere la siguiente tabla, que muestra los cálculos de una corrección Holm para una colección de cinco valores p:

crudo p	rango j	p×j	Holm p
.001	5	.005	.005
.005	4	.020	.020
.019	3	.057	.057
.022	2	.044	.057
.103	1	.103	.103

Ojalá eso deje las cosas claras.

Aunque es un poco más difícil de calcular, la corrección Holm tiene algunas propiedades muy bonitas: es más potente que Bonferroni (es decir, tiene una tasa de error tipo II más baja), pero —contradictorio por parecer— tiene la misma tasa de error Tipo I. Como consecuencia, en la práctica nunca hay razón para utilizar la corrección más simple de Bonferroni, ya que siempre es superada por la corrección Holm un poco más elaborada. Debido a esto, la corrección Holm es la predeterminada utilizada por pairwise.t.test () y posthocPairWiset (). Para ejecutar la corrección Holm en R, podrías especificar p.adjust.method = “Holm” si quisieras, pero como es el valor predeterminado solo puedes hacer esto:

posthocPairwiseT( my.anova )

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  mood.gain and drug 
## 
##          placebo anxifree
## anxifree 0.1502  -       
## joyzepam 9.1e-05 0.0011  
## 
## P value adjustment method: holm

Como puedes ver, el valor p más grande (correspondiente a la comparación entre Anxifree y el placebo) está inalterado: a un valor de .15, es exactamente el mismo que el valor que obtuvimos originalmente cuando no aplicamos ninguna corrección. En contraste, el valor p más pequeño (Joyzepam versus placebo) se ha multiplicado por tres.

Redactar la prueba post hoc

Finalmente, después de haber ejecutado el análisis post hoc para determinar qué grupos son significativamente diferentes entre sí, podrías escribir el resultado así:

Las pruebas post hoc (usando la corrección Holm para ajustar p) indicaron que Joyzepam produjo un cambio de humor significativamente mayor que Anxifree (p=.001) y el placebo (p=9.1×10 ⁻⁵). No encontramos evidencia de que Anxifree se desempeñara mejor que el placebo (p=.15).

O bien, si no te gusta la idea de reportar valores p exactos, entonces cambiarías esos números a p<.01, p<.001 y p>.05 respectivamente. De cualquier manera, lo clave es que indiques que usaste la corrección de Holm para ajustar los valores p. Y claro, estoy asumiendo que en otra parte de la redacción has incluido las estadísticas descriptivas relevantes (es decir, las medias grupales y las desviaciones estándar), ya que estos valores p por sí solos no son terriblemente informativos.