Powtórka

W pierwszym zadaniu będziemy kolejny raz losować serię prób statystycznych o ustalonym rozmiarze.
Dotychczas robiliśmy to generując ciąg zmiennych losowych i ustawiając wylosowane wartości w macierz.
Wygodniejszym sposobem jest użycie specjalnych funkcji do powtarzania określonych wartości. W R są dwie takie funkcje: rep() oraz replicate(). Różnicę pomiędzy nimi najlepiej wyjaśni poniższy przykład: - Funkcja rep(rnorm(10), 20) najpierw wygeneruje wektor rnorm(10), a następnie skopiuje go 20 razy i sklei kopie w wektor długości 200. - Funkcja replicate(20, rnorm(10)) wywoła 20 razy funkcję rnorm(10) i ustawi wynik w macierz wymiaru 10 x 20, w której każda kolumna zawiera wynik wywołania funkcji rnorm(10).

Do otrzymania wektorów z powtórzonymi wartościami najlepiej używać funkcji rep. Do otrzymania prób statystycznych korzystamy natomiast z funkcji replicate. Zwróć uwagę, że w rep() najpierw podajemy co chcemy powtarzać, a następnie liczbę powtórzeń, podczas gdy w replicate() robimy odwrotnie.

Zadanie 1. Zadanie przykładowe. Wygeneruj 1000 prób po 10 obserwacji z rozkładu normalnego o średniej 0.5 i wariancji 1 i kolejne 1000 prób o średniej -0.5 i takiej samej wariancji. Porównaj moc testu t Studenta i testu rang Wilcoxona na poziomie istotności \(\alpha = 0.05\). Jakie są założenia tych testów? Czy w tym przypadku są spełnione? Który test działa w tej sytuacji lepiej i dlaczego? Porównaj swoje wyniki dla testu t Studenta z kalkulatorem mocy tego dostępnym tutaj.

Wskazówka 1: Funkcje do przeprowadzania testów w R zwracają listy zawierające wiele przydatnych informacji. Na przykład, mając wektory X oraz Y, p-wartość z testu t Studenta na tych dwóch wektorach możemy otrzymać pisząc t.test(X, Y)$p.value lub wynik_testu <- t.test(X, Y), a następnie wynik_testu$p.value.

Przypomnienie z wykładu: Poziom istotności \(\alpha\) to próg p-wartości poniżej którego odrzucamy hipotezę zerową. Jest to zatem (z definicji p-wartości) prawdopodobieństwo odrzucenia \(H_0\) w sytuacji gdy jest ona prawdziwa, czyli popełnienia tak zwanego błędu pierwszego rodzaju lub otrzymania wyniku fałszywie pozytywnego. Ta ostatnia nazwa pochodzi stąd, że na ogół \(H_1\) oznacza zaobserwowanie jakiegoś interesującego nas zjawiska, na przykład różnicy pomiędzy badanymi populacjami.
Moc testu, oznaczana często jako \(1 - \beta\), to prawdopodobieństwo przyjęcia \(H_1\) w stytuacji gdy jest ona prawdziwa. Sama wartość \(\beta\) jest zatem prawdopodobieństwem, że nie odrzucimy fałszywej hipotezy zerowej, czyli nie wykryjemy istniejącego efektu. Jest to tak zwany błąd drugiego rodzaju lub wynik fałszywie negatywny.

Rozwiązanie:

Na początku wygenerujemy nasze próby statystyczne.

proba1 <- replicate(1000, rnorm(10, 0.5))
proba2 <- replicate(1000, rnorm(10, -0.5))

Teraz chcemy porównać odpowiadające sobie kolumny otrzymanych macierzy za pomocą testu t Studenta.
Ponieważ wiemy, że wariancje porównywanych prób są równe, to wykorzystamy tę informację.
Jest na to kilka sposobów, na przykład połączenie obu macierzy, napisanie odpowiedniej funkcji i wykorzystanie apply().
Ja proponuję poniższy sposób:

t_p_values <- sapply(1:ncol(proba1), function(i) t.test(proba1[,i], proba2[,i], var.equal = T)$p.value)

Powyżej, po kolei, dzieją się następujące rzeczy:

W ten sam sposób przeprowadzimy test rang Wilcoxona.

w_p_values <- sapply(1:ncol(proba1), function(i) wilcox.test(proba1[,i], proba2[,i])$p.value)

Mając dwa wektory p-wartości, w pierwszym kroku najlepiej jest utworzyć histogram obrazujący z jakimi wartościami mamy do czynienia. Przedstawimy histogramy obu wektorów na jednym wykresie. Estetyka fill sprawi, że histogramy zostaną pokolorowane zgodnie z rodzajem testu.
Domyślnie biblioteka ggplot2 ustawia jeden histogram na drugim. Trzeba bardzo uważać, żeby nie pomylić się przy interpretacji takiego wykresu. Na przykład, na wykresie poniżej, pierwszy czerwony słupek jest podobnej wysokości niż pierwszy niebieski słupek, mimo, że na pierwszy rzut oka można odnieść wrażenie że jest dwa razy wyższy.

dane_do_wykresu <- data.frame('Test'=rep(c('Student', 'Wilcox'), each=1000), 'p_wartosc' = c(t_p_values, w_p_values))
ggplot(dane_do_wykresu) + geom_histogram(aes(x=p_wartosc, fill=Test), binwidth=0.01) + theme_minimal()

Z wykresu wydaje się że p-wartości z testu t są raczej niższe niż te z testu Wilcoxona. Sprawdźmy, w ilu próbach test t dał niższą p-wartość niż test Wilcoxona:

sum(t_p_values < w_p_values)
## [1] 700

W komórce powyżej w pierwszym kroku, za pomocą komendy t_p_values < w_p_values, tworzony jest wektor logiczny długości 100. Następnie zostają wysumowane jego wartości, gdzie wartość TRUE jest liczona jako 1, a FALSE jako 0. Dla 700 prób na 1000 test t dał niższą p-wartość niż test rang Wilcoxona. Czy potrafisz wyjaśnić ten fakt?

Sprawdźmy teraz moce obu testów. Na początku stworzymy wektory które powiedzą nam w których próbach każdy z testów odrzucił hipotezę \(H_0\) na poziomie istotności alpha.

alpha <- 0.05
student_odrzucil <- t_p_values < alpha
wilcox_odrzucil <- w_p_values < alpha

A teraz zobaczmy proporcję odrzuconych hipotez \(H_0\), czyli poprawnie wykrytych różnic pomiędzy grupami.

mean(student_odrzucil)
## [1] 0.58
mean(wilcox_odrzucil)
## [1] 0.529

Moc testu t Studenta na poziomie istotności alpha wyniosła 0.58, podczas gdy moc testu Wilcoxona była niższa i wyniosła 0.529. Teoretyczna moc testu t Studenta dla tych danych wynosi 0.562.

Wykorzystany przez nas test t zakłada, że nasze próby pochodzą z rozkładów normalnych o identycznych wariancjach.
Test Wilcoxona nie ma żadnych założeń co do rozkładu z których pochodzą próby.
Oba testy zakładają niezależność zmiennych losowych w obrębie prób oraz pomiędzy próbami.
Wszystkie założenia są spełnione.

Dlaczego test t ma lepszą moc? Na to pytanie można udzielić odpowiedzi na dwa sposoby.
Matematycznie ścisła odpowiedź jest taka, że tak wynika ze wzorów i już.
Odpowiedź mniej ścisła, ale za to użyteczna w praktyce, jest taka, że test t ma do dyspozycji więcej informacji na temat naszych prób. Powiedzieliśmy testowi t nie tylko to, jakie wartości mamy w danych, ale też to, jakich wartości ogólnie możemy się spodziewać.
Na ogół procedury, które mają do dyspozycji więcej informacji, zwracają bardziej precyzyjne wyniki. Oczywiście pod warunkiem, że te informacje są prawdziwe. Jeśli zmylimy nasz test i damy mu próby które nie spełniają założeń, to oczywiście dostaniemy gorsze wyniki.

Uwaga praktyczna. Teoretycznie test t wymaga prób z rozkładu normalnego. W praktyce, w większości zastosowań, jeśli próba ma dostatecznie duży rozmiar (na ogół więcej niż 20 obserwacji), to daje on dobre wyniki również dla innych rozkładów prawdopodobieństwa. Wynika to z faktu, że jedynymi statystykami z których korzysta ten test jest średnia oraz wariancja. Tak zwane centralne twierdzenie graniczne mówi, że przy dość ogólnych założeniach oba te estymatory zbiegają do rozkładu normalnego przy rosnącej liczności próby. To sprawia, że dla dostatecznie dużych prób nie ma znaczenia czy same obserwacje mają rozkład normalny, ponieważ estymatory zachowują się w przybliżeniu tak, jak gdyby to założenie było spełnione.

Zadanie 2. (podchwytliwe). W danych Zarobki.csv znajdują się miesięczne zarobki losowo wybranych obywateli krajów A i B. Wykorzystując poznane dotychczas techniki analizy danych, spróbuj odpowiedzieć na pytanie: Obywatele którego kraju są bogatsi?.

Testowanie wielokrotne

Załóżmy że testujemy 1 000 potencjalnych leków na raka.
Każdą z substancji podajemy grupie testowej, i porównujemy rozmiar guza z grupą kontrolną.
Przyjmujemy, że dana substancja nadaje się do dalszych badań, jeśli p-wartość przy odpowiednim teście jest mniejsza niż 0.05.

Nawet, jeśli żadna z badanych substancji nie ma jakiejkolwiek aktywności, to przy tym poziomie istotności dostaniemy \(1000 \cdot 0.05 = 50\) wyników pozytywnych (mówimy, że wynik jest pozytywny, jeśli odrzucamy \(H_0\), a negatywny w przeciwnym przypadku). Wynika to wprost z definicji p-wartości, a więc nie zależy od tego, jaki test przeprowadzamy.

To oznacza, że firma farmaceutyczna zmarnuje miliony dolarów na badania nad substancjami, które nie mają żadnego sensu, zamiast rozsądnie przeznaczyć te pieniądze na premie dla zarządu (w ostateczności na badanie innych substancji).
No trochę słabo.

Opracowano kilka metod kontrolowania wyników fałszywie pozytywnych. Na wykładzie poznaliśmy metodę Bonferroniego, Holma i Benjamini-Hochberga. W dalszej części skryptu zajmiemy się tą pierwszą i tą ostatnią.
Szczegóły techniczne dotyczące wszystkich metod są opisane na slajdach do wykładu.

Metoda Bonferroniego polega na otrzymaniu poprawionych p-wartości, \(p^*\), które mają tę własność, że jeśli odrzucimy \(H_0\) gdy \(p^* < \alpha\), to prawdopodobieństwo otrzymania jakiegokolwiek wyniku fałszywie pozytywnego wynosi \(\alpha\).
Jest to poprawka równie popularna, co bezużyteczna w większości zastosowań. Jest to spowodowane tym, że jest bardzo konserwatywna, czyli bardzo niechętnie odrzuca \(H_0\). Powoduje to bardzo wysoki błąd drugiego rodzaju. Na ogół możemy sobie pozwolić na pewną liczbę wyników fałszywie pozytywnych, jeśli dzięki temu uzyskamy więcej rzeczywiście interesujących wyników.

Metoda Benjamini-Hochberga kontroluje false discovery rate (FDR), czyli proporcję wyników fałszywie pozytywnych wśród wszystkich wyników pozytywnych. Polega na otrzymaniu tak zwanych q-wartości, które mają tę własność, że jeśli odrzucimy \(H_0\) gdy \(q < \alpha\), to nasze FDR wyniesie \(\alpha\). Innymi słowy, wśród wszystkich potencjalnych leków na raka wybranych do dalszych badań, proporcja tych, które nie mają aktywności, wyniesie \(\alpha\).

Bardzo dobre wprowadzenie do powyższych metod, jak również opis kilku praktycznych problemów z nimi związanych, można znaleźć tutaj.

Do przeprowadzenia korekcji p-wartości w R stosujemy funkcję p.adjust(). Sposób jej zastosowania jest opisany w dokumentacji.

Zadanie 3. Zainstaluj i załaduj bibliotekę ISLR. Załaduj dane Khan za pomocą komendy data(Khan) (zwróć uwagę na brak cudzysłowu).
Dane dotyczą ekspresji genów w zależności od rodzaju nowotworu. Wyniki pomiaru ekspresji genów zapisane są w polu Khan$xtrain, a rodzaj nowotworu w polu Khan$ytrain. Więcej informacji na temat danych możesz znaleźć w dokumentacji, którą możesz wyświetlić w Rstudio wpisując w konsolę ?Khan lub help(Khan).

Zadanie 4. Jaką metodę korekcji p-wartości (lub jej brak) zastosowalibyśmy do poniższych zadań badawczych?

  1. Chcemy znaleźć geny, które są potencjalnym celem terapii genowej, i mamy środki na kilkadziesiąt dalszych eksperymentów
  2. Chcemy znaleźć geny, które mogą być potencjalnym celem niesamowicie drogiego eksperymentu mającego na celu przetestowanie nowej metody leczenia
  3. Chcemy wstępnie określić, które geny mogą mieć związek z danym typem nowotworu
  4. Chcemy ogólnie scharakteryzować, które geny mają różną ekspresję w różnych typach nowotworu
  5. Chcemy znaleźć geny markerowe, czyli takie, które pozwolą na efektywne rozróznienie typu nowotworu

W następnym zadaniu wykorzystamy kilka miar jakości procedur testowania.
Miary te są w powszechnym użyciu zarówno w statystyce, jak i uczeniu maszynowym.
Spotkamy się z nimi jeszcze wielokrotnie, między innymi na zajęciach dotyczących klasyfikatorów.
Ze względu na brak dobrych polskich odpowiedników, poniżej w większości przypadków stosuję nazwy angielskie.

Do opisu miar jakości testów wykorzystamy terminologię zapożyczoną z uczenia maszynowego, a podsumowaną w poniższej tabeli:

\(H_0\) fałszywa \(H_0\) prawdziwa
Odrzucamy \(H_0\) True Positive, TP False Positive, FP
Nie odrzucamy \(H_0\) False Negative, FN True Negative, TN

Miary, które wykorzystamy, to:

Większość z powyższych miar ma kilka nazw, które są stosowane w różnych dziedzinach nauki.
Więcej na ich temat zostanie powiedziane na jednym z następnych wykładów.
Osoby zainteresowane mogą przeczytać więcej na ten temat tutaj.

Zadanie 5. Stwórz macierz o wymiarach 10x1000, taką, że obserwacje z pierwszych stu kolumn są wylosowane z rozkładu \(\mathcal{N}(1, 1)\), a z pozostałych 900 z rozkładu \(\mathcal{N}(0, 1)\). Następnie: