Test de diferencia de proporciones bayesiano
import pymc3
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Diferencia de proporciones con PyMC3¶
Introducción¶
En (Pham-Gia et al 2017) se describe con detalle una solución cerrada para estimar la distribucion de la diferencia de proporciones binomiales usando estadística bayesiana.
Aquí hago un intento por realizar el análisis de comparar dos proporciones utilizando el paquete para estadística bayesiana PyMC3. Aunque la solución cerrada (con series infinitas pero que sí convergen) es bastante complicada, sería óptimo implementarla y obtener cálculos exactos, sin embargo, es un ejercicio bastante didáctico resolver este problema aproximadamente utilizando simulaciones de monte carlo con PyMC3.
Método¶
El modelo es bastante sencillo:
model = pm.Model()
with model:
th1 = pymc3.Beta('th1', alpha=0.5, beta=0.5)
p1 = pymc3.Binomial('p1', n1, p = th1, observed=k1)
th2 = pymc3.Beta('th2', alpha=0.5, beta=0.5)
p2 = pymc3.Binomial('p2', n2, p = th2, observed=k2)
d = pymc3.Deterministic("d", th1 - th2)
model
Se tienen dos distribuciones previas Beta, dos distribuciones binomiales que van a modelar las proporciones y el término de la diferencia de las proporciones:
$d = \theta_1 – \theta_2$
He hecho dos modelos alternativos. El primero tiene una distribución previa Beta(0.5, 0.5) conocida como Jeffrey’s Prior, que genera una creencia previa de que la proporción que se está modelando es 0 o 1. La segunda alternativa es con Beta(1,1) que genera una creencia previa de que la proporción está distribuida uniformemente entre $[0,1]$.
def plotTrace(trace, varname):
hist,bins = np.histogram(trace[varname], bins=np.linspace(np.min(trace[varname]),np.max(trace[varname]),20), normed=True)
width = 0.7 * (bins[1] - bins[0])
center = (bins[:-1] + bins[1:]) / 2
plt.bar(center, hist, align='center', width=width)
def propsDiff(n1, k1, n2, k2):
model = pm.Model()
with model:
th1 = pymc3.Beta('th1', alpha=0.5, beta=0.5)
p1 = pymc3.Binomial('p1', n1, p = th1, observed=k1)
th2 = pymc3.Beta('th2', alpha=0.5, beta=0.5)
p2 = pymc3.Binomial('p2', n2, p = th2, observed=k2)
d = pymc3.Deterministic("d", th1 - th2)
print(model)
map_estimate = pymc3.find_MAP(model=model)
print("MAP quick Estimates:", map_estimate)
print("Difference, according to data = ", k1/n1-k2/n2)
with model:
# draw 10000 posterior samples
p_trace = pymc3.sample(10000)
plt.rcParams["figure.figsize"] = (12,3)
plt.subplot(1,3,1)
plotTrace(p_trace, "th1")
plt.subplot(1,3,2)
plotTrace(p_trace, "th2")
plt.subplot(1,3,3)
plotTrace(p_trace, "d")
print("HPD: ", pymc3.stats.hpd(p_trace))
def propsDiffB(n1, k1, n2, k2):
model = pm.Model()
with model:
th1 = pymc3.Beta('th1', alpha=1, beta=1)
p1 = pymc3.Binomial('p1', n1, p = th1, observed=k1)
th2 = pymc3.Beta('th2', alpha=1, beta=1)
p2 = pymc3.Binomial('p2', n2, p = th2, observed=k2)
d = pymc3.Deterministic("d", th1 - th2)
print(model)
map_estimate = pymc3.find_MAP(model=model)
print("MAP quick Estimates:", map_estimate)
print("Difference, according to data = ", k1/n1-k2/n2)
with model:
# draw 10000 posterior samples
p_trace = pymc3.sample(10000)
plt.rcParams["figure.figsize"] = (12,3)
plt.subplot(1,3,1)
plotTrace(p_trace, "th1")
plt.subplot(1,3,2)
plotTrace(p_trace, "th2")
plt.subplot(1,3,3)
plotTrace(p_trace, "d")
print("HPD: ", pymc3.stats.hpd(p_trace))
Resultados¶
Se muestran a continuación dos pruebas para ambos modelos. Una población con pocos casos y una población con mayor cantidad de casos.
En la primera se comparan dos grupos, uno de 60 y otro de 45 casos, con proporciones 0.2 y 0.125 respectivamente. La diferencia es de 0.075.
En la segunda se comparan grupos con 600 y 450 casos, con las mismas proporciones que arriba y con la misma diferencia.
Como es de esperar, la diferencia estimada para el primer caso incluye el cero en la highest posterior density (HPD), por lo que no se puede rechazar la hipótesis nula de que las proporciones son iguales. Mientras que en la segunda prueba el cero sí queda fuera del HPD y por lo tanto sí se puede concluir que hay una diferencia creíble entre las proporciones.
En cada caso se muestran estimaciones rápidas, la diferencia teórica, se hace un sampling pequeño (de 10000 samples), y se muestra los intervalos HPD. Luego se muestran las distribuciones posteriores de $\theta_1$, $\theta_2$ y $d$.
Los cálculos con el prior de Jeffrey muestran valores ligeramente lejanos a los valores teóricos mientras que el prior plano muestra estimaciones bastante cercanas a los valores teóricos.
Prueba con tamaño de muestra pequeño y con prior Beta(0.5,0.5)¶
propsDiff(50,10, 40, 5)
Prueba con tamaño de muestra pequeño y con prior Beta(1,1)¶
propsDiff(500,100, 400, 50)
Prueba con tamaño de muestra grande y con prior Beta(0.5,0.5)¶
propsDiffB(50,10, 40, 5)
Prueba con tamaño de muestra grande y con prior Beta(1,1)¶
propsDiffB(500,100, 400, 50)
Referencias¶
Pham-Gia, T., Thin, N.V. and Doan, P.P. (2017) Inferences on the Difference of Two Proportions: A Bayesian Approach. Open Journal of Statistics, 7, 1-15. https://doi.org/10.4236/ojs.2017.71001