avoid bi-modal posterior distribution for Phi_approx(normal)?

[Lei Zhang]

Hi Stan,

I am modeling a discounting factor (hence lies in the range [0,1]) in a hierarchical fashion. I first used Matt Trick (non-centered reparameterization) for sampling from a normal distribution and then transform them into the [0 1] range with Phi_approx. Code as below.

real disc_mu_pr; 

real<lower=0> disc_sd; 

vector[nSubjects] disc_raw;

vector<lower=0,upper=1+1e-10>[nSubjects] disc;

for (s in 1:nSubjects) {

  disc[s] <- Phi_approx( disc_mu_pr + disc_sd * disc_raw[s] );

}

disc_mu_pr ~ normal(0,1);

disc_sd    ~ cauchy(0,5); 

disc_raw   ~ normal(0,1);

The fitting result shows a bi-modal posterior distribution of the disc parameter for nearly all the participants. I wonder it might because that the group level disc_sd is too large, such that the transformed disc parameter is pushed to both ends (0 and 1). I noticed that the R-hat values are actually not bad (all below 1.02, most are 1.00). I paste some of the summaries and I attach the density and pairs plot in case you may want to have a look.

My question is, is there a way to avoid such a bi-modal posterior distribution? Should I use a more restricted prior for disc_sd (eg. half-normal(0,2)), or should I model the disc parameter as beta distribution and discard the Matt Trick?

I appreciate any suggestions. Thanks!

Lei

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

          mean se_mean    sd 2.5%  25%  50%  75% 97.5% n_eff Rhat

disc_mu   0.55    0.01  0.18 0.21 0.43 0.56 0.68  0.87  1036 1.01

disc_sd   9.34    1.05 14.45 1.11 3.06 5.08 9.37 51.72   190 1.02

disc[1]   0.56    0.01  0.41 0.00 0.07 0.67 0.99  1.00  1814 1.00

disc[2]   0.65    0.01  0.39 0.00 0.27 0.84 1.00  1.00  1280 1.00

disc[3]   0.73    0.01  0.32 0.00 0.53 0.87 1.00  1.00  2721 1.00

disc[4]   0.54    0.01  0.41 0.00 0.08 0.61 0.98  1.00  4000 1.00

disc[5]   0.85    0.00  0.22 0.16 0.77 0.96 1.00  1.00  2088 1.00

disc[6]   0.57    0.01  0.41 0.00 0.09 0.69 0.99  1.00  2293 1.00

disc[7]   0.42    0.01  0.40 0.00 0.01 0.30 0.88  1.00  4000 1.00

disc[8]   0.41    0.01  0.42 0.00 0.00 0.23 0.89  1.00   814 1.00

disc[9]   0.64    0.01  0.38 0.00 0.30 0.80 0.99  1.00  4000 1.00

disc[10]  0.25    0.01  0.31 0.00 0.00 0.10 0.41  1.00  2022 1.00

[Krzysztof Sakrejda]

I can't tell from your question if these results are a problem: a) do they conflict with how you think this model should work? b) could you post the full code so it's easier to tell what the parameters are?  From a purely modeling perspective this could all be fine.

Krzysztof

[Lei Zhang]

Hi Krzysztof,

sorry for the late response. We have a cognitive task, and the model is a bit complex, with the 'disc' parameter being only part of the model. But I still attach it in case you may want to take a look. You could search the relevant lines of the disc parameter.

The concern I have is that, is such a bimodal posterior distribution a stable estimate of the parameter? If so, does this mean a parameter value of nearly 0 or 1 is more plausible than other parameter values? As we are also interested in the posterior mean (to correlate with fMRI data later), averaging across such a bimodal distribution doesn't seem good to us in this particular scenario. 

Hope to hear your (and of course others') advice.

Lei

[Bob Carpenter]

Point estimates are pretty much always lousy summaries of multimodal densities.

For example, if you take a beta(0.01, 0.01), there is a lot of mass around 0 and
1 and not much elsewhere. It's not multimodal per se, because there are no modes,
but clearly no point estimate gives you a good summary of the behavior.

Averaging across the modes using full Bayes is exactly what you want to do.
That is, integrate over the posterior.

The usual problem is that it's hard to sample from a bimodal posterior unless
you can get a single chain to move among the modes.

- Bob

> --
> 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.
> <RevLearn_RLbeta_alt4_c_w_v6_1lr.stan>

[Krzysztof Sakrejda]

On Tuesday, April 26, 2016 at 1:04:02 PM UTC-4, Lei Zhang wrote:

Hi Krzysztof,

sorry for the late response. We have a cognitive task, and the model is a bit complex, with the 'disc' parameter being only part of the model. But I still attach it in case you may want to take a look. You could search the relevant lines of the disc parameter.

That was helpful, thanks.  I don't understand the full model but I see that the disc parameters are transformed parameters so the parameters themselves are not biomodal and the sampler should have no trouble working correctly.

 

The concern I have is that, is such a bimodal posterior distribution a stable estimate of the parameter?

I'm not sure what stable means in this context but some models just have multi-modal posteriors so if the model

is right that's what you're stuck interpreting.   

 

If so, does this mean a parameter value of nearly 0 or 1 is more plausible than other parameter values?

That's the right interpretation.

 

As we are also interested in the posterior mean (to correlate with fMRI data later), averaging across such a bimodal distribution doesn't seem good to us in this particular scenario. 

Yes in this case a posterior mean will not be a useful summary if this is supposed to correlate with a measurement of this process later.  That said, the reason for the bi-modality is that disc_sd is large.  Have you thought about taking the disc_sd ~ cauchy(0,5) to something less extreme?  Given that this is a random effect plugged through Phi_approx I think a disc_sd ~ normal(0,1) would be a good start.  It might get rid of the bi-modality if the main issue is that disc_sd is just poorly informed from the data.  Of course you would then need to check that it's not driving your results later but I think it should be ok.

[Nathaniel Haines]

On Thursday, April 21, 2016 at 10:03:36 AM UTC-4, Lei Zhang wrote:

Hi all, 

I found this post after observing some strange behavior with this same transformation (i.e. Phi_approx transformation after using the Matt trick). 

In summary, I think this issue may be caused by the prior, which is placing a large mass on the extremes of the parameter bounds after the probit transformation is applied. I have attached code below (plot_priors.R) that can be used to draw samples from the prior distribution for a parameter being sampled in this way. The code can be modified by changing the distribution of the scale parameter from a cauchy to a normal (line 7), and you can also change the other values used in the sampling (e.g. the hyper_sd scale on line 11). Using the same parameterization as in the model you have attached, I get prior distributions of the form shown in image_1.png (attached). The image shows what the prior looks like both with (left) and without (right) the probit transformation being applied. On the left, you can see that the prior mass is heavily localized around the parameter bounds. On the right, the prior looks like what you would expect (Note that I am "zooming in" by limiting the x axis between -200 and 200, else you cannot see the form). I think this may explain why bimodal posteriors are common when using this type of parameterization? But I could be wrong–I wanted to post here to see what everyone else thinks. 

Interestingly, when the scale parameter on the cauchy distribution is set very low, the prior distribution over the probit transformed space becomes relatively uniform. image_2.png shows an example of the prior when the scale parameter in the cauchy is set to 0.001 as opposed to 5. Also, the right panel now appears more or less normally distributed. 

So, I guess my question is now this–what is the best way to approach this problem? Can the scale parameter for the cauchy be set low, or will this not allow for individual level differences in that parameter estimate (in the case of Lei's example code, differences among disc)? Another potential solution would be to truncate the parameter being estimated so that values under 0.01 and over 0.99 are not allowed (using the code example above; see image_3.png for a visual depiction):

 vector<lower=0.01,upper=0.99>[nSubjects] disc;

However, I understand that hard constraints like this are not ideal. I guess the problem then is coming up with a good soft constraint. 

If anyone has information on the best steps to take in this situation, it would be greatly appreciated. Please let me know if anything is unclear. 

Best, 

Nate

[Bob Carpenter]

Thanks---that's right. The Cauchy always has fat tails, but you
can keep most of the mass more central and avoid getting a U-shaped
result.

Truncating won't help as it'll just pile up against the truncation
points, and still be problematic. You need to control where the
parameters get estimated.

- Bob

> However, I am not sure if this is an appropriate way to do this.

>
> If anyone has information on the best steps to take in this situation, it would be greatly appreciated. Please let me know if anything is unclear.
>
> Best,
> Nate
>

> --
> 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.

> <plot_priors.R><image_1.png><image_2.png>

[Nathaniel Haines]

Hi Bob, 

Thanks a lot for the response! I will look into better soft constraints. I will post results back here if I find anything that works best in this situation. I posted this same discussion here (https://groups.google.com/forum/#!category-topic/stan-users/general/g3KEiNvHOdE) since I was unsure whether or not this older post would get attention, so I will post a follow-up there as well (if that is OK). Sorry about the double post. 

Best,

Nate