Modélisation bayésienne hiérarchique des taux d'incidence

8

Le livre de Kevin Murphy aborde un problème bayésien hiérarchique classique (initialement abordé dans Johnson and Albert, 1999, p24):

Supposons que nous essayons d'estimer le taux de cancer dans villes. Dans chaque ville, nous échantillonnons un certain nombre de personnes et mesurons le nombre de personnes atteintes de cancer , où est le véritable taux de cancer dans la ville.NNjeXjePoubelle(Nje,θje)θje

Nous aimerions estimer les tout en permettant aux villes pauvres en données d'emprunter des forces statistiques aux villes riches en données.θje

Pour ce faire, il modélise afin que toutes les villes partagent le même avant, de sorte que les modèles finaux se présentent comme suit:θjeBêta(une,b)

p(,θ,η|N)=p(η)je=1NPoubelle(Xje|Nje,θje)Bêta(θje|η)

où .η=(une,b)

La partie cruciale de ce modèle est bien sûr (je cite) "que nous déduisons des données, car si nous le simplement à une constante, le sera conditionnellement indépendant, et là il n'y aura aucun flux d'informations entre eux ".η=(une,b)θje


J'essaie de modéliser cela dans PyMC , mais pour autant que je comprends, j'ai besoin d'un a priori pour et (je crois que c'est ci-dessus). Quel serait un bon avant pour ce modèle?unebp(η)

Au cas où cela aiderait, le code, tel que je l'ai maintenant, est:

bins = dict()
ps   = dict()
for i in range(N_cities):
    ps[i]   = pm.Beta("p_{}".format(i), alpha=a, beta=b)
    bins[i] = pm.Binomial('bin_{}'.format(i), p=ps[i],n=N_trials[i],  value=N_yes[i], observed=True)

mcmc = pm.MCMC([bins, ps])

où je crois que j'ai besoin d'un préalable pour aet b. Comment en choisir un?

Amelio Vazquez-Reina
la source

Réponses:

9

Un problème similaire est discuté dans Gelman, Bayesian Data Analysis , (2e éd., P. 128; 3e édition p. 110). Gelman suggère un préalablep(une,b)(une+b)-5/2, ce qui limite efficacement la "taille d'échantillon antérieure" une+b, et donc l'hyperprior bêta ne devrait pas être très informatif en soi. (Comme la quantitéune+baugmente, la variance de la distribution bêta diminue; dans ce cas, une variance antérieure plus petite contraint le "poids" des données observées dans la partie postérieure.) De plus, cet a priori ne définit pas siune>b, ou l'inverse, donc des distributions appropriées de paires de (une,b) sont déduits de toutes les données ensemble, comme vous préférez dans ce problème.

Gelman suggère également de re-paramétrer le modèle en termes de logit de la moyenne de θet la "taille de l'échantillon" du prieur. Donc, au lieu de faire l'inférence sur(une,b) directement, le problème concerne l'inférence sur les quantités transformées logit(uneune+b) et Journal(une+b). Cela admet des valeurs antérieures transformées sur le plan réel, plutôt que des valeurs antérieures non transformées qui doivent être strictement positives. De plus, cela accomplit une densité postérieure qui est plus diffuse lorsqu'elle est tracée. Cela rend les graphiques d'accompagnement plus lisibles, ce que je trouve utile.

Sycorax dit de réintégrer Monica
la source
1
Merci @ user777. Malheureusement, je ne suis pas en mesure d'utiliser des a priori multivariés pour le moment, j'ai donc laissé une question de suivi ici
Amelio Vazquez-Reina