Laboratorium 11

Testowanie wielokrotne

Mamy \(m\) problemów testowania hipotezy \(H_0^i\), przeciwko hipotezie \(H_1^i\), \(i\in\{1, \dots, m\}\). Niech \(p_1, \dots, p_m\) oznaczają p-wartości dla \(m\) testów. Z definicji p-wartości zachodzi, że \[P_{H_0^i}(p_i\leq \alpha)\leq \alpha. \]

Przy testowaniu wielokrotnym pojawia się tak zwany problem testowania wielokrotnego. Otóż prawdopodobieństwo, że któryś z testów odrzuci prawdziwą hipotezę zerową wynosi \(1-(1-\alpha)^m\). Zauważmy, że \[1-(1-\alpha)^m \xrightarrow{m\rightarrow\infty}1.\]

Zadanie 1

Wygenerować \(m=500\) niezależnych prób prostych po \(20\) obserwacji z rozkładu \(N(d,1)\) gdzie dla ustalonej próby \(d\) jest stałe i \(d=1\) dla \(30\%\) prób oraz \(d=0\) dla \(70\%\) prób. Obliczyć wszystkie \(m\) p-wartości testu t-Studenta \(H_0: d=0\), przeciwko \(H_1: d\neq 0\) i sporządzić wykres punktowy \((i, p_{(i)})\), \(i=1, \dots, m\). Ile p-wartości przyjęło wartość poniżej \(\alpha=0.1\)?

n<-20
m<-500
a<-0.1
X<- matrix(0,nrow=n,ncol=m)
p<-numeric(m)
for (i in 1:m){
  if(i<=0.3*m){
    X[,i]<-rnorm(n,1,1)
  } else {
      X[,i]<-rnorm(n,0,1)
      }
  p[i]<-t.test(X[,i],mu=0,alternative="two.sided")$p.value
}

sp<-sort(p)
plot(sp)

A <- function(t){
  0*t+a
  }
plot(A,1,m,ylim=c(0,0.1))
points(sp,cex = 0.5)

r <- sum(p <= a)
r

## [1] 182

Niech \(V\) oznacza liczbę błędnie odrzuconych hipotez, a \(R\) liczbę wszystkich odrzuconych hipotez.

Rozważamy dwie miary dla procedur wielokrotnego testowania:

1. FWER (familywise error rate) \[FWER=P(V\geq 1).\]

2. FDR (false discovery rate) \[FDR=E\left(\frac{V}{\max\{R,1\}}\right).\]

Procedura Boniferroniego

Odrzucamy \(H_0^i\), gdy \(p_i\leq \frac{\alpha}{m}\). Dla procedury Boniferroniego zachodzi \(FWER\leq \alpha\).

Warto zauważyć, że moc testu maleje (najczęściej wykładniczo) wraz z poziomem istotności, czyli szansa na odrzucenie fałszywej hipotezy maleje. Dlatego też, procedury Boniferroniego nie używa się często w praktyce.

Zadanie 1 - cd.

B<-function(t){
  0*t+a/m
  }
plot(B,1,m,ylim=c(0,0.001),col="red")
points(sp,cex = 0.5)

rB <- sum(p <= a / m)

print(rB)

## [1] 77

Procedura Holma i Hochberga

Procedura Holma

Zachodzi dla niej \(FWER\leq \alpha\). Oznaczmy \(\left(p_{(1)}, \dots, p_{(m)}\right)\) - wektor statystyk pozycyjnych dla p-wartości oraz \(H_0^{(1)}, \dots, H_0^{(m)}\) odpowiadające im hipotezy zerowe.

W tym przypadku odrzucamy wszystkie hipotezy \(H_0^{(1)}, \dots, H_0^{(k)}\), gdy \(p_{(1)}\leq \frac{\alpha}{m}\), \(p_{(2)}\leq \frac{\alpha}{m-1}\), \(\dots\), \(p_{(k)}\leq \frac{\alpha}{m-k+1}\). Natomiast, nie odrzucamy, gdy \(p_{(k+1)}>\frac{\alpha}{m-k}\).

Procedura Holma jest lepsza od procedury Boniferroniego, ponieważ odrzuci ona wszystkie hipotezy, które odrzuciła procedura Boniferroniego, albo więcej.

Procedura Hochberga Dla procedury Hochberga zachodzi \(FWER\leq \alpha\) przy założeniu, że p-wartości \(p_1, \dots, p_m\) są niezależne. To znaczy, że testy wykonywane są dla różnych, niezależnych danych.

Dla tej procedury odrzucamy wszystkie hipotezy \(H_0^{(1)}, \dots, H_0^{(k)}\), gdy \(p_{(k)}\leq \frac{\alpha}{m-k+1}\). Natomiast, odrzucamy gdy \(p_{(k+1)}>\frac{\alpha}{m-k}\), \(p_{(k+2)}>\frac{\alpha}{m-k+1}\), \(\dots\), \(p_{(m)}>\alpha\).

Procedura Hochberga jest lepsza od procedury Holma, ponieważ nie zatrzymuje się na pierwszym odrzuconym wyniku, ale sprawdza dalej, co sprawia, że może odrzucić więcej fałszywych hipotez zerowych.

Zadanie 1 - cd.

Metoda Hochberga

H <- function(x){
  a/(m-x+1)
}
plot(B,1,m,ylim=c(0,0.001),col="red")
points(sp,cex = 0.5)
curve(H,1,m,add=TRUE,col="blue")

rHc<-0
for (i in 1:m){
  if(sp[i]<=a/(m-i+1)){
    rHc=rHc+1
    }
}
rHc

## [1] 79

Procedura Benjaminiego-Hochberga

Dla tej procedury zachodzi \(FDR\leq \alpha\). Dodatkowo zakładamy, że p-wartości odpowiadające prawdziwym hipotezom są niezależne od pozostałych. Dla p-wartości odpowiadającym fałszywym lub prawdziwym hipotezom możemy mieć dowolną zależność.

W tym przypadku, odrzucamy \(H_0^{(1)}, \dots, H_0^{(k)}\), gdy \(p_{(k)}\leq \frac{k\alpha}{m}\), natomiast nie odrzucamy, gdy \(p_{(k+1)}>\frac{(k+1)\alpha}{m}\), \(\dots\), \(p_{(m)}>\frac{m\alpha}{m}\).

Procedura Benjaminiego-Hochberga jest lepsza od Hochberga.

Zadanie 1 - cd

BH<-function(x){
  a*x/m
  }
plot(BH,1,m,ylim=c(0,0.1),col="green")
points(sp,cex = 0.5)
curve(B,1,m,add=TRUE,col="red")
curve(H,1,m,add=TRUE,col="blue")
curve(A,1,m,add=TRUE)

rBH<-0
for (i in 1:m){
  if(sp[i]<=a*i/m){
    rBH=rBH+1
    }
}
rBH

## [1] 162

Wszystkie powyższe metody możemy wywołać używając gotowej funkcji p.adjust.

sum(p.adjust(sp, method="bonferroni")<0.1)

## [1] 77

sum(p.adjust(sp, method="holm")<0.1)

## [1] 79

sum(p.adjust(sp, method="hochberg")<0.1)

## [1] 79

sum(p.adjust(sp, method="fdr")<0.1)

## [1] 162