Jerárquico Bayesiano de análisis en la diferencia de proporciones

Question

Por Qué Jerárquica? : He tratado de investigar este problema, y por lo que entiendo, esto es un "jerárquica" problema, porque están haciendo observaciones sobre las observaciones de una población, en lugar de hacer observaciones directas de la población. Referencia: http://www.econ.umn.edu/~bajari/iosp07/rossi1.pdf

Por Qué Bayesiana? : También, me ha etiquetado como Bayesiano porque un asintótica/frecuentista solución puede existir para un "diseño experimental", donde cada "celda" es disponer de suficientes observaciones, pero para los fines prácticos del mundo real/no de los conjuntos de datos experimentales (o al menos la mía) están escasamente pobladas. Hay un montón de datos agregados, pero las células individuales pueden ser en blanco o tiene sólo un par de observaciones.

El modelo en resumen:

Sea U una población de unidades de ${u_1, u_2, u_3 ... u_N}$ a cada uno de los cuales se puede aplicar un tratamiento, $T$, de $A$ o $B$, y de cada uno de los cuales observamos yo.yo.d. observaciones de ya sea 1 o 0.k.una. los éxitos y los fracasos. Deje $p_{iT}$ $i \in \{1...N\}$ la probabilidad de que una observación de objeto $i$ tratamiento $T$ resultados en un éxito. Tenga en cuenta que $p_{iA}$ puede estar correlacionada con $p_{iB}$.

Para hacer el análisis factible, nosotros (a) suponga que las distribuciones $p_A$ $p_B$ cada una instancia de una determinada familia de distribuciones, tales como la distribución beta, (b) y seleccione algunas de las distribuciones previas para hyperparameters.

Ejemplo de la modelo

Tengo una gran bolsa de la Magia 8 Bolas. Cada 8-ball, cuando se agitan, puede revelar "Sí" o "No". También, puedo sacudir con la bola de la derecha, o al revés (suponiendo que nuestra Magia 8 Pelotas de trabajar al revés...). La orientación de la pelota puede cambiar completamente la probabilidad de resultar en un "Sí" o un "No" (en otras palabras, en un principio que no tienen creencia de que $p_{iA}$ se correlaciona con $p_{iB}$).

Preguntas:

Alguien tiene una muestra aleatoria de un número, $n$, de las unidades de la población, y para cada unidad ha tomado y registrado un número arbitrario de las observaciones bajo tratamiento $A$ y un número arbitrario de las observaciones bajo tratamiento $B$. (En la práctica, en nuestro conjunto de datos, la mayoría de las unidades tienen observaciones sólo en virtud de un tratamiento)

Teniendo en cuenta estos datos, tengo que responder a las siguientes preguntas:

1. Si hago una nueva unidad de $u_x$ al azar de la población, ¿cómo puedo calcular (analíticamente o estocástica) de la articulación posterior distribución de $p_{xA}$$p_{xB}$? (Principalmente por lo que podemos determinar la diferencia esperada en proporciones, $\Delta=p_{xA}-p_{xB}$)
2. Para una unidad específica $u_y$, $y \in \{1,2,3...,n\}$, con las observaciones de $s_y$ éxitos y $f_y$ fallas, ¿cómo puedo calcular (analíticamente o estocástica) de la articulación posterior distribución por $p_{yA}$$p_{yB}$, de nuevo a construir una distribución de $\Delta_y$ de la diferencia en las proporciones $p_{yA}-p_{yB}$

Pregunta extra: a Pesar de que realmente esperan $p_A$ $p_B$ a ser muy correlacionadas, no estamos de modelar explícitamente que. En el caso probable de una estocástico solución, creo que esto podría causar que algunas muestras, incluyendo Gibbs, el ser menos eficaz en la exploración de la distribución posterior. Es este el caso, y si es así, debemos utilizar otro sampler, de alguna manera, el modelo de la correlación como una variable independiente y transformar la $p$ distribuciones para hacerlos no correlacionados, o simplemente ejecutar el sampler no?

Responder a los criterios de

Estoy en busca de una respuesta que:

• Tiene el código, utilizando preferiblemente Python/PyMC, o restricción que, ENTRECORTADO, que soy capaz de ejecutar

• Es capaz de manejar una entrada de miles de unidades

• Dada la suficiente cantidad de unidades y de las muestras, es capaz de salida de las distribuciones de $p_A$, $p_B$, y $\Delta$ como una respuesta a la Pregunta 1, que puede ser mostrado para que coincida con el subyacente de la población de las distribuciones (se comprobará en contra de hojas de excel que figuran en el "desafío de los conjuntos de datos" en la sección)

• Dada la suficiente cantidad de unidades y de las muestras, es capaz de emitir el derecho de las distribuciones de $p_A$, $p_B$, y $\Delta$ (voy a utilizar las hojas de excel que figuran en el "desafío de los conjuntos de datos" de la sección a de verificación) como una respuesta a la Pregunta 2, y ofrece algunas razones de por qué estas distribuciones son correctos

• Si la respuesta es similar a la de las últimas batallas contra la modelo que he publicado, explicar por qué funciona con dpar(0.5,1) priores, pero no con dgamma(2,20) de los priores. Gracias a Martyn Plummer en las EXIGENCIAS de foro para la captura de los errores en mi ENTRECORTADO modelo. En el intento de fijar una antes de Gamma(Shape=2,Escala=20), me estaba llamando dgamma(2,20) que, en realidad, un conjunto de antes de la Gamma(Shape=2,InverseScale=20) = Gamma(Shape=2,Escala=0.05).

Desafío de los conjuntos de datos

He generado una muestra de algunas bases de datos en Excel, con un par de diferentes escenarios posibles, el cambio de la opresión de la p distribuciones, la correlación entre ellos, y lo que es fácil cambiar a otros insumos. https://docs.google.com/file/d/0B_rPBjs4Cp0zLVBybU1nVnd0ZFU/edit?usp=sharing (~8 mb)

Mi intento/solución parcial(s) a la fecha

1) he descargado e instalado Python 2.7 y PyMC 2.2. Inicialmente, tengo una incorrecta modelo a correr, pero cuando traté de reformular el modelo, la extensión se congela. Mediante la adición/eliminación de código, he determinado el código que activa el congelamiento es el mc.Binomial(...) , aunque esta función no funcionó en el primer modelo, por lo que supongo que hay algo mal con cómo he especificado el modelo.

import pymc as mc
import numpy as np
import scipy.stats as stats
from __future__ import division
cases=[0,0]
for case in range(2):
    if case==0:
        # Taken from the sample datasets excel sheet, Focused Correlated p's
        c_A_arr, n_A_arr, c_B_arr, n_B_arr=np.loadtxt("data/TightCorr.tsv", unpack=True)

    if case==1:
        # Taken from the sample datasets excel sheet, Focused Uncorrelated p's
        c_A_arr, n_A_arr, c_B_arr, n_B_arr=np.loadtxt("data/TightUncorr.tsv", unpack=True)

    scale=20.0
    alpha_A=mc.Uniform("alpha_A", 1,scale)
    beta_A=mc.Uniform("beta_A", 1,scale)
    alpha_B=mc.Uniform("alpha_B", 1,scale)
    beta_B=mc.Uniform("beta_B", 1,scale)
    p_A=mc.Beta("p_A",alpha=alpha_A,beta=beta_A)
    p_B=mc.Beta("p_B",alpha=alpha_B,beta=beta_B)

    @mc.deterministic
    def delta(p_A=p_A,p_B=p_B):
        return p_A-p_B

    obs_n_A=mc.DiscreteUniform("obs_n_A",lower=0,upper=20,observed=True, value=n_A_arr)
    obs_n_B=mc.DiscreteUniform("obs_n_B",lower=0,upper=20,observed=True, value=n_B_arr)

    obs_c_A=mc.Binomial("obs_c_A",n=obs_n_A,p=p_A, observed=True, value=c_A_arr)
    obs_c_B=mc.Binomial("obs_c_B",n=obs_n_B,p=p_B, observed=True, value=c_B_arr)


    model = mc.Model([alpha_A,beta_A,alpha_B,beta_B,p_A,p_B,delta,obs_n_A,obs_n_B,obs_c_A,obs_c_B])
    cases[case] = mc.MCMC(model)
    cases[case].sample(24000, 12000, 2)

    lift_samples = cases[case].trace('delta')[:]

    ax = plt.subplot(211+case)
    figsize(12.5,5)
    plt.title("Posterior distributions of lift from 0 to T")
    plt.hist(lift_samples, histtype='stepfilled', bins=30, alpha=0.8,
             label="posterior of lift", color="#7A68A6", normed=True)
    plt.vlines(0, 0, 4, color="k", linestyles="--", lw=1)
    plt.xlim([-1, 1])

2) he descargado e instalado ENTRECORTADO 3.4. Después de conseguir una corrección en mi priores de las EXIGENCIAS de foro, ahora tengo este modelo, que se ejecuta con éxito:

Modelo

var alpha_A, beta_A, alpha_B, beta_B, p_A[N], p_B[N], delta[N], n_A[N], n_B[N], c_A[N], c_B[N];
model {
    for (i in 1:N) {
        c_A[i] ~ dbin(p_A[i],n_A[i])
        c_B[i] ~ dbin(p_B[i],n_B[i])
        p_A[i] ~ dbeta(alpha_A,beta_A)
        p_B[i] ~ dbeta(alpha_B,beta_B)
        delta[i] <- p_A[i]-p_B[i]
    }
    alpha_A ~ dgamma(1,0.05)
    alpha_B ~ dgamma(1,0.05)
    beta_A ~ dgamma(1,0.05)
    beta_B ~ dgamma(1,0.05)
}

Datos

"N" <- 60
"c_A" <- structure(c(0,6,0,3,0,8,0,4,0,6,1,5,0,5,0,7,0,3,0,7,0,4,0,5,0,4,0,5,0,4,0,2,0,4,0,5,0,8,2,7,0,6,0,3,0,3,0,8,0,4,0,4,2,6,0,7,0,3,0,1))
"c_B" <- structure(c(5,0,2,2,2,0,2,0,2,0,0,0,5,0,4,0,3,1,2,0,2,0,2,0,0,0,3,0,6,0,4,1,5,0,2,0,6,0,1,0,2,0,4,0,4,1,1,0,3,0,5,0,0,0,5,0,2,0,7,1))
"n_A" <- structure(c(0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3))
"n_B" <- structure(c(9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3,9,0,9,0,9,0,3,0,9,0,9,0,9,3))

Control

model in Try1.bug
data in Try1.r
compile, nchains(2)
initialize
update 400
monitor set p_A, thin(3)
monitor set p_B, thin(3)
monitor set delta, thin(3)
update 1000
coda *, stem(Try1)

La aplicación real para cualquier persona que prefiere escoger aparte de la modelo :)

En la web, típico de las pruebas a/B considera que el impacto sobre la tasa de conversión de una sola página o de la unidad de contenido, con variaciones posibles. Soluciones típicas incluyen un clásico de significación de la prueba en contra de la hipótesis nula o dos proporciones iguales, o, más recientemente, analítica Bayesiano soluciones de aprovechamiento de la distribución beta como un conjugado antes.

En lugar de esto de una sola unidad de contenido enfoque, que por cierto se requeriría una gran cantidad de visitantes a cada unidad de estoy interesado en las pruebas, queremos comparar las variaciones en un proceso que genera múltiples unidades de contenido (no un escenario inusual, en realidad...). Así que en conjunto, las unidades/páginas producidas por el proceso de a o B tienen una gran cantidad de visitas/datos, pero cada unidad individual puede tener sólo un par de observaciones.

Answer 1

Desde el momento de la recompensa caducado y no he recibido respuestas, voy a publicar la respuesta que yo era capaz de llegar, a pesar de mi limitada experiencia con inferencia Bayesiana sugiere que esto debe ser tomado con una saludable dosis de escepticismo.

I) El Programa De Instalación

1. He descargado y instalado EXIGENCIAS de 3.4 y R 3.0.1
2. He instalado los paquetes rjags y coda iniciando R y el uso de instalar.paquetes(pkgname)

II) Modelo de Datos e - Utiliza el modelo de datos y archivos ya se ha detallado en la pregunta. Para contestar la pregunta #1, he añadido una observación adicional sobre los datos con las cuatro variables como 0.

III) Responder a preguntas

1. Corrí ENTRECORTADO en el modelo de datos (abrir la línea de comandos, vaya al directorio con los archivos, y el tipo > entrecortado-terminal de Comandos.cmd . Se corrió y salida de un par de archivos
2. En R, he utilizado los siguientes comandos:
   • biblioteca("rjags") para cargar el paquete instalado (y sus paquete de coda)
   • setwd() para llegar al directorio donde el outputfiles fueron
   • resultados=leer.coda("STEMchain1.txt","STEMindex.txt")
3. Para responder a la primera pregunta:
   • Como un PDF de la trama, "plot(resultados[,3*N])"
   • Como cuantiles, "cuantil(resultados[,3*N],c(0.025,0.25,0.5,0.75,0.975))"
   • Donde N es el número de observaciones, y la última observación corresponde a la posición de "0 todos" de observación. (De 1 a n es variable p_A, n+1 2n es para p_B, y 2n+1 3n es para delta)
4. Para responder a la segunda pregunta, igual que el anterior, pero el cambio 3*N -> 2*N+y

No estoy seguro de que esta es la manera correcta para obtener la respuesta, o si de modelos más complejos podría producir mejores resultados, especialmente en el caso de la correlación, pero esperemos que finalmente alguien más experimentado va a meter la cuchara...