Análisis bayesiano jerárquico sobre la diferencia de proporciones

¿Por qué jerárquico? : He intentado investigar este problema, y ​​por lo que entiendo, este es un problema "jerárquico", porque estás haciendo observaciones sobre observaciones de una población, en lugar de hacer observaciones directas de esa población. Referencia: http://www.econ.umn.edu/~bajari/iosp07/rossi1.pdf

¿Por qué bayesiano? : Además, lo etiqueté como bayesiano porque puede existir una solución asintótica / frecuenta para un "diseño experimental" en el que a cada "celda" se le asignan suficientes observaciones, pero a efectos prácticos, conjuntos de datos del mundo real / no experimentales (o en menos mías) están escasamente pobladas. Hay muchos datos agregados, pero las celdas individuales pueden estar en blanco o tener solo algunas observaciones.

El modelo en abstracto:

Sea U una población de unidades a cada una de las cuales podemos aplicar un tratamiento, , de o , y de cada una de las cuales observamos observaciones de 1 o 0 también conocido como éxitos y fracasos. Sea para la probabilidad de que una observación del objeto bajo tratamiento resulte exitosa. Tenga en cuenta que puede estar correlacionado con .${u_1, u_2, u_3 ... u_N}$$T$$A$$B$$p_{iT}$$i \in \{1...N\}$$i$$T$$p_{iA}$$p_{iB}$

Para que el análisis sea factible, (a) asumimos que las distribuciones y son cada una una instancia de una familia específica de distribuciones, como la distribución beta, (b) y seleccionamos algunas distribuciones previas para hiperparámetros.$p_A$$p_B$

Ejemplo del modelo

Tengo una bolsa realmente grande de Magic 8 Balls. Cada bola 8, cuando se agita, puede revelar "Sí" o "No". Además, puedo sacudirlos con la pelota al revés o al revés (suponiendo que nuestras Magic 8 Balls funcionen 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, inicialmente no cree que esté correlacionado con ).$p_{iA}$$p_{iB}$

Preguntas:

Alguien ha muestreado al azar un número, , de unidades de la población, y para cada unidad ha tomado y registrado un número arbitrario de observaciones bajo tratamiento y un número arbitrario de observaciones bajo tratamiento . (En la práctica, en nuestro conjunto de datos, la mayoría de las unidades tendrán observaciones solo bajo un tratamiento)$n$$A$$B$

Teniendo en cuenta estos datos, necesito responder las siguientes preguntas:

1. Si tomo una nueva unidad al azar de la población, ¿cómo puedo calcular (analítica o estocásticamente) la distribución posterior conjunta de y ? (Principalmente para que podamos determinar la diferencia esperada en proporciones, )$u_x$$p_{xA}$$p_{xB}$$\Delta=p_{xA}-p_{xB}$
2. Para una unidad específica , , con observaciones de éxitos y fallas, ¿cómo puedo calcular (analítica o estocásticamente) la distribución posterior conjunta para y , nuevamente para construir una distribución de la diferencia en proporciones$u_y$$y \in \{1,2,3...,n\}$$s_y$$f_y$$p_{yA}$$p_{yB}$$\Delta_y$$p_{yA}-p_{yB}$

Pregunta : aunque realmente esperamos que y estén muy correlacionados, no estamos modelando eso explícitamente. En el caso probable de una solución estocástica, creo que esto causaría que algunos muestreadores, incluido Gibbs, sean menos efectivos para explorar la distribución posterior. ¿Es este el caso, y si es así, deberíamos usar una muestra diferente, modelar de alguna manera la correlación como una variable separada y transformar las distribuciones para que no estén correlacionadas, o simplemente ejecutar la muestra más tiempo?$p_A$$p_B$$p$

Criterio de respuesta

Estoy buscando una respuesta que:

• Tiene código, usando preferiblemente Python / PyMC, o excluyendo eso, JAGS, que puedo ejecutar

• Es capaz de manejar una entrada de miles de unidades.

• Dadas suficientes unidades y muestras, es capaz de generar distribuciones para , y como respuesta a la Pregunta 1, que se puede demostrar que coinciden con las distribuciones de población subyacentes (para verificar con las hojas de Excel proporcionadas en los "conjuntos de datos de desafío" sección)$p_A$$p_B$$\Delta$

• Dadas suficientes unidades y muestras, puede generar las distribuciones correctas para , y ( las hojas de Excel proporcionadas en la sección "conjuntos de datos de desafío" para verificar) como respuesta a la pregunta 2, y proporciona algunos justificación de por qué estas distribuciones son correctas$p_A$$p_B$$\Delta$

• Si la respuesta es similar al último modelo JAGS que publiqué, explique por qué funciona con anteriores dpar(0.5,1)pero no con anteriores dgamma(2,20). Gracias a Martyn Plummer en el foro JAGS por detectar el error en mi modelo JAGS. Al tratar de establecer un previo de Gamma (Forma = 2, Escala = 20), estaba llamando a lo dgamma(2,20)que realmente establece un prior de Gamma (Forma = 2, InverseScale = 20) = Gamma (Forma = 2, Escala = 0.05).

Conjuntos de datos de desafío

He generado algunos conjuntos de datos de muestra en Excel, con algunos escenarios posibles diferentes, cambiando la rigidez de las distribuciones p, la correlación entre ellas y facilitando el cambio de otras entradas. https://docs.google.com/file/d/0B_rPBjs4Cp0zLVBybU1nVnd0ZFU/edit?usp=sharing (~ 8Mb)

Mis intentos / soluciones parciales hasta la fecha

1) Descargué e instalé Python 2.7 y PyMC 2.2. Inicialmente, obtuve un modelo incorrecto para ejecutar, pero cuando intenté reformular el modelo, la extensión se congela. Al agregar / eliminar código, he determinado que el código que desencadena la congelación es mc.Binomial (...), aunque esta función funcionó en el primer modelo, así que supongo que hay algo mal con la forma en que he especificado 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) Descargué e instalé JAGS 3.4. Después de obtener una corrección de mis antecedentes del foro JAGS, 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))

Controlar

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 cualquiera que prefiera desarmar el modelo :)

En la web, las pruebas A / B típicas consideran el impacto en la tasa de conversión de una sola página o unidad de contenido, con posibles variaciones. Las soluciones típicas incluyen una prueba de significación clásica contra las hipótesis nulas o dos proporciones iguales, o más recientemente soluciones analíticas bayesianas que aprovechan la distribución beta como un conjugado previo.

En lugar de este enfoque de una sola unidad de contenido, que incidentalmente requeriría una gran cantidad de visitantes a cada unidad de la que estoy interesado en probar, queremos comparar variaciones en un proceso que genera múltiples unidades de contenido (no es un escenario inusual realmente ...) Entonces, en general, las unidades / páginas producidas por el proceso A o B tienen muchas visitas / datos, pero cada unidad individual puede tener solo unas pocas observaciones.

Como el tiempo para la recompensa expiró y no recibí respuestas, publicaré la respuesta que pude encontrar, aunque mi experiencia limitada con la inferencia bayesiana sugiere que esto debe tomarse con una dosis saludable de escepticismo.

I) Configuración

1. Descargué e instalé JAGS 3.4 y R 3.0.1
2. Instalé los paquetes rjags y coda iniciando R y usando install.packages (pkgname)

II) Modelo y datos: se utilizaron los archivos de modelo y datos ya detallados en la pregunta. Para responder la pregunta n. ° 1, agregué una observación adicional a los datos con las cuatro variables como 0.

III) Responder preguntas

1. Ejecuté JAGS en el modelo / datos (abra la línea de comando, vaya al directorio con sus archivos y escriba> jags-terminal Command.cmd. Se ejecutó y generó algunos archivos
2. En R, utilicé los siguientes comandos:
   • biblioteca ("rjags") para cargar el paquete instalado (y su paquete coda requerido)
   • setwd () para llegar al directorio donde estaban los archivos de salida
   • resultados = read.coda ("STEMchain1.txt", "STEMindex.txt")
3. Para responder la primera pregunta:
   • Como gráfico PDF, "gráfico (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 la observación de "todos los 0". (1 a n es para la variable p_A, n + 1 a 2n es para p_B y 2n + 1 a 3n es para delta)
4. Para responder la segunda pregunta, igual que arriba, pero cambie 3 * N -> 2 * N + y

No estoy seguro de que esta sea la forma correcta de obtener la respuesta, o si un modelo más complejo produciría mejores resultados, particularmente en el caso de la correlación, pero espero que eventualmente alguien más experimentado intervenga ...