Multiple testing
Problema:
Ipotiziamo una collezione di ipotesi nulle \(\mathcal{H}=\{H_1, \ldots, H_m\}\) che vogliamo rigettare. Un numero sconosciuto \(m_0\) di queste risultano essere vere mentre altre \(m_1= m-m_0\) sono false. Chiamiamo l’insieme di ipotesi vere \(\mathcal{T} \subseteq \mathcal{H}\) e la rimenente collezione di ipotesi false \(\mathcal{F} = \mathcal{H} \setminus \mathcal{T}\), la proporzione delle ipotesi vere sarà quindi \(\pi_0 = m_0/m\).Obbiettivo:
L’obbiettivo della multiple testing procedure è quello di scegliere \(\mathcal{R} \subseteq \mathcal{H}\) di ipotesi da rigettareRisoluzione teorica:
Se abbiamo dei \(p\)-values \(p_1,\ldots,p_m\) per ogni ipotesi \(H_1,\ldots,H_m\) ci verrebbe da selezionare in modo intuitivo quei test che: \[\mathcal{R} = \{H_i: p_{i}\leq T\}\] rigettano tutte le ipotesi con un \(p\)-value minore di una soglia pari a \(T\). A livello teorico il set delle ipotesi rifiutate \(\mathcal{R}\) dovrebbe essere coerente \(\mathcal{F}\) il piu possibileDue tipi di errore:
I falsi positivi, o gli errori di tipo I, sono le ipotesi rifiutate che non sono false, i.e. \(\mathcal{R} \cap \mathcal{T}\); falsi negativi o errori di tipo II sono le false ipotesi che non siamo riusciti a rifiutare, i.e. \(\mathcal{F} \setminus \mathcal{R}\). Le ipotesi rifiutate sono talvolta denominate scoperte, quindi i termini true discovery e false discovery vengono talvolta utilizzati per respingimenti corretti e scorretti.
Possiamo riassumere il numero di errori che si verificano in una procedura di verifica delle ipotesi in una tabella di contingenza. Possiamo osservare \(m\) e \(R\), ma tutte le quantità nelle prime due colonne della tabella non sono osservabili.
| true (Ipotesi nulla \(H_0\) è vera) | false (Ipotesi alternativa \(H_A\) è vera) | Totale | |
|---|---|---|---|
| rejected (Test viene dichiarato significativo) | \(V\) | \(U\) | \(R\) |
| not rejected (Test viene dichiarato non significativo) | \(m_0-V\) | \(m_1-U\) | \(m-R\) |
| Totale | \(m_0\) | \(m_1\) | \(m\) |
dove:
- \(m\) è il numero totale delle ipotesi testate
- \(m_0\) è il numero delle ipotesi nulle vere
- \(m_1=m-m_0\) è il numero delle ipotesi alternative vere
- \(V\) è il numero di falsi positivi (errore di tipo I) (chiamato anche “false scoperte”)
- \(U\) è il numero di veri positivi (chiamato anche “vere scoperte”)
- \(m_1-U\) è il numero di falsi negativi (errore di tipo II)
- \(m_0-V\) è il numero di veri negativi
- \(R = V + U\) è il numero di ipotesi di nulle rifiutate (chiamato anche “scoperte”, true o false)
Gli errori di tipo I e di tipo II sono in concorrenza diretta e occorre un compromesso tra i due. Se rifiutiamo più ipotesi, di solito abbiamo più errori di tipo I ma meno errori di tipo II e viceversa.
ESEMPIO:
\(Z_i \sim N(\mu_i,1)\) è la statistica test
\(H_i: \mu_i = 0\) vs \(\bar{H}_i: \mu_i > 0\)
Rifiutiamo \(H_i\) if \(Z_i \geq z_{1-\alpha}\) o equivalentemente quando \(p_i = 1-\Phi(Z_i) \leq \alpha\)
dove \(z_{1-\alpha}\) il quantile di ordine \(1-\alpha\) di una \(N(0,1)\) e \(\Phi\) è la cdf (funzione di ripartizione) di una \(N(0,1)\)
set.seed(123)
alpha = 0.05
m = 100
m0 = 80
effect = 3
setT = sample(1:m, size=m0, replace=F)
stats <- rnorm(m)
stats[-setT] <- stats[-setT] + effect
pvals = pnorm(stats, lower.tail = FALSE)
setR = which(pvals <= alpha)
# setR = which(stats >= qnorm(1-alpha)) # equivalently
table( rejected= 1:m %in% setR,
true = 1:m %in% setT)## true
## rejected FALSE TRUE
## FALSE 3 75
## TRUE 17 517 + 5 = 22 sono l’ipotesi rifiutata mentre 3 + 75 sono l’ipotesi accettata. Ci sono 5 errori di tipo I (5 sbagliati respinti su 22) e 3 errori di tipo II
Familywise error rate
Numerosi metodi di prova cercano di respingere il maggior numero possibile di ipotesi mantenendo nel contempo una certa misura degli errori di tipo I. Il modo più classico di controllo per test multipli è il Familywise error rate (FWER). Il FWER è la probabilità che il set rifiutato contenga errori:
\[\mathrm{FWER} = \mathrm{P}(V > 0)\]
Il FWER è controllato al livello $ $ quando il set $ $ (cioè la soglia $ T $) viene scelto in modo tale che: \[ \mathrm {FWER} \leq \alpha \]
Bonferroni
Il metodo di Bonferroni controlla FWER al livello \(\alpha\) rifiutando le ipotesi solo se hanno un valore p inferiore a \(\alpha / m\):
\[\mathcal{R} = \left\{H_i: p_{i}\leq \frac{\alpha}{m}\right\}\]
Questo aggiustamento a steps singoli della soglia di significatività è il metodo di test multiplo più semplice (non richiede alcuna ipotesi).
Il p-value aggiustato secondo la procedura di Bonferroni corrisponde:
\[\tilde{p}_i = \min(mp_i, 1)\]
e il metodo Bonferroni rifiuta le ipotesi nulle che hanno un p-value aggiustato più piccolo di \(\alpha\):
\[\mathcal{R} = \{H_i: \tilde{p}_{i}\leq \alpha \}\]
Bonferroni proof
Le assunzioni sui \(p\)-values valgono spesso solo sulle ipotesi vere che chiamiamo $ q_1, , q_{m_0} $, se le loro ipotesi corrispondenti sono vere, questi \(p\)-values sono distribuiti uniformemente tra 0 e 1, oppure possono essere stocasticamente maggiori di uniformi: \[
\mathrm{Pr}(q_i \leq u ) \leq u
\]
Per motivare correttamente il metodo Bonferroni, dovremmo considerarlo come un corollario alla disuguaglianza di Boole, che dice che per ogni raccolta di eventi $ E_1, ldots, E_k $, abbiamo:
\[\mathrm{P}\big(\bigcup_{i=1}^k E_i\big) \leq \sum_{i=1}^k \mathrm{P}(E_i)\]
Dalla diseguaglianza di Boole deriva la probabilità che esiste qualche \(i\) per il quale \(q_i \leq \alpha/m\) is given by \[\Pr\big(\bigcup_{i=1}^{m_0} \{ q_i \leq \alpha/m \}\big) \leq \sum_{i=1}^{m_0} \mathrm{P}( q_i \leq \alpha/m ) \leq m_0\frac{\alpha}{m} \leq \alpha\]
Questo metodo commette errori solo del I tipo se \(q_i \leq \alpha/m\) per qualche \(i\), questo dimostra il controllo FWER al livello \(\alpha\) per il metodo Bonferroni.
set.seed(123)
alpha = 0.05
m = 100
m0 = 80
effect = 3
setT = sample(1:m, size=m0, replace=F)
stats <- rnorm(m)
stats[-setT] <- stats[-setT] + effect
pvals = pnorm(stats, lower.tail = FALSE)
setR = which(pvals <= alpha/m)
# setR = which(stats >= qnorm(1-alpha/m)) # equivalently
table( rejected= 1:m %in% setR,
true = 1:m %in% setT)## true
## rejected FALSE TRUE
## FALSE 16 80
## TRUE 4 0Ora abbiamo errore di tipo I 0 ma 16 errori di tipo II, questi ultimi sono aumentati.
Magnitude of Bonferroni’s threshold (Attrazione della soglia di Bonferroni)
\(Z_i \sim N(\mu_i,1)\) is the test statistic
\(H_i: \mu_i = 0\) vs \(\bar{H}_i: \mu_i > 0\)
Bonferroni rifiuta \(H_i\) se \(Z_i \geq z_{1-\alpha/m}\).
Ponendo \(\alpha\) fisso, possiamo mostrare che per un \(m > > 0\) \[z_{1-\alpha/m} \approx \sqrt{2\log m}\]
Non c'è asintoticamente alcuna dipendenza da \(\alpha\)
alpha = 0.05
m = 10^(2:12)
threshold = qnorm(1-alpha/m)
approx = sqrt(2*log(m))
plot(m,threshold,type="l",log="x")
lines(m,approx, col=2)