Ciągłe uogólnianie ujemnego rozkładu dwumianowego

Ujemny rozkład dwumianowy (NB) jest zdefiniowany na nieujemnych liczbach całkowitych i ma funkcję masy prawdopodobieństwa

f (k; r, p) = (\binom{k + r - 1}{k}) p^{k} (1 - p)^{r} .

$f(k;r,p)={\binom {k+r-1}{k}}p^{k}(1-p)^{r}.$ Czy ma sens rozważenie ciągłego rozkładu na

rzeczywistych nieujemnych zdefiniowanych przez tę samą formułę (zastępując

k \in N_{0}

$k\in \mathbb N_0$ przez

x \in R_{\geq 0}

$x\in\mathbb R_{\ge 0}$ )? Współczynnik dwumianowy można przepisać jako iloczyn

(k + 1) \cdot \dots \cdot (k + r - 1)

$(k+1)\cdot\ldots\cdot(k+r-1)$ , który jest dobrze zdefiniowany dla dowolnego rzeczywistego

k

$k$ . Mamy więc plik PDF

f (x; r, p) \propto \prod_{i = 1}^{r - 1} (x + i) \cdot p^{x} (1 - p)^{r} .

$f(x;r,p)\propto\prod_{i=1}^{r-1}(x+i)\cdot p^{x}(1-p)^{r}.$ Mówiąc bardziej ogólnie, możemy zastąpić współczynnik dwumianowy funkcjami Gamma, pozwalając na wartości nie całkowite

r

$r$ :

f (x; r, p) \propto \frac{Γ (x + r)}{Γ (x + 1) Γ (r)} \cdot p^{x} (1 - p)^{r} .

$f(x;r,p)\propto\frac{\Gamma(x+r)}{\Gamma(x+1)\Gamma(r)}\cdot p^{x}(1-p)^{r}.$

Czy to jest prawidłowa dystrybucja? Czy to ma imię? Czy to ma jakieś zastosowania? Czy to może jakiś związek lub mieszanina? Czy istnieją zamknięte wzory na średnią i wariancję (i stałą proporcjonalności w pliku PDF)?

(Obecnie studiuję artykuł, który wykorzystuje model mieszanki NB (ze stałym $r=2$ ) i pasuje do niego za pomocą EM. Jednak dane są liczbami całkowitymi po pewnej normalizacji, tj. Nie liczbami całkowitymi. Niemniej autorzy stosują standardową formułę NB do obliczenia prawdopodobieństwo i uzyskać bardzo rozsądne wyniki, więc wszystko wydaje się działać dobrze. Uważam, że to bardzo zagadkowe. Pamiętaj, że to pytanie nie dotyczy NB GLM).

— ameba mówi Przywróć Monikę
źródło

Czy nie byłoby to połączenie gamma z parametrem skali

- \log p

$-\log p$ ? Jeśli rozwiniesz wielomian

Π_{i = 1}^{r - 1} (x + i)

$\Pi_{i=1}^{r-1}(x+i)$ , otrzymasz po prostu

\sum_{i = 2}^{r} a_{i} x^{i - 1}

$\sum_{i=2}^ra_ix^{i-1}$ , a następnie pomnożymy przez

p^{x}

$p^x$ jest taki sam jak

\exp {x \log p}

$\exp\{x\log p\}$ , gdzie

a_{i}

$a_i$ jest współczynnikiem

x^{i - 1}

$x^{i-1}$ w wielomianu i

\log p < 0

$\log p < 0$ oczywiście, więc wygląda na to, że przekonwertowałby na średnia ważona rozkładów gamma, tj. mieszaniny.

— jbowman

... właściwie powinno być

i = 1

$i=1$ w powyższej sumie.

— jbowman

Ponieważ zależy tylko od parametrów, jest to stała, którą można zaabsorbować proporcjonalnością. Ponadto również ma stałą która może być ignorowanym. Pisząc dla , pytasz o gęstość proporcjonalną doTo identyfikuje jako współczynnik skali, a jako parametr kształtu. W przypadku całki jest to wyraźnie mieszanina rozkładów gamma. Jednak nie ma sensu ograniczać do liczb całkowitych.

(1 - p)^{r}

$(1-p)^r$

(\binom{x + r - 1}{x}) = Γ (x + r) / (Γ (r) Γ (x + 1))

$\binom{x+r-1}{x}=\Gamma(x+r)/(\Gamma(r)\Gamma(x+1))$

1 / Γ (r)

$1/\Gamma(r)$

p^{k} = e^{- k ρ}

$p^k=e^{-k\rho}$

ρ = - \log (p) \geq 0

$\rho=-\log(p) \ge 0$

f (x; r, ρ) = \frac{Γ (x + r)}{Γ (x + 1)} e^{- ρ x} .

$f(x;r,\rho)=\frac{\Gamma(x+r)}{\Gamma(x+1)}\,e^{-\rho x}.$

ρ

$\rho$

r

$r$

r

$r$

r

$r$

— whuber

@whuber Right. Właściwie używam rozkładu, który jest ciągły na wartościach dodatnich i ma masę punktową równą zero. Uważam, że jest to właściwe podejście. Ale zasugerowano mi, aby stosować ciągłe uogólnianie NB, które miałoby niezerowe prawdopodobieństwo przy zera, a zatem pozornie pozwalało poradzić sobie z dokładnymi zerami. Stąd moje pytanie.

— ameba mówi Przywróć Monikę

Wydaje mi się, że w tej sugestii może być pewne zamieszanie: wydaje się, że łączy prawdopodobieństwo (czyli to, co ma masa punktowa lub rozkład NB na zero) z gęstością prawdopodobieństwa (która jest wartością byłoby). Niezerowa gęstość nie pozwala ci radzić sobie z dokładnymi zerami, ponieważ nadal przewiduje zerową szansę, że pojawi się dowolna wartość !

f (0, θ)

$f(0,\theta)$

0

$0$

— whuber

Odpowiedzi:

To interesujące pytanie. Moja grupa badawcza korzysta z dystrybucji, o której mówisz od kilku lat, w naszym publicznie dostępnym oprogramowaniu bioinformatycznym. O ile mi wiadomo, dystrybucja nie ma nazwy i nie ma na jej temat literatury. Chociaż artykuł Chandra i wsp. (2012) cytowany przez Aksakala jest ściśle powiązany, ich rozkład wydaje się ograniczony do wartości całkowitych dla i nie wydaje się, aby zawierały wyraźne wyrażenie dla pliku pdf. $r$

Aby dać ci trochę tła, rozkład NB jest bardzo intensywnie wykorzystywany w badaniach genomowych do modelowania danych dotyczących ekspresji genów wynikających z sekwencji RNA i powiązanych technologii. Dane zliczania powstają, gdy liczba odczytanych sekwencji DNA lub RNA wyekstrahowanych z próbki biologicznej, którą można zmapować do każdego genu. Zazwyczaj z każdej próbki biologicznej są dziesiątki milionów odczytów zmapowanych do około 25 000 genów. Alternatywnie można mieć próbki DNA, z których odczyty są mapowane na okna genomowe. My i inni spopularyzowaliśmy podejście, w którym NB glms są dopasowywane do odczytów sekwencji dla każdego genu, a empiryczne metody Bayesa są wykorzystywane do moderowania genówowych estymatorów dyspersji (dyspersja $\phi=1/r$ ). Takie podejście zostało przytoczone w dziesiątkach tysięcy artykułów w czasopiśmie w literaturze genomicznej, dzięki czemu można zorientować się, jak bardzo się przyzwyczai.

Moja grupa utrzymuje pakiet oprogramowania edgeR R.. Kilka lat temu zmieniliśmy cały pakiet, aby działał z licznikami ułamkowymi, używając ciągłej wersji NB pmf. Po prostu przekonwertowaliśmy wszystkie współczynniki dwumianowe w NB pmf na stosunki funkcji gamma i zastosowaliśmy je jako (mieszany) ciągły plik pdf. Motywacją tego było to, że zliczanie odczytów sekwencji może czasami być ułamkowe z powodu (1) niejednoznacznego mapowania odczytów na transkryptom lub genom i / lub (2) normalizacji zliczeń w celu skorygowania efektów technicznych. Tak więc liczby są czasami oczekiwanymi lub szacowanymi, a nie obserwowanymi. I oczywiście liczba odczytów może wynosić dokładnie zero z prawdopodobieństwem dodatnim. Nasze podejście zapewnia, że wyniki wnioskowania z naszego oprogramowania są ciągłe w zliczeniach, dokładnie dopasowane do dyskretnych wyników NB, gdy szacowane liczby są liczbami całkowitymi.

O ile mi wiadomo, nie ma zamkniętej formy stałej normalizującej w pliku pdf, ani też nie ma zamkniętych postaci średniej lub wariancji. Gdy ktoś uważa, że nie ma formy zamkniętej dla całki (stała Fransena-Robinsona), jasne jest, że nie może być całki ciągłej NB pdf albo. Wydaje mi się jednak, że tradycyjne wzory średnich i wariancji dla NB powinny nadal być dobrym przybliżeniem dla ciągłej NB. Ponadto stała normalizująca powinna zmieniać się powoli wraz z parametrami, a zatem może być ignorowana jako mająca znikomy wpływ na obliczenia maksymalnego prawdopodobieństwa.

\int_{0}^{\infty} \frac{1}{Γ (x)} d z

$\int_0^\infty \frac{1}{\Gamma(x)}dz$

Można potwierdzić te hipotezy poprzez całkowanie numeryczne. Rozkład NB powstaje w bioinformatyce jako mieszanina gamma rozkładów Poissona (patrz Wikipedia dwumianowy artykuł negatywny lub McCarthy et al poniżej). Ciągły rozkład NB powstaje po prostu przez zastąpienie rozkładu Poissona ciągłym analogiem pdf dla gdzie jest stałą normalizującą, aby zapewnić integrację gęstości z 1. Załóżmy na przykład, że . Rozkład Poissona ma pmf równe powyższemu pdf na liczbach całkowitych nieujemnych, a

f (x; λ) = a (λ) \frac{e^{- λ} λ^{x}}{Γ (x + 1)}

$f(x;\lambda)=a(\lambda)\frac{e^{-\lambda}\lambda^x}{\Gamma(x+1)}$

x \geq 0

$x\ge 0$

a (λ)

$a(\lambda)$

λ = 10

$\lambda=10$

λ = 10

$\lambda=10$ , średnia Poissona i wariancja są równe 10. Całkowanie numeryczne pokazuje, że a średnia i wariancja ciągłego rozkładu są równe 10 do około 4 cyfr znaczących. Zatem stała normalizująca wynosi praktycznie 1, a średnia i wariancja są prawie dokładnie takie same, jak w przypadku dyskretnego rozkładu Poissona. Przybliżenie poprawia się jeszcze bardziej, jeśli dodamy korektę ciągłości, integrując od do zamiast od 0. Przy korekcji ciągłości wszystko jest poprawne (stała normalizująca wynosi 1, a momenty zgadzają się z dyskretnym Poissonem) do około 6 liczby.

a (10) = 1 / 0.999875

$a(10)=1/0.999875$

- 1 / 2

$-1/2$

\infty

$\infty$

W naszym pakiecie edgeR nie musimy dokonywać żadnych korekt ze względu na fakt, że masa jest równa zeru, ponieważ zawsze pracujemy z warunkowymi prawdopodobieństwami logarytmicznymi lub z różnicami wiarygodności logarytmicznymi, a wszelkie funkcje delta anulują obliczenia. Jest to typowa wartość BTW dla mieszanych rozkładów prawdopodobieństwa. Alternatywnie, moglibyśmy rozważyć rozkład, który nie ma masy w punkcie zerowym, ale ma wsparcie rozpoczynające się od -1/2 zamiast od zera. Każda z perspektyw teoretycznych prowadzi do tych samych obliczeń w praktyce.

Chociaż aktywnie korzystamy z ciągłej dystrybucji NB, nie opublikowaliśmy niczego na jej temat. Artykuły cytowane poniżej wyjaśniają podejście NB do danych genomowych, ale nie omawiają wyraźnie ciągłego rozkładu NB.

Podsumowując, nie dziwię się, że artykuł, który studiujesz, uzyskał rozsądne wyniki z kontynuowanej wersji NB pdf, ponieważ takie jest również nasze doświadczenie. Kluczowym wymaganiem jest prawidłowe modelowanie średnich i wariancji, co będzie w porządku, pod warunkiem, że dane, zarówno całkowite, jak i nie, wykazują tę samą formę kwadratowej zależności średniej wariancji, co rozkład NB.

Referencje

Robinson, M., i Smyth, GK (2008). Szacowanie małej próby ujemnej dyspersji dwumianowej, z zastosowaniem danych SAGE . Biostatistics 9, 321-332.

Robinson, MD i Smyth, GK (2007). Moderowane testy statystyczne do oceny różnic w liczności znaczników . Bioinformatics 23, 2881-2887.

McCarthy, DJ, Chen, Y, Smyth, GK (2012). Analiza różnicowa ekspresji eksperymentów wieloczynnikowych RNA-Seq w odniesieniu do zmienności biologicznej . Nucleic Acids Research 40, 4288-4297.

Chen, Y, Lun, ATL i Smyth, GK (2014). Analiza różnicowa ekspresji złożonych eksperymentów z sekwencją RNA z wykorzystaniem edgeR. W: Analiza statystyczna danych sekwencji następnej generacji, Somnath Datta i Daniel S Nettleton (red.), Springer, Nowy Jork, strony 51--74. Przedruk

Lun, ATL, Chen, Y i Smyth, GK (2016). Jest DE-licious: przepis na analizy ekspresji różnicowej eksperymentów z sekwencją RNA z wykorzystaniem metod quasi-prawdopodobieństwa w EdgeR. Methods in Molecular Biology 1418, 391-416. Przedruk

Chen Y, Lun ATL i Smyth, GK (2016). Od odczytów przez geny do ścieżek: analiza ekspresji różnicowej eksperymentów RNA-Seq przy użyciu Rsubread i potoku quasi-prawdopodobieństwa edgeR . F1000 Badanie 5, 1438.

— Gordon Smyth
źródło

Jest to niezwykle pomocne, @Gordon; wielkie dzięki za poświęcenie czasu na napisanie tego. Pracuję również z danymi o sekwencji RNA, więc odpowiedź z tej perspektywy jest szczególnie cenna (dodałem teraz do pytania znacznik [bioinformatyka]). Wasza praca dotyczy ekspresji różnicowej, podczas gdy moja obecna praca dotyczy grupowania (artykuł, który czytałem to Harris i wsp. Na temat interneuronów CA1; biorxiv ). W każdym razie pozwól, że zadam ci kilka drobnych pytań / wyjaśnień. [cd.]

— ameba mówi Przywróć Monikę

(1) Powiedziałeś, że ciągła NB jest mieszanką gamma ciągłych Poissonów. Czy możesz go trochę rozwinąć, a może pokazać to bardziej jednoznacznie? Myślę, że będzie to przydatne dla ogółu odbiorców. W związku z tym w komentarzach do mojego pytania dwie osoby napisały, że ciągła NB powinna być mieszanką gamma z parametrem skali , ale tylko dla liczby całkowitej . Czy oba poglądy są prawdziwe? (2) Powiedziałeś, że funkcja delta na zero nie ma znaczenia dla GLM. Jednocześnie istnieje duża literatura na temat GLM o zerowych napompowaniach. Jak to do siebie pasuje?

- \log (p)

$-\log(p)$

r

$r$

— ameba mówi Przywróć Monikę

(3) W swojej praktycznej pracy, należy użyć ML oszacować wszystkie parametry, w tym , czy też naprawić do pewnej wartości określonej z góry (być może ta sama wartość wspólna dla wszystkich genów?), A potem trzymać go na stałym poziomie? Sądzę, że powinno to być znacznie łatwiejsze. (Np. Sama NB jest rodziną dyspersyjną wykładniczą, ale tylko ze stałym .)

r

$r$

r

$r$

r

$r$

— mówi ameba Przywróć Monikę

@amoeba Dzięki za biorxiv ref. (1) Pochodzenie NB jako mieszaniny Poissons jest dość dobrze znane i znajduje się w naszych pracach, np. McCarthy i in. Wyprowadzenie ciągłego NB następuje po prostu przez zastąpienie ciągłego Poissona Poissonem. Czy powinienem to dodać do mojej odpowiedzi? Długo by to trwało. Nie rozumiem, w jaki sposób ciągła NB mogłaby być użytecznie reprezentowana jako mieszanka gamma. (2) Nie, inflacja zerowa jest inną dodatkową komplikacją. Unikamy tych komplikacji w naszej pracy.

— Gordon Smyth,

@amoeba (3) Szacujemy wszystkie parametry. Aby uzyskać kontrolę poziomu błędu, należy oszacować rozproszenie genów, a należy to zrobić ze szczególną ostrożnością, ponieważ rozmiary próbek są często małe, a rozmiar danych jest ogromny. Stosujemy złożoną procedurę, która obejmuje skorygowane prawdopodobieństwo profilu (pomyśl REML) w obrębie każdego genu połączonego z empiryczną procedurą Bayesa między genami o ważonym prawdopodobieństwie. Genomowe NB glms są następnie dopasowywane przez ML z ustalonymi dyspersjami. Na koniec współczynniki są testowane przy użyciu testów F quasi-prawdopodobieństwa.

— Gordon Smyth,

Spójrz na ten artykuł: Chandra, Nimai Kumar i Dilip Roy. Ciągła wersja ujemnego rozkładu dwumianowego. Statistica 72, no. 1 (2012): 81 .

W artykule zdefiniowano ją jako funkcję przeżycia, co jest naturalnym podejściem od czasu wprowadzenia dwumianu ujemnego w analizie niezawodności:

S_{r} (x) = {\begin{cases} q^{x} & for r = 1 \\ \sum_{k = 0}^{r - 1} (\binom{x + k - 1}{k}) p^{k} q^{x} & for r = 2, 3, \dots \end{cases}

$S_r(x)=\begin{cases}q^x & \text{for}\ r=1 \\ \sum_{k=0}^{r-1}\binom {x+k-1}{k}p^kq^x & \text{for}\ r=2,3,\dots \end{cases}$ gdzie i .

q = e^{- λ}, λ \geq 0, p + q = 1

$q=e^{-\lambda},\lambda\ge 0,p+q=1$

r \in N, r > 0

$r\in\mathbb N,r>0$

— Aksakal
źródło

Dzięki! Rzucę okiem na ten artykuł. (To nie ja przegłosowałem.)

— Ameba mówi Przywróć Monikę

@amoeba, nie martwię się o przegłosowanie, to internet :)

— Aksakal

(Dziwne, że odpowiedź została odrzucona ...) +1

— whuber

Dobrze jest mieć to odniesienie, ale idealnie chciałbym zobaczyć tutaj bardziej szczegółową dyskusję. Czy ta funkcja przeżycia definiuje taki sam rozkład jak plik PDF w moim pytaniu? (Nawiasem mówiąc, wydaje mi się to trochę dziwne, że autorzy używają współczynników dwumianowych dla niecałkowitych wartości .) Kilka powyższych uwag wskazuje, że jest to mieszanina rozkładów gamma (nie widzę w tym żadnej dyskusji papier); jakie są parametry tych gamma, jakie są masy mieszanin? Czy formuły NB dla średniej i wariancji dotyczą wersji ciągłej?

x

$x$

— ameba mówi Przywróć Monikę

@amoeba, gazeta ma chwile, nie są takie same jak w NB, niestety

— Aksakal