Stalling NUTS

[fonnesbeck]

The following model, which is a relatively simple network meta-analysis, works pretty well in PyMC 2:

# Difference means
mu = pm.Normal('mu', 0, 0.001, value=[0]*3)
tau = pm.Wishart('tau', 4, np.eye(3), value=np.eye(3)*0.001)
#sigma = pm.Uniform('sigma', 0, 100)
#tau = sigma**-2

# Study random effects (one for each outcome)
m = [pm.MvNormal('m_{}'.format(i), mu, tau, value=[0]*3) for i in range(len(unique_studies))]
mi = [m[i][j] for i,j in zip(study_id, outcome_id)]

# Covariates for treatment effects:
# high_intensity_eibi, low_intensity_eibi, low_intensity_parent_training, defined_school_program, 
# high_intensity_community, adjunctive_community_service, eclectic_school
beta = pm.Normal('beta', 0, 0.001, value=[0]*7)

# Index appropiate study outcome, add intervention indicators
theta = pm.Lambda('theta', lambda mi=mi, beta=beta: 
                  mi + np.dot(intervention_ind, beta))

# Likelihood
d = pm.Normal('d', theta, change_se**-2, observed=True, value=change_mean)

The idea is that there are 7 different treatments and 3 different outcome variables, and I am modeling them jointly. However, when I port this model to PyMC 3, the NUTS sampler has problems with it. Specifically, it runs for about 100 iterations and then stalls. I have no success getting a decent sample with Slice and Metropolis as well. So, I'm wondering if I've simply made a mistake in implementation:

with pm.Model() as model:

    # Difference means
    mu = pm.Normal('mu', 0, 0.001, shape=3)
    sigma = pm.Uniform('sigma', 0, 100)
    tau = sigma**-2

    # Study random effects (one for each outcome)
    m = pm.MvNormal('m', mu, TT.eye(3)*tau, shape=(len(unique_studies), 3))

    # Covariates for treatment effects
    beta = pm.Normal('beta', 0, 0.001, shape=7)

    # Index appropiate study outcome, add intervention indicators
    theta = m[study_id, outcome_id] + TT.dot(intervention_ind, beta)

    # Likelihood
    d = pm.Normal('d', theta, sd=TT.constant(change_se), observed=change_mean) 

It seems like a pretty straightforward one, so I'm hoping someone spots something I'm doing wrong. (Note I am using a simple diagonal matrix for the variance-covariance here, as I cannot get the Wishart to work).

[jsalvatier]

What if you replace MvNormal with Normal? does that work? That's the part
where I have the least confidence.

On Sun, May 4, 2014 at 6:53 PM, Chris Fonnesbeck
[email protected]:

The following model, which is a relatively simple network meta-analysis,
works pretty well in PyMC 2:

Difference means

mu = pm.Normal('mu', 0, 0.001, value=[0]_3)
tau = pm.Wishart('tau', 4, np.eye(3), value=np.eye(3)_0.001)
#sigma = pm.Uniform('sigma', 0, 100)
#tau = sigma**-2

Study random effects (one for each outcome)

m = [pm.MvNormal('m_{}'.format(i), mu, tau, value=[0]*3) for i in range(len(unique_studies))]
mi = [m[i][j] for i,j in zip(study_id, outcome_id)]

Covariates for treatment effects:

high_intensity_eibi, low_intensity_eibi, low_intensity_parent_training, defined_school_program,

high_intensity_community, adjunctive_community_service, eclectic_school

beta = pm.Normal('beta', 0, 0.001, value=[0]*7)

Index appropiate study outcome, add intervention indicators

theta = pm.Lambda('theta', lambda mi=mi, beta=beta:
mi + np.dot(intervention_ind, beta))

Likelihood

d = pm.Normal('d', theta, change_se**-2, observed=True, value=change_mean)

The idea is that there are 7 different treatments and 3 different outcome
variables, and I am modeling them jointly. However, when I port this model
to PyMC 3, the NUTS sampler has problems with it. Specifically, it runs for
about 100 iterations and then stalls. I have no success getting a decent
sample with Slice and Metropolis as well. So, I'm wondering if I've simply
made a mistake in implementation:

with pm.Model() as model:

# Difference means
mu = pm.Normal('mu', 0, 0.001, shape=3)
sigma = pm.Uniform('sigma', 0, 100)
tau = sigma**-2

# Study random effects (one for each outcome)
m = pm.MvNormal('m', mu, TT.eye(3)*tau, shape=(len(unique_studies), 3))

# Covariates for treatment effects
beta = pm.Normal('beta', 0, 0.001, shape=7)

# Index appropiate study outcome, add intervention indicators
theta = m[study_id, outcome_id] + TT.dot(intervention_ind, beta)

# Likelihood
d = pm.Normal('d', theta, sd=TT.constant(change_se), observed=change_mean)

It seems like a pretty straightforward one, so I'm hoping someone spots
something I'm doing wrong. (Note I am using a simple diagonal matrix for
the variance-covariance here, as I cannot get the Wishart to work).

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/540
.

[fonnesbeck]

Thanks. I did notice some issues with fitting models that use the Wishart, but did not know there were MVN concerns. I'm sure an independent normal model would work, but wanted to explicitly include correlation. I suppose I could use the conditional forms.

[jsalvatier]

That at least narrows down the problem to the MVN, we should take a look
and see what's going on.

On Mon, May 5, 2014 at 5:32 PM, Chris Fonnesbeck
[email protected]:

Thanks. I did notice some issues with fitting models that use the Wishart,
but did not know there were MVN concerns. I'm sure an independent normal
model would work, but wanted to explicitly include correlation. I suppose I
could use the conditional forms.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/540#issuecomment-42256799
.

[twiecki]

I'm sure this is related to the Wishart. Closing here and moving discussion over to #538.