Dlaczego i jak dostosować wartości P w testowaniu wielu hipotez
Testowanie wielu hipotez ma miejsce, gdy wielokrotnie testujemy modele na wielu cechach, ponieważ prawdopodobieństwo uzyskania jednego lub więcej fałszywych odkryć rośnie wraz z liczbą testów. Na przykład w dziedzinie genomiki naukowcy często chcą sprawdzić, czy którykolwiek z tysięcy genów ma znacząco różną aktywność w wyniku będącym przedmiotem zainteresowania. Lub czy żelki powodują trądzik .
W tym poście na blogu omówimy kilka popularnych metod używanych do wyjaśniania testowania wielu hipotez poprzez dostosowanie wartości p modelu:
- Odsetek wyników fałszywie dodatnich (FPR)
- Wskaźnik błędów rodzinnych (FWER)
- Współczynnik fałszywych odkryć (FDR)
Ten dokument można podsumować następującym obrazem:
Utwórz dane testowe
Stworzymy symulowany przykład, aby lepiej zrozumieć, w jaki sposób różne manipulacje wartościami p mogą prowadzić do różnych wniosków. Aby uruchomić ten kod, potrzebujemy Pythona z zainstalowanymi bibliotekami pandas, numpy, scipyi .statsmodels
Na potrzeby tego przykładu zaczynamy od utworzenia Pandas DataFrame zawierającej 1000 funkcji. 990 z nich (99%) będzie miało swoje wartości wygenerowane z rozkładu normalnego ze średnią = 0, zwanego modelem zerowym. (W funkcji norm.rvs()użytej poniżej średnia jest ustawiana za pomocą locargumentu.) Pozostały 1% cech zostanie wygenerowany na podstawie rozkładu normalnego = 3, zwanego modelem Non-Null. Użyjemy ich jako reprezentujących interesujące funkcje, które chcielibyśmy odkryć.
import pandas as pd
import numpy as np
from scipy.stats import norm
from statsmodels.stats.multitest import multipletests
np.random.seed(42)
n_null = 9900
n_nonnull = 100
df = pd.DataFrame({
'hypothesis': np.concatenate((
['null'] * n_null,
['non-null'] * n_nonnull,
)),
'feature': range(n_null + n_nonnull),
'x': np.concatenate((
norm.rvs(loc=0, scale=1, size=n_null),
norm.rvs(loc=3, scale=1, size=n_nonnull),
))
})
Wartości P można obliczyć z skumulowanego rozkładu ( norm.cdf()z scipy.stats), który reprezentuje prawdopodobieństwo uzyskania wartości równej lub mniejszej niż obserwowana. Następnie, aby obliczyć wartość p, obliczamy 1 - norm.cdf(), aby znaleźć prawdopodobieństwo większe niż obserwowane:
df['p_value'] = 1 - norm.cdf(df['x'], loc = 0, scale = 1)
df
Pierwsza koncepcja nosi nazwę współczynnika wyników fałszywie dodatnich i jest zdefiniowana jako ułamek hipotez zerowych, które oznaczamy jako „istotne” (zwane również błędami typu I). Wartości p, które obliczyliśmy wcześniej, można interpretować jako odsetek wyników fałszywie dodatnich z samej ich definicji: są to prawdopodobieństwa uzyskania wartości co najmniej tak dużej jak określona wartość, gdy próbkujemy rozkład zerowy.
Dla celów ilustracyjnych zastosujemy wspólny (magiczny ) próg wartości p wynoszący 0,05, ale można zastosować dowolny próg:
df['is_raw_p_value_significant'] = df['p_value'] <= 0.05
df.groupby(['hypothesis', 'is_raw_p_value_significant']).size()
hypothesis is_raw_p_value_significant
non-null False 8
True 92
null False 9407
True 493
dtype: int64
Główny problem z FPR polega na tym, że w rzeczywistym scenariuszu nie wiemy a priori, które hipotezy są zerowe, a które nie. Wtedy sama surowa wartość p (odsetek wyników fałszywie dodatnich) ma ograniczone zastosowanie. W naszym przypadku, gdy ułamek cech niepustych jest bardzo mały, większość cech oznaczonych jako znaczące będzie zerowa, ponieważ jest ich znacznie więcej. Konkretnie, z 92 + 493 = 585 cech oznaczonych jako prawda („pozytywne”), tylko 92 pochodzi z naszej dystrybucji innej niż zerowa. Oznacza to, że większość lub około 84% zgłoszonych znaczących cech (493/585) to fałszywe alarmy!
Więc, co możemy z tym zrobić? Istnieją dwie popularne metody rozwiązania tego problemu: zamiast współczynnika fałszywych trafień możemy obliczyć współczynnik błędów rodzinnych (FWER) lub współczynnik fałszywych odkryć (FDR). Każda z tych metod przyjmuje zestaw nieprzetworzonych, nieskorygowanych wartości p jako danych wejściowych i tworzy nowy zestaw „skorygowanych wartości p” jako dane wyjściowe. Te „skorygowane wartości p” reprezentują szacunki górnych granic FWER i FDR. Można je uzyskać z multipletests()funkcji, która jest częścią statsmodelsbiblioteki Pythona:
def adjust_pvalues(p_values, method):
return multipletests(p_values, method = method)[1]
Współczynnik błędów rodzinnych to prawdopodobieństwo fałszywego odrzucenia jednej lub więcej hipotez zerowych lub innymi słowy oznaczanie prawdziwej wartości zerowej jako niezerowej lub prawdopodobieństwo napotkania jednego lub więcej wyników fałszywie dodatnich.
Gdy testowana jest tylko jedna hipoteza, jest ona równa surowej wartości p (odsetek wyników fałszywie dodatnich). Jednak im więcej hipotez zostanie przetestowanych, tym bardziej prawdopodobne jest, że otrzymamy jeden lub więcej fałszywych alarmów. Istnieją dwa popularne sposoby szacowania FWER: procedury Bonferroniego i Holma. Chociaż ani procedury Bonferroniego, ani Holma nie zakładają zależności przeprowadzanych testów od poszczególnych cech, będą one zbyt konserwatywne. Na przykład w skrajnym przypadku, gdy wszystkie cechy są identyczne (ten sam model powtórzony 10 000 razy), korekta nie jest potrzebna. Będąc na drugim biegunie, gdzie żadne cechy nie są ze sobą skorelowane, wymagany jest pewien rodzaj korekty.
Procedura Bonferroniego
Jedną z najpopularniejszych metod korygowania przy testowaniu wielu hipotez jest procedura Bonferroniego. Powodem, dla którego ta metoda jest popularna, jest to, że jest bardzo łatwa do obliczenia, nawet ręcznie. Ta procedura mnoży każdą wartość p przez całkowitą liczbę przeprowadzonych testów lub ustawia ją na 1, jeśli to mnożenie spowodowałoby przekroczenie 1.
df['p_value_bonf'] = adjust_pvalues(df['p_value'], 'bonferroni')
df.sort_values('p_value_bonf')
Procedura Holma zapewnia korektę, która jest silniejsza niż procedura Bonferroniego. Jedyna różnica polega na tym, że nie wszystkie wartości p są mnożone przez całkowitą liczbę testów (tutaj 10000). Zamiast tego każda posortowana wartość p jest mnożona progresywnie przez malejącą sekwencję 10000, 9999, 9998, 9997, …, 3, 2, 1.
df['p_value_holm'] = adjust_pvalues(df['p_value'], 'holm')
df.sort_values('p_value_holm').head(10)
Jeśli ponownie zastosujemy nasz próg wartości p równy 0,05, przyjrzyjmy się, jak te skorygowane wartości p wpływają na nasze prognozy:
df['is_p_value_holm_significant'] = df['p_value_holm'] <= 0.05
df.groupby(['hypothesis', 'is_p_value_holm_significant']).size()
hypothesis is_p_value_holm_significant
non-null False 92
True 8
null False 9900
dtype: int64
Jednak to podejście ma wadę: nie oznaczono innych 92 niepustych funkcji jako istotnych. Chociaż bardzo rygorystycznie starano się upewnić, że żadna z pustych funkcji nie wślizgnęła się, udało się znaleźć tylko 8% (8 ze 100) niepustych funkcji. Można to postrzegać jako przyjęcie innej skrajności niż podejście fałszywie dodatnie.
Czy istnieje bardziej pośrednie podłoże? Odpowiedź brzmi „tak”, a tym środkiem jest współczynnik fałszywych odkryć.
Wskaźnik fałszywych odkryć
A co, jeśli nie mamy nic przeciwko dopuszczaniu niektórych fałszywych trafień, ale przechwytywaniu więcej niż jednocyfrowego procentu prawdziwych trafień? Być może nie mamy nic przeciwko temu, by mieć kilka fałszywych trafień, ale nie tak wiele, aby przytłoczyły one wszystkie funkcje, które oznaczyliśmy jako znaczące — jak miało to miejsce w przykładzie z FPR.
Można to zrobić, kontrolując wskaźnik fałszywego wykrywania (zamiast FWER lub FPR) na określonym poziomie progowym, powiedzmy 0,05. Współczynnik fałszywych odkryć to ułamek fałszywych trafień spośród wszystkich cech oznaczonych jako pozytywne: FDR = FP / (FP + TP), gdzie FP to liczba fałszywych trafień, a TP to liczba prawdziwych trafień. Ustawiając próg FDR na 0,05, mówimy, że jesteśmy w porządku z 5% (średnio) fałszywych alarmów wśród wszystkich naszych funkcji, które oznaczyliśmy jako pozytywne.
Istnieje kilka metod kontrolowania FDR i tutaj opiszemy, jak wykorzystać dwie popularne: procedury Benjaminiego-Hochberga i Benjaminiego-Yekutieli. Obie te procedury są podobne, chociaż bardziej zaangażowane niż procedury FWER. Nadal polegają na sortowaniu wartości p, mnożeniu ich przez określoną liczbę, a następnie zastosowaniu kryterium odcięcia.
Procedura Benjaminiego-Hochberga
Procedura Benjaminiego-Hochberga (BH) zakłada, że każdy z testów jest niezależny . Testy zależne występują na przykład wtedy, gdy testowane cechy są ze sobą skorelowane. Obliczmy wartości p skorygowane o BH i porównajmy je z naszym wcześniejszym wynikiem z FWER przy użyciu poprawki Holma:
df['p_value_bh'] = adjust_pvalues(df['p_value'], 'fdr_bh')
df[['hypothesis', 'feature', 'x', 'p_value', 'p_value_holm', 'p_value_bh']] \
.sort_values('p_value_bh') \
.head(10)
df['is_p_value_holm_significant'] = df['p_value_holm'] <= 0.05
df.groupby(['hypothesis', 'is_p_value_holm_significant']).size()
hypothesis is_p_value_holm_significant
non-null False 92
True 8
null False 9900
dtype: int64
df['is_p_value_bh_significant'] = df['p_value_bh'] <= 0.05
df.groupby(['hypothesis', 'is_p_value_bh_significant']).size()
hypothesis is_p_value_bh_significant
non-null False 67
True 33
null False 9898
True 2
dtype: int64
Zauważ, że w tym przypadku mamy stawkę FDR 6%, mimo że naszym celem było kontrolowanie jej na poziomie 5%. FDR będzie kontrolowany średnio na poziomie 5% : czasem może być niższy, a czasem wyższy.
Procedura Benjaminiego-Yekutielego
Procedura Benjamini-Yekutieli (BY) kontroluje FDR niezależnie od tego, czy testy są niezależne, czy nie. Ponownie warto zauważyć, że wszystkie te procedury mają na celu ustalenie górnych granic FDR (lub FWER), więc mogą być mniej lub bardziej konserwatywne. Porównajmy procedurę BY z powyższymi procedurami BH i Holm:
df['p_value_by'] = adjust_pvalues(df['p_value'], 'fdr_by')
df[['hypothesis', 'feature', 'x', 'p_value', 'p_value_holm', 'p_value_bh', 'p_value_by']] \
.sort_values('p_value_by') \
.head(10)
df['is_p_value_by_significant'] = df['p_value_by'] <= 0.05
df.groupby(['hypothesis', 'is_p_value_by_significant']).size()
hypothesis is_p_value_by_significant
non-null False 93
True 7
null False 9900
dtype: int64
Streszczenie
Ostatecznie wybór procedury pozostawia się użytkownikowi i zależy od tego, co próbuje zrobić analiza. Cytując Benjaminiego, Hochberg (Royal Stat. Soc. 1995):
Często kontrola FWER nie jest do końca potrzebna. Kontrola FWER jest ważna, gdy wniosek z różnych indywidualnych wniosków może być błędny, gdy przynajmniej jeden z nich jest błędny.
Może tak być na przykład w przypadku, gdy kilka nowych zabiegów konkuruje ze standardem, a jeden zabieg jest wybierany z zestawu zabiegów, które są deklarowane jako znacznie lepsze od standardu.
W innych przypadkach, w których możemy mieć pewne wyniki fałszywie dodatnie, metody FDR, takie jak korekcja BH, zapewniają mniej rygorystyczne korekty wartości p i mogą być preferowane, jeśli przede wszystkim chcemy zwiększyć liczbę wyników prawdziwie dodatnich, które przekraczają określoną wartość p próg.
Istnieją inne metody regulacji, o których tutaj nie wspomniano, w szczególności wartość q , która jest również używana do sterowania FDR iw chwili pisania tego tekstu istnieje tylko jako pakiet R.
![Czym w ogóle jest lista połączona? [Część 1]](https://post.nghiatu.com/assets/images/m/max/724/1*Xokk6XOjWyIGCBujkJsCzQ.jpeg)
