¿Cuáles son las diferencias prácticas entre los procedimientos de tasa de falsos descubrimientos de Benjamini y Hochberg (1995) y de Benjamini y Yekutieli (2001)?

Question

Mi programa estadístico implementa los procedimientos de tasa de falsos descubrimientos (FDR) de Benjamini & Hochberg (1995) y Benjamini & Yekutieli (2001). He hecho todo lo posible por leer el último documento, pero es bastante denso desde el punto de vista matemático y no estoy razonablemente seguro de entender la diferencia entre los procedimientos. Puedo ver en el código subyacente de mi programa de estadística que, en efecto, son diferentes y que este último incluye una cantidad q a la que he visto referirse en relación con el FDR, pero que tampoco entiendo del todo.

¿Hay alguna razón para preferir el procedimiento de Benjamini y Hochberg (1995) frente al de Benjamini y Yekutieli (2001)? ¿Tienen supuestos diferentes? ¿Cuáles son las diferencias prácticas entre estos enfoques?

Benjamini, Y., y Hochberg, Y. (1995). Control de la tasa de falsos descubrimientos: un enfoque práctico y potente de las pruebas múltiples. Journal of the Royal Statistical Society Series B, 57, 289-300.

Benjamini, Y., y Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29, 1165-1188.

El documento de 1999 al que se hace referencia en los comentarios más abajo: Yekutieli, D., & Benjamini, Y. (1999). Procedimientos de control de la tasa de falsos descubrimientos basados en el remuestreo para pruebas estadísticas correlacionadas. Journal of Statistical Planning and Inference, 82(1), 171-196.

Answer 1

Benjamini y Hochberg (1995) introdujeron la tasa de falso descubrimiento. Benjamini y Yekutieli (2001) demostraron que el estimador es válido bajo algunas formas de dependencia. La dependencia puede surgir de la siguiente manera. Considere la variable continua utilizada en una prueba t y otra variable correlacionada con ella; por ejemplo, probar si el IMC difiere en dos grupos y si la circunferencia de la cintura difiere en estos dos grupos. Como estas variables están correlacionadas, los valores p resultantes también estarán correlacionados. Yekutieli y Benjamini (1999) desarrollaron otro procedimiento de control de la FDR, que puede utilizarse bajo dependencia general remuestreando la distribución nula. Como la comparación es con respecto a la distribución nula de permutación, a medida que aumenta el número total de verdaderos positivos, el método se vuelve más conservador. Resulta que BH 1995 también es conservador a medida que aumenta el número de verdaderos positivos. Para mejorar esto, Benjamini y Hochberg (2000) introdujeron el procedimiento FDR adaptativo. Esto requirió la estimación de un parámetro, la proporción nula, que también se utiliza en el estimador pFDR de Storey. Storey ofrece comparaciones y argumenta que su método es más potente y destaca la naturaleza conservadora del procedimiento de 1995. Storey también tiene resultados y simulaciones bajo dependencia.

Todas las pruebas anteriores son válidas con independencia. La cuestión es qué tipo de desviación de la independencia pueden afrontar estas estimaciones.

Mi opinión actual es que si no se esperan demasiados verdaderos positivos, el procedimiento BY (1999) es bueno porque incorpora las características de distribución y la dependencia. Sin embargo, no conozco una implementación. El método de Storey fue diseñado para muchos verdaderos positivos con cierta dependencia. BH 1995 ofrece una alternativa a la tasa de error por familias y sigue siendo conservadora.

Benjamini, Y y Y Hochberg. On the Adaptive Control of the False Discovery Rate in Multiple Testing with Independent Statistics. Journal of Educational and Behavioral Statistics, 2000.

Answer 2

P.adjust no es misciting para BY. La referencia es al teorema 1.3 (prueba en la sección 5 de la página 1182) del documento:

Benjamini, Y., y Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29, 1165-1188.

Como en este documento se discuten varios ajustes diferentes, la referencia en la página de ayuda (en el momento de escribir este artículo) para p.adjust() es algo oscura. El método está garantizado para controlar la FDR, en la tasa indicada, bajo la estructura de dependencia más general. Hay comentarios informativos en las diapositivas de Christopher Genovese en www.stat.cmu.edu/~genovese/talks/hannover1-04.pdf Obsérvese el comentario de la diapositiva 37, que se refiere al método del Teorema 1.3 en el documento BY 2001 [method='BY' con p.adjust()] que: "Desafortunadamente, esto es típicamente muy conservador, a veces incluso más que Bonferroni".

Ejemplo numérico: method='BY' vs method='BH'

A continuación se compara el método='BY' con el método='BH', utilizando la función p.adjust() de R, para los valores p de la columna 2 de la Tabla 2 del artículo de Benjamini y Hochberg (2000):

> p <-    c(0.85628,0.60282,0.44008,0.41998,0.3864,0.3689,0.31162,0.23522,0.20964,
0.19388,0.15872,0.14374,0.10026,0.08226,0.07912,0.0659,0.05802,0.05572,
0.0549,0.04678,0.0465,0.04104,0.02036,0.00964,0.00904,0.00748,0.00404,
0.00282,0.002,0.0018,2e-05,2e-05,2e-05,0)
> pmat <- rbind(p,p.adjust(p, method='BH'),p.adjust(p, method='BY'))
> rownames(pmat)<-c("pval","adj='BH","adj='BY'")
> round(pmat,4)

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] pval 0.8563 0.6028 0.4401 0.4200 0.3864 0.3689 0.3116 0.2352 0.2096 adj='BH 0.8563 0.6211 0.4676 0.4606 0.4379 0.4325 0.3784 0.2962 0.2741 adj='BY' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] pval 0.1939 0.1587 0.1437 0.1003 0.0823 0.0791 0.0659 0.0580 0.0557 adj='BH 0.2637 0.2249 0.2125 0.1549 0.1332 0.1332 0.1179 0.1096 0.1096 adj='BY' 1.0000 0.9260 0.8751 0.6381 0.5485 0.5485 0.4856 0.4513 0.4513 [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] pval 0.0549 0.0468 0.0465 0.0410 0.0204 0.0096 0.0090 0.0075 0.0040 adj='BH 0.1096 0.1060 0.1060 0.1060 0.0577 0.0298 0.0298 0.0283 0.0172 adj='BY' 0.4513 0.4367 0.4367 0.4367 0.2376 0.1227 0.1227 0.1164 0.0707 [,28] [,29] [,30] [,31] [,32] [,33] [,34] pval 0.0028 0.0020 0.0018 0e+00 0e+00 0e+00 0 adj='BH 0.0137 0.0113 0.0113 2e-04 2e-04 2e-04 0 adj='BY' 0.0564 0.0467 0.0467 7e-04 7e-04 7e-04 0

Nota: El multiplicador que relaciona los valores BY con los valores BH es $\sum_{i=1}^m (1/i)$ , donde $m$ es el número de valores p. Los multiplicadores son, por ejemplo, los valores m = 30, 34, 226, 1674, 12365:

> mult <- sapply(c(11, 30, 34, 226, 1674, 12365), function(i)sum(1/(1:i)))

setNames(mult, paste(c('m =',rep('',5)), c(11, 30, 34, 226, 1674, 12365)) m = 11 30 34 226 1674 12365 3.020 3.995 4.118 6.000 8.000 10.000

Comprueba que para el ejemplo anterior, donde $m$ =34, el multiplicador es 4,118