1 El problema de hipotesis multiple Edgar Acuna Universidad de Puerto Rico en Mayaguez

2 EL problema de pruebas de hipotesis multiples Prueba simultanea de m hipotesis nulas, una por cada gene j (j=1,…m) H j : No hay relacion entre el nivel de expresion del gene j y las distintas condiciones (gene no expresados) Debido a que los experimentos con microarreglos monitorean simultaneamente niveles de expresion de miles de genes hay que hacer ajuste a los p-values de las hipotesis individuales.

3 El p-value en hipotesis simples Resultado de la prueba (p>  ): No se rechaza Ho Resultado de la prueba (p<  ). Se rechaza Ho La Ho realmente es verdadera Especificidad Error Tipo I(  ) (probabilidad  ) La Ho no es verdadera Error tipo II (probabilidad  ) Sensitividad El p-value (nivel de significacion observado) controla la tasa de error tipo I. La probabilidad de rechazar una hipotesis por error no es mayor que  (=.05). De 100 pruebas solo 5 serian significativas

4 Posibilidades en Hipotesis multiples # no rechazadas # rechazadas total # H ciertasUV (F +)m0m0 # H no ciertas T(F-)Sm1m1 totalm - RRm Verdad Decision m es fijo, m o y m 1 son fijas pero desconocidas, R es aleatorio y observado, V es aleatorio y no observado.

5 El problema de multiplicidad Cuando miles de hipotesis son probadas simultaneamente se incrementa la posibilidad de los falsos positivos Por ejemplo, en un microarray con 10,000 genes al 5% de significacion se puede esperar que 500 de ellos sean identificados como diferencialmente expresados ( es decir sus p-values serian menores que.05) simplemente por el azar. Se ha perdido control del error tipo I y los p-values individuales no se pueden usar directamente para concluir que el gen es expresado. Necesitan ser ajustados (corregidos)

6 Tasas de error tipo I (Falsos Positivos) a controlar Tasa de error por familia PFER = E(V): numero esperado de falsos positivos Tasa de error por comparacion PCER = E(V)/m: proporcion esperada de falsos positivos Family-wise Error Rate FWER = p(V ≥ 1) (probabilidad de al menos un falso positivo) False Discovery Rate: FDR = E(V/R)P(R>0) Proporcion esperada de falsos positivos entre las pruebas que fueron significativas

7 Comparacion de tasas de error Tipo I En general, dado un procedimiento de prueba de hipotesis multiple PCER  FWER  PFER, y FDR  FWER,

8 Tipos de control del error en hipotesis multiples Control debil. Se controla el error tipo I asumiendo que todas las hipotesis nulas ( ) son ciertas. Es decir, m o =m. No bria genes diferncialemente expresados. Control fuerte. Se controla el error tipo I asumiendo cualquier combinacion de hipotesis nulas y falsas. Algunos genes serian diferencialmente expresados

9 p-values ajustados (p*) Objetivo: dado una tasa de error tipo I , hay que usar un procedimiento para selecionar el conjunto de genes significantes de tal manera que error tipo I sea  Si el interes es controlar el FWER, entonces el p-value ajustado para la hipotesis H j is: p j * = inf {  ’: H j is rechazado al FWER  } Hipotesis H j es rechazada al FWER  si p j *  

10 Notacion Para hipotesis H j, j = 1, …, m prueba estadistica calculada: t j p-value no ajustado: p j Ordernamiento de t j (absoluto) observado: {r j } tal que |t r 1 |  |t r 2 |  …  |t r m | Ordernamiento de p j observado: {r j } tal que |p r 1 |  |p r 2 |  …  |p r m |

11 Control de la FWER Metodo de Sidak Rechazar H j con p j  1-(1-  ) 1/m, p-value ajustado p j * = 1-(1-p j ) m Metodo de Bonferroni Rechazar H j con p j  /m, p-value ajustado p j * = min (mp j, 1) Los metodos de Sidak y Bonferroni son faciles de calcular pero son muy conservativos. Para un numero grande de genes, el nivel individual de significancia por gene se vuelve bien pequeno bien rapidamente. Estos metodos son llamados de un solo paso porque usan el mismo ajuste de multiplicidad para todas las hipotesis y no usan el ordenamiento de los p-values observados Los p-values ajustados pueden ser obtenidos usando la funcion mt.rawp2adjp de la libreria multtest

12 Control de la FWER Metodo step-down de Holm (1979) Los p-values ajustados se definen como p* r1 =mp r1 p* ri =max(p* ri-1,(m-i+1)pr i ) para 2  i  m con p* ri  1 tomados como 1. O sea si j*=min{j:P rj  >  /(m-j+1)}, rechazar H j para j=1,…j*-1. Los p-values ajustados estan dados por: p r j * = max k = 1…i {min ((m-k+1)p r k, 1)} Los p-values ajustados pueden ser obtenidos usando la funcion mt.rawp2adjp de la libreria multtest

13 Control de la FWER Por ejemplo supongamos que los 4 primeros p-values ordenados son: 0.0006950, 0.0008659, 0.00149758, 0.00417016 de un total de m=1000 p- values. Entonces los p-values ajustados correspondientes seran: p* 1 =.6950, p* 2 =max(.6950, 999*0.0008659)=0.86510, p* 3 =max(.86510,998*0.00194758)=1.494585. Luego se toma p* 3 =1 pues el p-value no debe exceder a 1. P* 4 =max(1,997*0.004170)==1 y todos los que siguen tambien se toman como 1.

14 Control de la FWER Metodo step up de Hochberg (1988) Los p-values ajustados se definen como p* rm =p rm p* ri =min(p* ri+1,(m-i+1)p ri ) para i=m-1,…1 con p-values  1 tomados como 1. O sea si j*=min{j:P rj  /(m-j+1)}, rechazar H j para j=1,…j*. Los p-values ajustados vienen dados por p r j * = min k = j..m {min ((m-k+1)p r k, 1) } Los p-values ajustados pueden ser obtenidos usando la funcion mt.rawp2adjp de la libreria multtest

15 Control de la FWER P-value ordenado SidakBonferroni Holm  Hochberg  p r1 1-(1-  ) 1/m  /m p r2 1-(1-  ) 1/m  /m  /(m-1) ….……………. p rj 1-(1-  ) 1/m  /m  /(m-j+1) ….………..…… p rm 1-(1-  ) 1/m  /m 

16 Ejemplo; Aplicacion a datos de Golub 38 pacientes de Leucemia, 27 con leucemia aguda linfoblastica (ALL) y 11 con leucemia aguda meloide (AML). Inicialmente hay 6817 genes que se reducen a 3051 aplicando cierto criterios de exclusion para valores de expresion

17 Ejemplo (cont) 1045 genes expresados con p-values

18 Ejemplo (cont) En ambos casos se detectan 98 genes expresados

19 Top 10 genes twilcoxsidakbonferroniHolmHochberg 829896829 3782124378 21248292124 8082670808 248929392489 394 26707662670 10098081009 199518341995 9372600937

20 Pruebas basadas en permutaciones Estimar la distribucion conjunta de los estadisticos de prueba T 1,..,T m, donde m es el número de genes, mediante permutaciones de las columnas de la matriz de expression genética Labels originales de los grupos estadistico t Labels permutados de los grupos 1.45 1.34 1.89 -1.33 -0.78 2.17

21 Algoritmo de Permutación para p- values-no ajustados Calcular las pruebas tj para cada una de las hipótesis Hj. Hacer B permutaciones (B aprox 1000) de las columnas de la matriz de expresión genética Para la b-ésima permutacion b=1,2,….B A) Calcular las pruebas estadisticas t 1,b,……t m,b B) El p-value estimado por permutacion para la prueba de hipotesis Hj está dado por la proporcion de |t ib | ‘s que son mayores que |t i | donde I (.) representa la función indicadora La libreria multtest tiene la funcion mt.sample.teststat que calcula el test estadistico por permutaciones pero lo hace vector por vector

22 Westfall & Young (1993) p- values ajustados En cada paso hace pequenos ajustes Toma en cuenta la distribucion conjunta (dependencia ) de la pruebas estadisticas Menos conservativo que los anteriores metodos de ajustar p-values. Puede ser estimado por remuestreo resampling pero tarda bastante (especialmente la version minP)

23 Metodo maxT de Westfall & Young Ordenar los t-values observados Para la b-esima (b=1…B) permutacion de las columnas del conjunto de datos calcular a) los t j,b values para cada H j (j=1,…m) b) Calcular los maximos consecutivos de la prueba estadistica. u m,b =|t rm,b | u j,b =max(u j+1,b,|t rj,b |) para j=m-1,……1 Calcular los pvalues ajustados usando

24 Metodo maxT de Westfall & Young Tambien es conocido como step-down maxT Un formula equivalentemente para los p-values ajustados es p r j * = max k = 1…j { p(max l  { r k… r m} |T l | ≥ |t r k |  H 0 C )} T l es el valor de la prueba estadistica correspondiente a la l-esima hipotesis.

25 Ejemplo del metodo MaxT gene|t| 10.1t r5 40.2t r4 52.8t r3 23.4t r2 37.1t r1 Gene|t b ||u b |I(u b >|t|) 11.3 1 40.81.31 53.0 1 22.13.00 31.83.00  P*=  /B 935.935 876.876 138.138 145.145 48.048 Genes ordenados segun su valor t Valores u para la b-esima permutacion P-Values ajustados

26 Explicacion de los valores de la tabla De la primera tabla hay que guardar los indices de los ordenamientos. Esto es, r 1 =3, r 2 =2,r 3 =5,r 4 =4,r 5 =1. Luego, de la segunda tabla u r5 =|t r5, b| =1.3, Tambien U r4 =max(u r5,|t r4 |)=max(1.3,0.8)=1.3 U r3 =max(u r4,|t r3 |)=max(1.3,3.0)=3.0 U r2 =max(u r3,|tr2|)=max(3.0,2.1)=3.0 y U r1 =max(u r2,|tr1|)=max(3.0,1.8)=2.0

27 Metodo minP de Westfall & Young Tambien llamado step-down minP En este caso los p-values ajustados vienen dados por: p* r1 = p(min l  { r 1… r m} P l  p r 1  H 0 C )} p r j * = max(p rj-1,p(min l  { r j… r m} P l  p r j  H 0 C )} para j=2,…m P l es la variable aleatoria para el p-value de la l- esima hipotesis. Por ejemplo, P l ~U(0,1)

28 maxT vs. minP Los p-values ajustados por maxT y minP son los mismos cuando las pruebas estadisticas son identicamente distribuidas (id) maxT es mas rapida computacionalmente que minP maxT es mas poderosa en los casos en que el numero de genes m es grande y numero de arreglos n es pequeno.

29 Top 10 genes maxTminP 2124 829 896 766 2600 2939 1995 2386 717 2489

30 El criterio de controlar el FWER es demasiado conservativo esto significa que muchos genes que son diferencialmente expresados podrian no ser detectados. Para el ejemplo de Golub solo se detectan 98 genes como expresados con los ajusten de bonferoni y Holm. El maxT de Westfall detecta 91, pero el minP no detecta ninguno. El criterio FDR trata de resolver este problema

31 Control de la FDR Benjamini & Hochberg (1995): step-up Asumiendo que el FDR es de nivel , se rechaza H j para j=1,…j*, donde j*=max{j: p rj

32 Top 10 genes BHBY 829 378 2124 808 2489 394 2670 1009 1995 937 Bejamini y Hocberg detecta 681 genes expresados y Bejnjamini y Yuketilei detecta 269. Tambien se puede ajustar los p- values obtenidos con maxT o minP res1=mt.maxT(golub,golub.cl) Rawp=res1$rawp[order(res1$index)] # Permutation adjusted p-values for simple multiple testing procedures Res2=mt.rawp2adjp(rawp,”BH”)

33 Otras propuestas para multiple testing ‘Significance Analysis of Microarrays (SAM)’ (2 versions) –Tusher et al. (2001) –Efron et al. (2001), basada en empirical Bayes SAM tambien estima el ‘FDR’, pero este es definido como E(V|H 0 C )/R y no como E(V/R) La libreria siggenes de Bioconductor encuentra los genes diferenciados por SAM

34 Analisis de significancia de Microarrays (SAM) Tusher et al. (2001) Does not assume normal distribution – Instead, p-values computed via values computed via permutation Test statistic: similar to t Test statistic: but with modified with modified variance estimate – Improved for small experiments.

35 Analisis de significancia de Microarrays (SAM) First, compute test statistic per gene for the observed data. donde s* p es un estimado de la desviacion estandar combinada y s o es el coeficiente de variacion minimo de la expresion de todos los genes. Compute average of test statistics distribution, over all statistics distribution, over all permutations – This gives estimate of distribution, if treatment has no effect. Un gen es considerado significamente diferenciado si su distancia con respecto a la media de su distribucin excede un threshold  SAM es facil de implementar El estimado usa todos los genes, luego si uno deellos se afecta por el tratamiento tambie se afectara el estimado. No es confiable en experimentos pequenos.

36 Ejemplo de SAM Row d.value stdev rawp q.value R.fold 1 829 8.165222 0.2958251 0 0 7.2771792 2 2124 7.964784 0.1778697 0 0 3.3953035 3 2600 6.102371 0.1911219 0 0 2.6686992 4 2664 5.975750 0.3918749 0 0 4.7229540 5 766 5.970848 0.1731333 0 0 2.4972299 6 2489 -5.726212 0.2154975 0 0 0.3449532 7 717 -5.704438 0.2068956 0 0 0.3450401 8 1995 -5.696514 0.1933259 0 0 0.3735687 9 2939 -5.576921 0.1650727 0 0 0.4132925 10 2663 5.547021 0.4178283 0 0 5.4551523 11 378 5.408458 0.3024200 0 0 4.4230278 12 1778 5.336856 0.2215924 0 0 2.7817426 13 1911 5.170084 0.1897508 0 0 2.3574207 14 1413 5.168875 0.2809704 0 0 3.2926668 15 808 5.139462 0.1819399 0 0 2.4278709

37 SAM usando empirical Bayes (EBAM) Aqui se estima la probabilidad posterior p 1 (Z) =1-p o f o (Z)/f(Z) de que un gen con score Z sea expresado. La razon f o (Z)/f(Z) es estimada usando una regresion logistica basada en las densidades relativas de los scores Z i

38 Ejemplo de SAM (empirical Bayes) # El valor optimo del factor a0 is determinado, donde # los posibles valores de a0 son 0 y los quantiles 0, 0.05 y 0.1 # de la desviacion estandard de los genes. Hacer rand=123 find.out=find.a0(golub,golub.cl,alpha=c(0,0.05,0.1),rand=123) # Una vez que se establece el valor optimo de a0,se efectua #un analisis por Empirical Bayes. >ebam.out=ebam(find.out,gene.names=golub.gnames[,3]) >cat("\n el numero de genes diferenciados es",length(ebam.out$row.sig.genes),"\n") el numero de genes diferenciados es 714 > cat("\n los top 10 genes son\n") los top 10 genes son > ebam.out$row.sig.genes[1:10] [1] 2489 394 1995 2939 717 1042 2702 523 1811 849

39 References Alizadeh et al. (2000) Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature 403: 503-511 Benjamini and Hochberg (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. JRSSB 57: 289-200 Benjamini and Yuketieli (2001) The control of false discovery rate in multiple hypothesis testing under dependency. Annals of Statistics Efron et al. (2000) Microarrays and their use in a comparative experiment. Tech report, Stats, Stanford Golub et al. (1999) Molecular classification of cancer. Science 286: 531-537

40 References Hochberg (1988) A sharper Bonferroni procedure for multiple tests of significance. Biometrika 75: 800-802 Holm (1979) A simple sequentially rejective multiple testing procedure. Scand. J Statistics 6: 65-70 Tusher et al. (2001) Significance analysis of microarrays applied to transcriptional responses to ionizing radiation. PNAS 98: 5116 -5121 Westfall and Young (1993) Resampling-based multiple testing: Examples and methods for p-value adjustment. New York: Wiley Yuketieli and Benjamini (1999) Resampling based false discovery rate controlling multiple test procedures for correlated test statistics. J Stat Plan Inf 82: 171-196