Création d'un modèle de régression logistique à trois niveaux dans pymc3

Question

Je tente de créer un modèle de régression logistique à trois niveaux dans pymc3. Il y a un niveau supérieur, un niveau intermédiaire et un niveau individuel, où les coefficients de niveau intermédiaire sont estimés à partir des coefficients de niveau supérieur. Cependant, j'ai de la difficulté à spécifier la structure de données appropriée pour le niveau intermédiaire.

Voici mon code:

with pm.Model() as model: 
    # Hyperpriors 
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.) 
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)  

    # Priors 
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top) 
    mid_level = [pm.Normal('mid_level_{}'.format(j), 
          mu=top_level[mid_to_top_idx[j]], 
          tau=mid_level_tau) 
       for j in range(k_mid)] 

    intercept = pm.Normal('intercept', mu=0., sd=100.) 

    # Model prediction 
    yhat = pm.invlogit(mid_level[mid_to_bot_idx] + intercept) 

    # Likelihood 
    yact = pm.Bernoulli('yact', p=yhat, observed=y) 

Je reçois l'erreur "only integer arrays with one element can be converted to an index" (ligne 16), qui je pense est lié au fait que la variable mid_level est une liste, pas un conteneur approprié de pymc. (Je ne vois pas non plus la classe Container dans le code source pymc3.)

Toute aide serait appréciée.

Modifier: Ajout des données fictives

y = np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0]) 
mid_to_bot_idx = np.array([0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 3, 2, 3, 3, 3, 3, 2, 2, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 3, 2, 3, 2, 3, 3, 2, 2, 3, 2, 2, 3, 3, 3, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 2]) 
mid_to_top_idx = np.array([0, 0, 1, 1]) 
k_top = 2 
k_mid = 4 

Modifier n ° 2:

Il semble y avoir plusieurs façons de résoudre ce problème, même si aucun sont tout à fait satisfaisant:

1) On peut recadrer le modèle comme:

with pm.Model() as model: 
    # Hyperpriors 
    top_level_tau = pm.HalfNormal('top_level_tau', sd=100.) 
    mid_level_tau = pm.HalfNormal('mid_level_tau', sd=100.)  

    # Priors 
    top_level = pm.Normal('top_level', mu=0., tau=top_level_tau, shape=k_top) 
    mid_level = pm.Normal('mid_level', mu=0., tau=mid_level_tau, shape=k_top) 
    intercept = pm.Normal('intercept', mu=0., sd=100.) 

    # Model prediction 
    yhat = pm.invlogit(top_level[top_to_bot_idx] + mid_level[mid_to_bot_idx] + intercept) 

    # Likelihood 
    yact = pm.Bernoulli('yact', p=yhat, observed=y) 

Cela semble fonctionner, bien que Je ne peux pas comprendre comment l'étendre au cas où la variance de niveau intermédiaire n'est pas constante pour tous les groupes de niveau intermédiaire.

2) On peut envelopper les coefficients de niveau intermédiaire dans un tenseur Théano utilisant theano.tensor.stack: à savoir,

import theano.tensor as tt 
mid_level = tt.stack([pm.Normal('mid_level_{}'.format(j), 
          mu=top_level[mid_to_top_idx[j]], 
          tau=mid_level_tau) 
       for j in range(k_mid)]) 

Mais cela semble fonctionner très lentement sur mon jeu de données réelles (30k observations) et cela rend le tracé peu pratique (chacun des coefficients mid_level obtient sa propre trace en utilisant pm.traceplot).

De toute façon, quelques conseils/contributions des développeurs seraient appréciés.

Answer 1

Transforme la solution est simple: il semble que toute distribution (comme pm.Normal) peut accepter un vecteur de moyens comme un argument, donc remplacer cette ligne

mid_level = [pm.Normal('mid_level_{}'.format(j), 
         mu=top_level[mid_to_top_idx[j]], 
         tau=mid_level_tau) 
      for j in range(k_mid)] 

avec cette

mid_level = pm.Normal('mid_level', 
         mu=top_level[mid_to_top_idx], 
         tau=mid_level_tau, 
         shape=k_mid) 

œuvres. La même méthode peut également être utilisée pour spécifier des écarts types individuels pour chacun des groupes de niveau intermédiaire.

EDIT: Correction d'une faute

Answer 2

Peu de changements (note j'ai changé yhat à thêta):

theta = pm.Deterministic('theta', pm.invlogit(sum(mid_level[i] for i in mid_to_bot_idx)+intercept)) 
yact = pm.Bernoulli('yact', p=theta, observed=y)