Difference between monitor means and get_posterior_mean?

[Stephen Martin]

Hey all,

I am trying to make an R package that uses Stan as the backend (first package I've ever made, so hopefully the code isn't miserable).

Anyway, I noticed that when I get means for a parameter using, e.g.,:
monitor(extract(bDifOut,pars='parameterName',permuted=FALSE)[,1,,drop=FALSE])[,1], I get a different estimate than when I use get_posterior_mean(bDifOut, pars='parameterName') and look at chain 1's mean.

What's up with that?

Example output:

> monitor(extract(bDifOut,pars='pi_sigma',permuted=FALSE)[,1,,drop=FALSE],print=TRUE,digits=3) #monitor(chain1)
Inference for the input samples (1 chains: each with iter=500; warmup=250):

          mean se_mean    sd  2.5%   25%   50%   75% 97.5% n_eff  Rhat
pi_sigma 0.667   0.065 0.487 0.043 0.271 0.597 0.917 1.905    57 0.999

For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
> get_posterior_mean(bDifOut,pars='pi_sigma')
         mean-chain:1 mean-chain:2 mean-chain:3 mean-chain:4 mean-all chains
pi_sigma    0.7357284    0.6942247    0.6349823    0.7632244         0.70704

[Ben Goodrich]

On Thursday, November 17, 2016 at 5:10:04 PM UTC-5, Stephen Martin wrote:

Anyway, I noticed that when I get means for a parameter using, e.g.,:
monitor(extract(bDifOut,pars='parameterName',permuted=FALSE)[,1,,drop=FALSE])[,1], I get a different estimate than when I use get_posterior_mean(bDifOut, pars='parameterName') and look at chain 1's mean.

What's up with that?

monitor() by default throws away have the draws, so you should call it with warmup = 0 if you are already passing an array that is strictly post-warmup.

Ben

[Stephen Martin]

Ah, thank you. That is... a frustrating default, since extract() defaults to no warmups extracted.

[Stephen Martin]

Really, I wish there were an easier way to get posterior summaries of one or multiple chains, by choice.
E.g., to get posterior summaries for chains 1:2, but not 3:4, you have to type:
monitor(extract(stanfit, pars='parameters',permuted=FALSE)[,1:2,],print=FALSE,warmup=0)
That is crazy verbose in my opinion. Ideally, you could just do: summary(stanfit, chains=1:2).

For it to be programmable, it looks even worse. To get a particular chain or multiple chains with a column of your choice:

monitor(extract(stanfit, pars=parameters, permuted=FALSE)[,chains,,drop=FALSE],print=FALSE,warmup=0)[,column]

I think I may just write a wrapper function around that monstrosity for the purpose of my package, unless there is better way I simply don't know about?

[Ben Goodrich]

On Friday, November 18, 2016 at 12:27:05 AM UTC-5, Stephen Martin wrote:

Really, I wish there were an easier way to get posterior summaries of one or multiple chains, by choice.

We usually don't want people summarizing a subset of chains. Either all the chains have converged and it makes sense to summarize all of them or some of the chains haven't converged and it doesn't make sense to summarize any of them.

[Stephen Martin]

Sorry, but it's necessary for some problems. In mixture models, there can be arbitrary between-chain label switching. You can tell it's arbitrary visibly in the trace plot. For cases like that, you *do* want to examine similar chains for your inference. It's not ideal, but it most certainly is the correct decision to make. Some parameters are invariant to the swapped labels, but many aren't, and for those that aren't invariant, examining similar chains will yield the correct inference.

This is theoretically the correct thing to do, too. The identical, but swapped modes do indeed exist (so of course, Stan is correct), but one's decision is not 'what is the expected value assuming that the labels will swap in my population', but 'taking one of the arbitrary labelling solutions, what should my inference be?'.

Other times, there can be authentic multimodality, or true multiple solutions, and the inferences *would* depend on which mode you picked, so to speak. Even in this case, you may want to pick apart the separate modes to better understand why so much density is in separate regions.
E.g., I have a model that predicts latent group membership from demographic data, these groups which differ in how certain items correlate to a latent ability estimate (literally, it's latent group DIF detection for 2PL models and others). Usually there are only label-swapped modes, and your inferences are essentially the same no matter the mode. However, sometimes the predictors can produce true multimodality, such that one predictor separates people into two groups that implicates set A of items to be DIF, and another mode finds that a second predictor separates people into two other groups that implicate set B to be DIF. It is useful to examine both of these modes as well as the full posterior, because it's not arbitrarily multimodal and the posterior is actually expressing real uncertainty about the generative model.

Making per-chain summaries harder to come by on principle doesn't discourage people from doing it, it just makes using Stan unnecessarily more difficult. Put up a big warning when people select chains about the potential issues with non-convergence and be done with it. Otherwise, trust that the user is picking apart chains for a justifiable or at *least* a diagnostic reason.

[Michael Betancourt]

This is a common misunderstanding.  When properly implementing full Bayes,

all of the chains would explore all of the modes in the right proportions and

they would indeed look the same.  ANY application of a Markov chain where

multiple chains behave differently, however, is formally not a “correct thing to 

do” theoretically as it will necessarily be only an approximation to full Bayes.

For example, if you can’t isolate a unique labeling in the model specification

then you can’t guarantee that you can identify different labelings with different

chains.  Eventually some chains will transition between modes and you’ll be

left with a convolution of multiple labels that can’t be separated out so easily,

if at all.

This is even more important when the multi modality is not from a labeling

degeneracy.  In that case the modes will have different weights, but having

different chains explore different modes will tell you _nothing_ about what

those weights should be so you have no way of combining the individual

chains in a well-posed manner.

I have been studying Markov chains intensely for the past four years and

I still don’t really trust myself to pull off an approximate multimodal analysis 

like this (the more you learn, the more you see what can go wrong).  I think 

it’s fair not to expect most users to be able to pull it off, either.  Many of our

design choices have been motivated to facilitate good practices and

frustrate bad practices.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Stephen Martin]

I agree with everything you say, of course, and you know far more about markov chains than I do.

Until there is a better method to handle multimodality in mixture models though, basing one's inference on the full posterior (even *if* all chains explored all modes) would be pointless.
On the other hand, if there *is* true multimodality that would change one's decisions/inference, having all chains explore all modes in correct proportion would help you understand the probability of solutions; the expected value would be difficult, but examining the posteriors could still tell you something about the probability of certain generative models.

Currently, we don't have that, and the best we can do is to examine the multiple modes separately. They still tell you something about the solution, but nothing about the relative probability of each. Not sure if lp__ would be useful there though (2 chains hit one mode, 2 hit another, but the latter two have a greater lp__, so perhaps that is the 'more probable' solution? I'm not sure how useful that metric would be).

The good thing is that the model I'm implementing is really more an exploratory procedure. It finds combinations of items that latent groups seem to operate differently on, so it detects biased items among latent groups. Sometimes it can find two modal solutions and one doesn't necessarily need to pick one or the other, but to just see 1) which items seem to be detected for that mode and 2) what are the characteristics of that mode's latent groups?
For that purpose, I think it's fair to pick apart the modes. One may decide which solution is more *important* (say, if 12 items are detected as biased when male and female seem to be the latent groups, and the error in the latent assignment seems low, it may be more important than if 4 items are detected as biased when age seems to divide individuals into latent groups, but the assignment is error-prone).

Thoughts on this? Yes, it's not using the full marginal posterior, but it still seems to lead one to make the 'correct' decisions in the context of my model (verifiably so, because I'm generating the data).

Also, are there any plans for improving mixture modeling in Bayes, such that this multimodality isn't so much an issue? I know this is a relatively unsolved problem, even in ML estimation (EM has the same multiple-solution problem).

[Bob Carpenter]

For the pure label switching problem in symmetric mixture models,
you want to concentrate on inferences that are not sensitive to label.
That is, inferences that are exchangeable under component index permutation.
So where you might want to say item i is in group k with some probability
and item j is in group k with some probability to try to tell which items
are in which groups, that's not well formulated because the group k isn't
identified in the model. But what you can do is look at the probability
that i and j are in the same group---that event is exchangeable under
component re-indexing.

Having said that, I would sometimes like to be able to explore different chains
easily just to visualize what's going wrong in some cases. So a wrapper
to make it easier to extract chains out of a fit object would be helpful for
that or a way to visualize or summarize with only some of the chains.

- Bob

[Stephen Martin]

For the first paragraph -- Indeed, that's what is invariant.
Currently, the likelihood is a mixture (with one component itself being a mixture):


So in this case, \Delta_j is invariant to label switching and will simply tell you whether the item's parameters do indeed seem to differ between some latent groups.
However, we can also predict \pi_{ik} (on a log odds scale via softmax regression) from other characteristics about the individual. Because the labels can switch, the coefficients for that regression can also arbitrarily switch signs or labels. This is usually obvious, so this is why I just look at one mode or another for the inference, e.g., "what seems to predict membership to group k=2, in a set of groups that seem to differ on items with high \Delta_j values?"

I wound up just writing a wrapper around summary() called summary_chain() that takes the same arguments as summary() but with an additional chains argument (e.g., summary_chains(..., chains = c(1,2)) would give you the posterior summary over chains 1 and 2).

[Andrew Gelman]

My current favored approach is to get rid of the aliasing using constraints on the parameters. This does become more difficult in multiple dimensions, and there's more work to be done here. I've been thinking about the general problem for many years, and it seems to me that some restriction is the way to go. Also, it's not necessary that the restriction be a hard constraint; a soft constraint could work too, maybe.
A

[Stephen Martin]

For my model, the best option is to set a positivity constraint (half-cauchy, half-normal) on the beta coefficient predicting latent group membership, but it's not generalizable (not everyone using my model would use latent regressors).

I've tried an ordering constraint on item coefficients (e.g., difficulty of item 1 is ordered across latent groups), but this seems to prevent within-chain switching moreso than between-chain modes. I'm guessing the issue is that the ordering constraint doesn't affect the overall likelihood enough to prevent two modes from occurring. One mode will just have a slightly higher lp__ than the other, since only one item has a constraint (out of, e.g., 20 items).

Ignoring that multimodality (which, again, is easy to detect either through trace plots, summary stats, or looking at the latent regressors' magnitude and signs), the model *can* have authentic multimodality, which I think I can fix by using a more informative prior on the latent regression 'error term' (categorization error). Rather than allow the error to be small or fairly large, I can put a gamma(1,1) or exp(1) prior on the error to favor models that use the latent regressors more for latent group assignment than allowing people to freely be assigned to each cluster. Testing this out now.

I hope to publish this model soon (along with an R package for fitting it), I'm just trying to figure out what the best 'default' priors might be on realistic datasets. The only 'problem' with this model is that it works great with MCMC, but optimizing fails miserably and vb does fairly poorly as well. But to be fair, it's a weird model --- If there are J items, then the model is fitting J + 1 mixture models, something that optimizing() should fail at, and it isn't remotely mvn, so vb fails too.

[Bob Carpenter]

> On Nov 20, 2016, at 10:58 PM, Stephen Martin <hwki...@gmail.com> wrote:
>
> For my model, the best option is to set a positivity constraint (half-cauchy, half-normal) on the beta coefficient predicting latent group membership, but it's not generalizable (not everyone using my model would use latent regressors).
>
> I've tried an ordering constraint on item coefficients (e.g., difficulty of item 1 is ordered across latent groups), but this seems to prevent within-chain switching moreso than between-chain modes. I'm guessing the issue is that the ordering constraint doesn't affect the overall likelihood enough to prevent two modes from occurring. One mode will just have a slightly higher lp__ than the other, since only one item has a constraint (out of, e.g., 20 items).

I think you're right that this can happen---we've seen behavior
like this in IRT models before.

The prior might help, but it's hard to tell when you've truly
gotten rid of multimodality other than the label-switching
exchangeability.

- Bob