Generalización continua de la distribución binomial negativa.

La distribución binomial negativa (NB) se define en enteros no negativos y tiene una función de masa de probabilidad

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}.$ ¿Tiene sentido considerar una distribución continua en reales no negativos definidos por la misma fórmula (reemplazando

k \in N_{0}

$k\in \mathbb N_0$ por

x \in R_{\geq 0}

$x\in\mathbb R_{\ge 0}$ )? El coeficiente binomial se puede reescribir como un producto de

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

$(k+1)\cdot\ldots\cdot(k+r-1)$ , que está bien definido para cualquier

real

k

$k$ . Entonces tendríamos un 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}.$ De manera más general, podemos reemplazar el coeficiente binomial con funciones Gamma, lo que permite valores no enteros de

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}.$

¿Es una distribución válida? Eso tiene un nombre? ¿Tiene algún uso? ¿Es quizás algún compuesto o una mezcla? ¿Existen fórmulas cerradas para la media y la varianza (y la constante de proporcionalidad en el PDF)?

(Actualmente estoy estudiando un artículo que utiliza el modelo de mezcla NB (con fijo $r=2$ ) y lo ajusta a través de EM. Sin embargo, los datos son enteros después de cierta normalización, es decir, no enteros. Sin embargo, los autores aplican la fórmula estándar NB para calcular la probabilidad y obtener resultados muy razonables, por lo que todo parece funcionar bien. Me pareció muy desconcertante. Tenga en cuenta que esta pregunta no se trata de NB GLM).

distributions negative-binomial count-data continuous-data bioinformatics ameba dice Reinstate Monica
fuente

¿No sería una mezcla de Gammas con parámetro de escala

- \log p

$-\log p$ ? Si expande el polinomio

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

$\Pi_{i=1}^{r-1}(x+i)$ obtendrá

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

$\sum_{i=2}^ra_ix^{i-1}$ , y luego multiplique por

p^{x}

$p^x$ es lo mismo que by

\exp {x \log p}

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

a_{i}

$a_i$ es el coeficiente de

x^{i - 1}

$x^{i-1}$ en el polinomio y

\log p < 0

$\log p < 0$ , por supuesto, por lo que parece que se convertiría en un promedio ponderado de distribuciones Gamma, es decir, una mezcla.

jbowman

... debería ser

i = 1

$i=1$ en la suma anterior, en realidad.

jbowman

Como depende solo de los parámetros, es una constante que puede absorberse en la proporcionalidad. Además, también tiene una constante que puede ser ignorado. Al escribir para , está preguntando acerca de una densidad proporcional aEso identifica a como un factor de escala y como un parámetro de forma. Para integral es claramente una mezcla de distribuciones gamma. Sin embargo, no tiene sentido restringir a enteros.

(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 Derecha. De hecho, estoy usando una distribución que es continua en valores positivos y tiene un punto de masa en cero. Creo que este es el enfoque correcto. Pero me han sugerido usar una generalización continua de NB que tendría una probabilidad distinta de cero en cero y, por lo tanto, aparentemente permitiría tratar con ceros exactos. De ahí mi pregunta.

ameba dice Reinstate Monica

Creo que puede haber cierta confusión en esa sugerencia: parece combinar una probabilidad (que es lo que tiene una masa de punto o una distribución NB tiene en cero) con una densidad de probabilidad (que es lo que el valor de sería). Una densidad distinta de cero no le permite tratar con ceros exactos, ¡porque aún predice cero posibilidades de que surja cualquier valor de !

f (0, θ)

$f(0,\theta)$

0

$0$

whuber

Respuestas:

Esa es una pregunta interesante. Mi grupo de investigación ha estado utilizando la distribución a la que hace referencia durante algunos años en nuestro software de bioinformática disponible al público. Hasta donde yo sé, la distribución no tiene un nombre y no hay literatura al respecto. Si bien el artículo de Chandra et al (2012) citado por Aksakal está estrechamente relacionado, la distribución que consideran parece estar restringida a valores enteros para y no parecen dar una expresión explícita para el pdf. $r$

Para darle algunos antecedentes, la distribución NB se usa mucho en la investigación genómica para modelar datos de expresión génica que surgen de RNA-seq y tecnologías relacionadas. Los datos del recuento surgen a medida que el número de lecturas de secuencias de ADN o ARN extraídas de una muestra biológica que se puede mapear a cada gen. Por lo general, hay decenas de millones de lecturas de cada muestra biológica que se asignan a unos 25,000 genes. Alternativamente, uno podría tener muestras de ADN de las cuales las lecturas se asignan a ventanas genómicas. Nosotros y otros hemos popularizado un enfoque mediante el cual NB glms se ajustan a las lecturas de secuencia para cada gen, y se utilizan métodos empíricos de Bayes para moderar los estimadores de dispersión genewise (dispersion $\phi=1/r$ ) Este enfoque ha sido citado en decenas de miles de artículos de revistas en la literatura genómica, por lo que puede tener una idea de cuánto se usa.

Mi grupo mantiene el paquete de software edgeR R. Hace algunos años revisamos todo el paquete para que funcione con recuentos fraccionarios, utilizando una versión continua del NB pmf. Simplemente convertimos todos los coeficientes binomiales en el NB pmf a relaciones de funciones gamma y lo usamos como un pdf continuo (mixto). La motivación para esto fue que los recuentos de lectura de secuencia a veces pueden ser fraccionados debido a (1) mapeo ambiguo de lecturas al transcriptoma o genoma y / o (2) normalización de recuentos para corregir los efectos técnicos. Entonces, los recuentos a veces son recuentos esperados o recuentos estimados en lugar de conteos observados. Y, por supuesto, los recuentos de lectura pueden ser exactamente cero con probabilidad positiva. Nuestro enfoque garantiza que los resultados de inferencia de nuestro software sean continuos en los recuentos, coincidiendo exactamente con resultados discretos de NB cuando los recuentos estimados son enteros.

Hasta donde sé, no hay una forma cerrada para la constante de normalización en el pdf, ni hay formas cerradas para la media o la varianza. Cuando se considera que no existe una forma cerrada para la integral (la constante de Fransen-Robinson), está claro que no puede existir para la integral de la continua NB pdf tampoco. Sin embargo, me parece que las fórmulas tradicionales de media y varianza para el NB deberían seguir siendo buenas aproximaciones para el NB continuo. Además, la constante de normalización debe variar lentamente con los parámetros y, por lo tanto, puede ignorarse por tener una influencia insignificante en los cálculos de máxima verosimilitud.

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

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

Uno puede confirmar estas hipótesis por integración numérica. La distribución NB surge en bioinformática como una mezcla gamma de distribuciones de Poisson (ver el artículo binomial negativo de Wikipedia o McCarthy et al más abajo). La distribución continua de NB surge simplemente reemplazando la distribución de Poisson con su análogo continuo con pdf para donde es una constante de normalización para garantizar que la densidad se integre a 1. Supongamos, por ejemplo, que . La distribución de Poisson tiene pmf igual al pdf anterior en los enteros no negativos y, con

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$ , la media y la varianza de Poisson son iguales a 10. La integración numérica muestra que y la media y la varianza de la distribución continua son iguales a 10 a aproximadamente 4 cifras significativas. Por lo tanto, la constante de normalización es prácticamente 1 y la media y la varianza son casi exactamente las mismas que para la distribución discreta de Poisson. La aproximación se mejora aún más si agregamos una corrección de continuidad, integrando de a lugar de 0. Con la corrección de continuidad, todo es correcto (la constante de normalización es 1 y los momentos coinciden con Poisson discreto) a aproximadamente 6 figuras

a (10) = 1 / 0.999875

$a(10)=1/0.999875$

- 1 / 2

$-1/2$

\infty

$\infty$

En nuestro paquete edgeR, no necesitamos hacer ningún ajuste por el hecho de que hay masa en cero, porque siempre trabajamos con probabilidades de registro condicionales o con diferencias de probabilidad de registro y cualquier función delta se cancela de los cálculos. Esto es típico por cierto para glms con distribuciones de probabilidad mixtas. Alternativamente, podríamos considerar que la distribución no tiene masa en cero sino que tiene un soporte que comienza en -1/2 en lugar de en cero. Cualquiera de las perspectivas teóricas conduce a los mismos cálculos en la práctica.

Aunque hacemos un uso activo de la distribución NB continua, no hemos publicado nada explícitamente. Los artículos citados a continuación explican el enfoque NB de los datos genómicos, pero no discuten explícitamente la distribución continua de NB.

En resumen, no me sorprende que el artículo que está estudiando haya obtenido resultados razonables de una versión continua del pdf de NB, porque esa es también nuestra experiencia. El requisito clave es que deberíamos modelar las medias y las variaciones correctamente y eso estará bien siempre que los datos, sean enteros o no, exhiban la misma forma de relación cuadrática de media-varianza que la distribución NB.

Referencias

Robinson, M. y Smyth, GK (2008). Estimación de muestra pequeña de dispersión binomial negativa, con aplicaciones para datos SAGE . Bioestadística 9, 321-332.

Robinson, MD, y Smyth, GK (2007). Pruebas estadísticas moderadas para evaluar las diferencias en la abundancia de etiquetas . Bioinformática 23, 2881-2887.

McCarthy, DJ, Chen, Y, Smyth, GK (2012). Análisis de expresión diferencial de experimentos multifactoriales de RNA-Seq con respecto a la variación biológica . Nucleic Acids Research 40, 4288-4297.

Chen, Y, Lun, ATL y Smyth, GK (2014). Análisis de expresión diferencial de experimentos complejos de RNA-seq usando edgeR. En: Análisis estadístico de datos de secuencia de próxima generación, Somnath Datta y Daniel S Nettleton (eds), Springer, Nueva York, páginas 51-74. Preimpresión

Lun, ATL, Chen, Y y Smyth, GK (2016). Es DE-licious: una receta para análisis de expresión diferencial de experimentos de RNA-seq utilizando métodos de cuasi-verosimilitud en edgeR. Methods in Molecular Biology 1418, 391-416. Preimpresión

Chen Y, Lun ATL y Smyth, GK (2016). De las lecturas a los genes a las vías: análisis de expresión diferencial de experimentos de RNA-Seq usando Rsubread y la tubería de cuasi-verosimilitud edgeR . F1000 Investigación 5, 1438.

Gordon Smyth
fuente

Esto es extremadamente útil, @Gordon; Muchas gracias por tomarse el tiempo para escribirlo. También estoy trabajando con datos de RNA-seq, por lo que una respuesta desde esta perspectiva es especialmente valiosa (ahora he agregado la etiqueta [bioinformática] a la pregunta). Su trabajo trata sobre la expresión diferencial, mientras que mi trabajo actual trata sobre la agrupación (el documento que estaba leyendo es Harris et al. Sobre interneuronas CA1; biorxiv ). De todos modos, déjame hacerte un par de pequeñas preguntas / aclaraciones. [cont.]

ameba dice Reinstate Monica

- \log (p)

$-\log(p)$

r

$r$

r

$r$

r

$r$

r

$r$

@amoeba Gracias por la referencia biorxiv. (1) La derivación de NB como una mezcla de Poissons es bastante conocida, y está en nuestros documentos, por ejemplo, McCarthy et al. La derivación del NB continuo sigue simplemente sustituyendo Poisson continuo por Poisson. ¿Debo agregar esto a mi respuesta? Lo haría largo. No veo cómo el NB continuo podría ser útilmente representado como una mezcla de gammas. (2) No, la inflación cero es una complicación adicional diferente. Evitamos esa complicación en nuestro trabajo.

Gordon Smyth

@amoeba (3) Estimamos todos los parámetros. Es crucial estimar las dispersiones genéticas para lograr el control de la tasa de error, y esto debe hacerse con especial cuidado porque los tamaños de muestra son a menudo pequeños y la dimensión de los datos es enorme. Utilizamos un procedimiento complejo que implica la probabilidad de perfil ajustada (piense en REML) dentro de cada gen vinculado con un procedimiento empírico de Bayes de probabilidad ponderada entre genes. Los genes NB glms son luego equipados por ML con las dispersiones fijadas. Finalmente, los coeficientes se prueban usando pruebas F de cuasi-verosimilitud.

Gordon Smyth

Mira este artículo: Chandra, Nimai Kumar y Dilip Roy. Una versión continua de la distribución binomial negativa. Statistica 72, no. 1 (2012): 81 .

Se define en el documento como la función de supervivencia, que es un enfoque natural desde que se introdujo neg binomial en el análisis de confiabilidad:

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}$

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
fuente

¡Gracias! Echaré un vistazo a este artículo. (No fui yo quien votó negativamente.)

ameba dice Reinstate Monica el

@amoeba, no me preocupa el voto negativo, es internet :)

Aksakal

(Es extraño que esta respuesta haya sido rechazada ...) +1

whuber

x

$x$

@amoeba, el periódico tiene momentos, desafortunadamente no son los mismos que en NB

Aksakal