Statistical learning

P(A P(A) P(B) P(A

quindi ∪ ≤

B) +

P(A P(A) P(B)

m E , ..., E

Per eventi m

1 m m

[ X

≤

E )

P P(E

i i

i=1 i=1

abbiamo l’uguaglianza se gli eventi sono disgiunti. 90

sim <- function(n,p){

X <- cbind(1,matrix(runif(n*(p-1)),ncol=(p-1)))

XX <- as.data.frame(X)

beta <- c(1,0,0)

y <- X %*% beta + rnorm(n)

anova(lm(y~1),lm(y~.,XX))$"Pr(>F)"[2]

}

res <- replicate(1000,sim(n=10,p=3))

hist(res,freq=F)

abline(h=1,col=2) Histogram of res

1.2

1.0

0.8

Density 0.6

0.4

0.2

0.0 0.0 0.2 0.4 0.6 0.8 1.0

res

plot(ecdf(res))

curve(punif,0,1,add=T,col=2) 91

ecdf(res)

1.0

0.8

0.6

Fn(x) 0.4

0.2

0.0 0.0 0.2 0.4 0.6 0.8 1.0

x

Bonferroni α

α̃ =

Rifiutare tutte le ipotesi nulle con p-value inferiore a m

≤

F W ER α

Garantisce senza nessuna assunzione.

Dimostrazione: [

α α α

X

≤ ≤ ≤ · ≤

≤ p m α

p P

P i i 0

m m m

i∈T i∈T

Tre disuguaglianze:

1) Usa la disuguaglianza di Boole. ≤ ≤ ∈ T

t) t i

2) Usa la proprietà dei p-values relativi a ipotesi nulle vere: per .

P(p

i

≤

m m

3) Usa 0

P-values aggiustati: α.

Rifiutare tutte le ipotesi nulle con p-value aggiustato inferiore a

Il p-value aggiustato per la procedura Bonferroni è: ·

p̃ min(m p ,

= 1)

i i

92 α

α̃

Bonferroni è sempre valido e rifiuta tutte le ipotesi nulle con un p-value inferiore ad = .

m

p̃ corrisponde al p-value aggiustato e dalla sua forma risulta essere sempre superiore al p-value non aggiustato,

i ·

m p

poichè prende il minimo tra e 1.

i

Prostate cancer study dataset

n

• Studio effettuato su = 102 uomini, 50 controls contro 52 cases.

m

• = 6033 geni.

H i

• : attività del gene è lo stesso in casi e controlli.

i

• Dopo il codice per importare, possiamo osservare un boxplot dei livelli di espressione dell’i-esimo gene.

prostmat <- Learning/Dataset/prostmat.csv")

read.csv("F:/Statistical

dim(prostmat)

## [1] 101 102

X <- t(prostmat)

Holm = Bonferroni sequenziale

α

α̃

Iniziare con = m

Ripetere i passi 1) e 2) ≤ α̃

1) Rifiutare tutte le ipotesi con p-value

α

α̃

2) Aggiornare = m−r

r

dove è il numero di ipotesi nulle rifiutate finora.

HOLM rifiuta sempre di più che bonferroni, quindi è uniformemente più potente di Bonferroni.

93

###Holm (procedura step-down) 94

1) Ordinare i p-values dal più piccolo al più grande

≤ ≤ ≤

p p ... p

(1) (2) (m)

e le ipotesi corrispondenti ≤ ≤ ≤

H H ... H

(1) (2) (m)

{1,

j, ..., m}

2) Determinare il più piccolo indice tra 2, tale che

α

p >

(j) −

m j +1

H i < j

3) Rifiutare tutte le ipotesi tali che

(i)

###Hochberg (procedura step-up) 95

1) Ordinare i p-values dal più piccolo al più grande

≤ ≤ ≤

p p ... p

(1) (2) (m)

{1,

j, ..., m}

2) Determinare il più grande indice tra 2, tale che

α

≤

p (j) −

m j +1

≤

H i j

3) Rifiutare tutte le ipotesi tali che

(i)

La procedura di Hochberg rifiuta come o più di Holm, ma bisogna assumere che i p-values

{p ∈ T }

, i soddisfano la diseguaglianza di Simes (dipendenza positiva).

i vado verso destra: se sto sotto, continuo. Nel momento in cui vado sopra, mi fermo e rifiuto

Step-down

tutto quello che ho prima.

Nella invece vado da destra a sinistra e se sto sopra, allora continuo. Quando rifiuto, andrò sotto,

step-up

allora mi fermo e rifiuto tutto quello che ho a sinistra. Quindi nella secondo rifiuterò di più.

96

Relazioni logiche

ANOVA ad una via: 2 2 2

∼ ∼ ∼

Y N , σ Y N , σ Y N , σ

(µ ), (µ ), (µ )

1 1 2 2 3 3

Confronti a coppie: H µ µ H µ µ H µ µ

: = : = : =

12 1 2 13 1 3 23 2 3

Relazioni logiche tra ipotesi nulle:

H H H

Ad esempio, se è falsa, allora e non possono essere contemporaneamente vere.

12 13 23

Shaffer = Holm con relazioni logiche

α

α̃

Iniziare con = m

Ripetere i passi 1) e 2) ≤ α̃

1) Rifiutare tutte le ipotesi con p-value

αs

α̃

2) Aggiornare =

s

dove è il numero massimo di ipotesi nulle che possono essere contemporaneamente vere, assumendo che le

ipotesi rifiutate finora siano ipotesi false.

Shaffer è uniformemente più potente di Holm.

Confronti a coppie:

H µ µ p µ µ p µ µ p

: = = 0.01H : = = 0.54H : = = 0.04

12 1 2 12 13 1 3 13 23 2 3 23

α̃

Valore critico

m α

1) = 3, = 5% 0.05 ≈

α̃ H

= 0.0167 e rifiuto 12

3 →

H H H s

2) Se è falsa, allora al massimo una tra e può essere vera = 1

12 13 23

0.05 H

α̃ e rifiuto

= 23

1 →

H H H s

3) Se e sono false, allora può essere vera = 1

12 23 13 0.05

α̃ = ma non rifiuto alcuna ipotesi

1

Ipotesi nulle dei confronti a coppie

H µ µ H µ µ H µ µ

: = : = : =

12 1 2 13 1 3 23 2 3

Ipotesi nulla globale ∩ ∩

H H H H µ µ µ

= : = =

123 12 13 23 1 2 3

Relazioni logiche

H H , H H

Se è falsa, al massimo una tra e può essere vera.

123 12 13 23

Vediamo graficamente ANOVA ad una via con 3 gruppi

97

Vado a fare test d’ipotesi sulla congiunta, ad esempio il test se rifiuto, allora vado ad analizzare le singole

F,

ipotesi e almeno una sarà vera.

Osserviamo un altro dataset:

prostz <- Learning/Dataset/prostz.txt", quote="\"", comment.char="")

read.table("F:/Statistical

pvals <- = F)*2

pnorm(abs(prostz$V1),lower.tail

# Histogram col="lightgray", main="", xlab="p-value", freq=F)

hist(pvals,20,

abline(h=1,lty=2) 98

1.5

1.0

Density 0.5

0.0 0.0 0.2 0.4 0.6 0.8 1.0

p−value

#grafico che ordina i p-values

Simulazione (FWER control)

m = length(pvals)

alpha = 0.05

<= alpha)

sum(pvals

## [1] 478

# Bonferroni <= alpha)

sum(p.adjust(pvals,"bonferroni")

## [1] 3

# Holm and Hochberg <= alpha)

sum(p.adjust(pvals,"holm")

## [1] 3 <= alpha)

sum(p.adjust(pvals,"hochberg")

## [1] 3 xlim=c(1,20), ylim = xlab="i", ylab="p-value")

plot(1:m, sort(pvals), c(0,.0001),

alpha/(m- 1:m +1))

lines(1:m, 99

8e−05

p−value 4e−05

0e+00 5 10 15 20

i

FDR (procedure per il controllo)

Benjamini & Hochberg

1) Ordinare i p-values dal più piccolo al più grande

≤ ≤ ≤

p p ... p

(1) (2) (m)

{1,

j, ..., m}

2) Determinare il più grande indice tra 2, che soddisfa

jα

≤

p (j) m

Diversamente da bonferroni ho un valore critico differente.

≤ α̃ p

3) Rifiutare tutte le ipotesi con p-value = (j)

• La procedura assume che i p-values sono indipendenti o più generalmente con dipendenza positiva

BH

- positive regression dipendence on the subset of p-values of true null hypotheses (PRDS).

m

≤

F DR α

• Se è vera l’assunzione PRDS, garantisce

BH 0

m

• è una procedura

BH step-up. 100

101

#wilcoxon.test

#pvals2 <- apply(X,2, function(x) wilcox.test(x=x[1:50],y=x[-c(1:50)])$p.value)

#hist(pvals2,20)

#sum(pvals2 <= alpha)

#sum(p.adjust(pvals2,"BH") <= alpha)

In questo caso rifiutiamo di più, quindi controllo la percentuale di falsi positivi in valore atteso. Mi aspetterò

·

circa 3 falsi positivi, dato dal seguente calcolo: 61 0.05.

Benjamini & Yekutieli

1) E’ una variante di valida per qualsiasi forma di dipendenza tra i p-values.

BH iα , i ..., m

2) Rende più stringenti i valori critici = 1, moltiplicandoli per

m 1

1 m→∞

=

m 1

P log(m)

i=1 i

102

#wilcoxon.test

#pvals2 <- apply(X,2, function(x) wilcox.test(x=x[1:50],y=x[-c(1:50)])$p.value)

#hist(pvals2,20)

#sum(pvals2 <= alpha)

#sum(p.adjust(pvals2,"BY") <= alpha)

Qui non funziona molto bene, poichè il metodo è troppo stringente.

Simulazione (FDR control)

<= alpha)

sum(p.adjust(pvals,"BH")

## [1] 21 xlim=c(1,50), ylim = xlab="i", ylab="p-value")

plot(1:m, sort(pvals), c(0,.002),

alpha/(m- 1:m +1))

lines(1:m, alpha*(1:m)/m, col=2)

lines(1:m,

0.0020

p−value 0.0010

0.0000 0 10 20 30 40 50

i

#Benjamini & Yekutieli <= alpha)

sum(p.adjust(pvals,"BY")

## [1] 2 103

Lezione 10

FDP (stima e confidenza)

Stima FDP di Storey

Impostazione rifiutata:

R {H ≤

p t} R

Supponiamo di rifiutare = : con = #R.

i i 104

Intuizione: m t

≈

Dall’uniformità dei p-values sotto FDP nulla 0

R

m

Stima di 0 ≥ −

> λ}) m λ),

Idea: (1 così

E(#{p i 0 > λ}

#{p + 1

i

m̂ =

0 − λ

1

Stima di FDP (q-value) m̂ t t > λ}

#{p + 1

i

0

ˆ

F DP (t) = = − ≤

R λ t}

1 #{p

i

Figure 2:

Storey e Benjamini & Hochberg

Relazioni vicine:

Una via alternativa di costruzione di impostazione rifiutata:

BH

m̂

1) Stimare = 1 invece di una stima di Storey.

0 ˆ ≤

t F DP α.

2) Prendere il valore più grande come (t)

Un look alternativo di Storey:

Metodo Storey = controllo adattato FDR.

Un look alternativo di BH:

Stima conservativa di FDP.

Storey e dipendenza

Metodo delle stime dei momenti:

⇒

Solo dipendenzain media Non affetta da una struttura correlata.

Però:

La variabilità della stima può essere grande se i p-values sono correlati.

105

Standard errors non disponibili:

Disponibili solo per p-values indipendenti.

Simulazione (Prostate data)

#dataset già importato in precedenza

# Sorted p-values plot type="l", xlab=&