Analyse bayésienne hiérarchique sur la différence de proportions

Pourquoi hiérarchique? : J'ai essayé de rechercher ce problème, et d'après ce que je comprends, c'est un problème "hiérarchique", parce que vous faites des observations sur les observations d'une population, plutôt que de faire des observations directes de cette population. Référence: http://www.econ.umn.edu/~bajari/iosp07/rossi1.pdf

Pourquoi bayésien? : De plus, je l'ai étiqueté comme bayésien car une solution asymptotique / fréquentiste peut exister pour un "plan expérimental" dans lequel chaque "cellule" se voit attribuer suffisamment d'observations, mais à des fins pratiques, des ensembles de données du monde réel / non expérimentaux (ou à moins la mienne) sont peu peuplées. Il y a beaucoup de données agrégées, mais les cellules individuelles peuvent être vides ou avoir seulement quelques observations.

Le modèle en résumé:

Soit U une population d'unités ${u_1, u_2, u_3 ... u_N}$ à chacun desquels on peut appliquer un traitement, $T$ , de n'importe quel $A$ ou $B$ , et à partir de chacun desquels nous observons iid observations de 1 ou 0 aka succès et échecs. Laisser $p_{iT}$ pour $i \in \{1...N\}$ être la probabilité qu'une observation à partir d'un objet $i$ sous traitement $T$ se traduit par un succès. Notez que $p_{iA}$ peut être corrélé avec $p_{iB}$ .

Pour rendre l'analyse possible, nous (a) supposons que les distributions $p_A$ et $p_B$ sont chacun une instance d'une famille spécifique de distributions, telles que la distribution bêta, (b) et sélectionnent certaines distributions antérieures d'hyperparamètres.

Exemple de modèle

J'ai un très gros sac de Magic 8 Balls. Chaque 8 balles, lorsqu'elle est secouée, peut révéler "Oui" ou "Non". De plus, je peux les secouer avec la balle à l'endroit ou à l'envers (en supposant que nos Magic 8 Balls fonctionnent à l'envers ...). L'orientation de la balle peut complètement changer la probabilité d'aboutir à un "Oui" ou à un "Non" (en d'autres termes, au départ, vous ne pensez pas que $p_{iA}$ est corrélé avec $p_{iB}$ ).

Des questions:

Quelqu'un a échantillonné un nombre au hasard, $n$ , des unités de la population, et pour chaque unité a pris et enregistré un nombre arbitraire d'observations en cours de traitement $A$ et un nombre arbitraire d'observations sous traitement $B$ . (Dans la pratique, dans notre ensemble de données, la plupart des unités n'auront des observations que sous un seul traitement)

Compte tenu de ces données, je dois répondre aux questions suivantes:

Si je prends une nouvelle unité $u_x$ au hasard à partir de la population, comment puis-je calculer (analytiquement ou stochastiquement) la distribution postérieure conjointe de $p_{xA}$ et $p_{xB}$ ? (Principalement pour pouvoir déterminer la différence de proportions attendue, $\Delta=p_{xA}-p_{xB}$ )
Pour une unité spécifique $u_y$ , $y \in \{1,2,3...,n\}$ , avec observations de $s_y$ succès et $f_y$ échecs, comment puis-je calculer (analytiquement ou stochastiquement) la distribution postérieure commune pour $p_{yA}$ et $p_{yB}$ , encore une fois pour construire une distribution $\Delta_y$ de la différence de proportions $p_{yA}-p_{yB}$

Question bonus: bien que nous nous attendions vraiment $p_A$ et $p_B$ pour être très corrélés, nous ne modélisons pas explicitement cela. Dans le cas probable d'une solution stochastique, je pense que cela rendrait certains échantillonneurs, dont Gibbs, moins efficaces pour explorer la distribution postérieure. Est-ce le cas et, dans l'affirmative, devrions-nous utiliser un échantillonneur différent, modéliser en quelque sorte la corrélation en tant que variable distincte et transformer $p$ distributions pour les rendre non corrélées, ou simplement exécuter l'échantillonneur plus longtemps?

Critères de réponse

Je cherche une réponse qui:

A du code, en utilisant de préférence Python / PyMC, ou sauf JAGS, que je suis capable d'exécuter
Est capable de gérer une entrée de milliers d'unités
Étant donné suffisamment d'unités et d'échantillons, est capable de produire des distributions pour $p_A$ , $p_B$ , et $\Delta$ comme réponse à la question 1, qui peut être montrée pour correspondre aux distributions de population sous-jacentes (à vérifier par rapport aux feuilles Excel fournies dans la section "ensembles de données de défi")
Avec suffisamment d'unités et d'échantillons, est capable de produire les bonnes distributions pour $p_A$ , $p_B$ , et $\Delta$ (Je vais utiliser les feuilles Excel fournies dans la section "jeux de données de défi" pour vérifier) comme réponse à la question 2, et explique pourquoi ces distributions sont correctes
~~Si la réponse est similaire au dernier modèle JAGS que j'ai publié, expliquez pourquoi cela fonctionne avec les dpar(0.5,1)priors mais pas avec les dgamma(2,20)priors.~~ Merci à Martyn Plummer sur le forum JAGS pour avoir détecté l'erreur dans mon modèle JAGS. En essayant de définir un a priori de Gamma (forme = 2, échelle = 20) dgamma(2,20), j'appelais qui définissait en fait un a priori de Gamma (forme = 2, InverseScale = 20) = Gamma (forme = 2, échelle = 0,05).

Jeux de données de défi

J'ai généré quelques exemples de jeux de données dans Excel, avec différents scénarios possibles, modifiant l'étanchéité des distributions p, la corrélation entre elles et facilitant la modification d'autres entrées. https://docs.google.com/file/d/0B_rPBjs4Cp0zLVBybU1nVnd0ZFU/edit?usp=sharing (~ 8 Mo)

Une visualisation des 4 jeux de données inclus dans le fichier Excel

Ma ou mes solutions partielles à ce jour

1) J'ai téléchargé et installé Python 2.7 et PyMC 2.2. Initialement, j'ai eu un modèle incorrect à exécuter, mais lorsque j'ai essayé de reformuler le modèle, l'extension se bloque. En ajoutant / supprimant du code, j'ai déterminé que le code qui déclenche le gel est mc.Binomial (...), bien que cette fonction ait fonctionné dans le premier modèle, je suppose donc que la façon dont j'ai spécifié le modèle.

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) J'ai téléchargé et installé JAGS 3.4. Après avoir obtenu une correction sur mes priors du forum JAGS, j'ai maintenant ce modèle, qui fonctionne avec succès:

Modèle

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

Les données

"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))

Contrôle

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)

L'application réelle pour tous ceux qui préfèrent choisir le modèle :)

Sur le Web, les tests A / B typiques prennent en compte l'impact sur le taux de conversion d'une seule page ou unité de contenu, avec des variations possibles. Les solutions typiques incluent un test de signification classique contre les hypothèses nulles ou deux proportions égales, ou plus récemment des solutions bayésiennes analytiques exploitant la distribution bêta comme un conjugué préalable.

Au lieu de cette approche à une seule unité de contenu, qui, incidemment, nécessiterait beaucoup de visiteurs pour chaque unité de je suis intéressé par les tests, nous voulons comparer les variations d'un processus qui génère plusieurs unités de contenu (ce n'est pas vraiment un scénario inhabituel ...). Donc, dans l'ensemble, les unités / pages produites par le processus A ou B ont beaucoup de visites / données, mais chaque unité individuelle peut n'avoir que quelques observations.

bayesian mcmc jags hierarchical-bayesian pymc Fabio Beltramini
la source

Vous avez ignoré des aspects cruciaux dans votre description. Après avoir dessiné uniformément un

x

$x$ , que se passe-t-il ensuite? Quelque chose doit, sinon vous n'avez absolument aucune information sur aucun des

p_{x}

$p_x$ ! Alors, qu'observez-vous? Observez-vous

x

$x$ et, disons, un tirage indépendant de

o_{x}

$o_x$ ? Ou observez-vous simplement le tirage de

o_{x}

$o_x$ et ne pas enregistrer

x

$x$ ? Pour estimer toutes ces distributions de densité, vous devez sûrement répéter le processus plusieurs fois - combien? Ou, par "estimation", supposez-vous que vous connaissez

p_{x}

$p_x$ et vous souhaitez calculer la distribution du résultat de votre expérience?

whuber

Désolé, peut-être l'expression "

x

$x$ est uniformément tiré de 1 ... N "était trompeur. Étant donné un grand nombre d'observations existantes sur divers

o_{x}

$o_x$ , Je veux déterminer la distribution de

p_{x}

$p_x$ . C'est quelque peu analogue à l'échantillonnage de 50 personnes d'une population, en observant leurs hauteurs, puis en demandant la distribution de la hauteur pour une personne «tirée au hasard de la population». Je viens d'utiliser l'uniforme puisque la population dans ce cas a été dénombrée, et je vois que c'est déroutant. (Poursuivant l'analogie, dans ce scénario, nous ne pouvons pas mesurer les gens, observons simplement quelques réponses binaires)

Fabio Beltramini

J'ai ajouté un exemple à la fin de la question qui, je l'espère, indique clairement ce que j'essaie de faire.

Fabio Beltramini

Merci. La question reste cependant confuse. Dans l'exemple, vous observez une paire ordonnée de valeurs: les temps roulés et les succès. Alors, quelle est exactement la "distribution" que vous cherchez? À quoi exactement correspondrait le gros sac de dés?

whuber

La paire de valeurs (rouleaux et succès) transmet des informations sur une variable latente, p. La distribution que je veux est la distribution des p dans le sac. Bien sûr, si nous ne supposons pas une classe de distributions à laquelle il appartient, il existe arbitrairement de nombreuses possibilités, mais si nous disons que p est par exemple bêta-distribué, il s'agit simplement de sélectionner les meilleurs paramètres d'ajustement pour cela. classe de distributions. Si cela aide, considérez maintenant au lieu de dés avec p = 1 / 4,1 / 6,1 / 8, que chaque objet a un p de Beta (2,2) ... ou Beta (4,4), etc.

Fabio Beltramini

Analyse bayésienne hiérarchique sur la différence de proportions

Réponses: