---
title: "Comparaciones múltiples"
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
format: 
    revealjs:  
      scrollable: true  
      echo: true  
      embed-resources: true
    pdf:  
      toc: true  
      number-sections: false
      colorlinks: true
---

##  Comparaciones múltiples

- Una formulación ($i = 1, \ldots, N$): 
    - $H_i$: El gen i **no tiene** una **expresión diferencial** entre 
      las condiciones consideradas.
    - $K_i$: El gen i **tiene** una **expresión diferencial** entre 
      las condiciones consideradas.
 
- Quizás es mejor:
    -  $H_i$: La expresión del  gen i no tiene asociación con la condición.
    -  $K_i$: La expresión del  gen i tiene asociación con la condición.


## Benjamini and Hochberg (1995)

--------------------------------------------------------------
Hipótesis nula        No rechazadas      Rechazadas     Total 
-------------------  ---------------    ------------   -------------
    Verdadera           U                 V              $N_0$    
    
    Falsa               T                 S              $N-N_0=N_1$  

    Total              N - R              R                  N         
-------------------------------------------------------------------
Table: 

- ¿Qué conocemos? Solamente la variable R y el número de hipótesis N.


# Tasas de error

## FWER: Family wise error rate
- Se define como: 
 $$FWER = P(V > 0) = P(V \geq 1).$$
- Es la tasa de error de uso clásico en Estadística. 
- Correcta para un número pequeño de hipótesis. 
- Es un criterio muy exigente si tenemos un número de hipótesis muy
     grande.  
- Notemos que nos fijamos en cometer al menos un error cuando
     con frecuencia tendremos decenas de miles de contrastes.   
- Como (mala) contrapartida muchos genes
     que tienen una expresión diferencial realmente no serán detectados.  

## FDR: false discovery rate 

- O tasa de falsamente rechazados
- Definimos: 
$Q = \frac{V}{R}$ si $R >0$ y $Q = 0$ en otro caso.
- Proporción de tests erróneamente rechazados.
- $$FDR = E (Q) = E \bigg ( \frac{V}{R} \bigg | R >0 \bigg ) P(R > 0).$$

- Una modificación importante: **pFDR** (Positive false discovery rate)
$$pFDR = E \bigg ( \frac{V}{R} \bigg | R >0 \bigg ).$$


## Relación entre las tasas de error tipo I

- Se verifican las siguientes desigualdades:
$$
FDR \leq FWER. 
$$

## p-valor ajustado

- Para cada contraste $H_i$ tendremos un p-valor $p_i$.
-  Podemos contrastar el contraste $H_i$ con un nivel de significación 
$\alpha_i$ rechazando la hipótesis nula $H_i$ si $p_i \leq \alpha_i$ y 
no rechazando en otro caso.
- Cuando consideramos simultáneamente todos los tests los valores
$\alpha_i$ serán distintos y por ello tendríamos que ir comprobando 
si la desigualdad $p_i \leq \alpha_i$ se verifica o no.
- La idea del p-valor ajustado es transformar el p-valor $p_i$ en otro 
valor $\tilde{p}_i$, el i-ésimo p-valor ajustado, de modo que sean
equivalentes: 
$$p_i \leq \alpha_i$$ y $$\tilde{p}_i \leq \alpha$$ siendo
$\alpha$ el valor que especificamos para controlar alguna de las
tasas de error tipo I.


## Método de Bonferroni

- Rechazamos $H_i$ si 
$$p_i \leq \frac{\alpha}{N}.$$
- El p-valor ajustado sería:
$$
\tilde{p}_i = \min\{N p_i, 1 \}
$$
-  Rechazamos $H_i$ si 
$$ \tilde{p}_i \leq \alpha.$$

## Método de Benjamini-Hochberg

1. Fijamos la tasa de error $\alpha$.
2. Para cada i (gen) aplicamos un test y obtenemos un p-valor $p_i$.
3. Ordenamos los p-valores
$$p_{r_1} \leq \ldots \leq  p_{r_N}.$$
4. Sea 
$$i^* = \max\{i: p_{r_i} \leq \frac{i}{N} \alpha\}$$
5. Rechazamos $H_{r_i}$ para $i=1, \ldots,i^*$. 
6. Si no existe $i^*$ entonces no rechazamos ninguna hipótesis. 

## Método de Benjamini y Yekutieli

-  Como en el anterior, $p_{r_1} \leq \ldots \leq  p_{r_N}$ 
son los p-valores originales ordenados.
- Los p-valores ajustados se definen como
  $$
  \tilde{p}_{r_i} = \min_{k=i,\ldots,N} \bigg \{ \min \big \{ 
  \frac{N \sum_{j=1}^N 1/j}{k} p_{r_k},1 \big \} \bigg \}.
  $$

## Lo aplicamos con dos grupos

- Leemos datos.

```{r}
#| label: cm1
pacman::p_load(Biobase)
data(gse21942,package="tamidata")
y = pData(gse21942)[,"FactorValue..DISEASE.STATE."]
head(exprs(gse21942))
dim(gse21942)
```

- Aplicamos los t-tests para cada gen.  
```{r}
#| label: cm2
tt = genefilter::rowttests(gse21942,y)
head(tt)
```
- Los p-valores originales los guardamos en p0.
```{r}
#| label: cm3
p0 = tt$p.value
head(p0)
```
o bien

```{r}
#| label: cm4
p0 = tt[,"p.value"]
```
Con un nivel de significación de $\alpha=0.01$ tendríamos
el siguiente número de características significativas.
```{r}
#| label: cm5
table(p0<=.01)
```

## 

-  Utilizamos el método de Benjamini-Hochberg. 

```{r}
#| label: cm6
p.BH = p.adjust(p0,method = "BH")
```

- Los p-valores ajustados son los siguientes

```{r}
#| label: cm7
View(cbind(p0,p.BH))
head(p.BH)
table(p.BH<.01)
```

# q-valor

## Definición 

- Propuesto por Storey (2002)
- Si consideramos un test determinado nos da la proporción esperada de 
falsos positivos en la que incurrimos cuando declaramos significativo ese test. 

## Estimación del q-valor 

- Fijamos un valor $t$ y consideremos que rechazamos $H_i$ cuando $P_i \leq t$.
- Definimos:
$$
V(t) = | \{P_i : P_i \leq t;  H_i \text{ es cierta}; i =1,\ldots,N \}|
$$
y 
$$
R(t) = | \{P_i : P_i \leq t;  i =1,\ldots,N\}|
$$
- pFDR se  puede aproximar con
$$
E \bigg [ \frac{V(t)}{R(t)} \bigg ] \approx \frac{E V(t)}{E R(t)}.
$$

##
- Estimamos 
$$
\hat{R}(t) = | \{p_i : p_i \leq t \}| 
$$
- Además 
$$
E V(t) = N_0 t. 
$$
- Estimamos $\pi_0 = N_0 /N$ con 
$$
\hat{\pi}_0 = \frac{|\{p_i: p_i > \lambda; i = 1,\ldots,N\}|}{N (1- \lambda)}.
$$

## 

- Podemos pues estimar pFDR con
$$
\widehat{pFDR} =  \frac{\hat{\pi}_0 N t}{|\{p_i: p_i \leq t\}|}.
$$
- El q-valor asociado a un contraste sería el mínimo valor de pFDR 
que se alcanza cuando el contraste es rechazado.
- El q-valor asociado al test $i$-ésimo sería
$$
q(p_i) = \min_{t \geq p_i} pFDR(t)
$$
y su estimador sería 
$$
\hat{q}(p_i) = \min_{t \geq p_i} \widehat{pFDR}(t).
$$

## q-valor y datos tamidata::gse21942

- Calculamos p-valores originales.

```{r}
#| label: cm8
tt = genefilter::rowttests(gse21942,
                           pData(gse21942)[,"FactorValue..DISEASE.STATE."])
pvalue = tt$p.value
```

- Calculamos los q-valores

```{r}
#| label: cm9
library(qvalue)
aa = qvalue(pvalue)
```

- Si simplemente queremos los q-valores los obtenemos con

```{r}
#| label: cm10
q.value = qvalue(pvalue)$qvalues
```


## Dibujos asociados al q-valor


```{r}
#| label: cm11
plot(aa)
```


