Test HBMNL model on simulated data

[marcdotson]

Here's the principled workflow with simulated data using the generated quantities block: https://betanalpha.github.io/assets/case_studies/principled_bayesian_workflow.html.

[dgmiller]

HBMNL.stan is tested on simulated data and is able to recover parameters. This issue is fixed.

[marcdotson]

Compare parameterization and model output to what is being implemented with flipChoice: https://www.r-bloggers.com/how-to-use-hierarchical-bayes-for-choice-modeling-in-r/.

[dgmiller]

This issue has been resolved.

[dgmiller]

@marcdotson In response to your model comparison I was bugged enough to run some tests on the centered vs non-centered parameterizations. I varied the adapt_delta parameter. Behold the images. I ran them on a very small data set of 2 respondents, 15 choice tasks, 2 alternatives, 3 feature-levels.

Histograms of the betas (the row corresponds to the respondent, the column to the feature)

NON-CENTERED adapt_delta = .8, .9, .95

CENTERED adapt_delta = .8, .9, .95

The histograms look roughly similar which is not too interesting. Now let's look at the "trace plots". The x-axis is the iteration or draw (out of 4000, 1000 from each of 4 chains) and the y-axis is the value for the beta drawn.

NON-CENTERED adapt_delta = .8, .9, .95

CENTERED adapt_delta = .8, .9, .95

Now that's interesting. Both models seem to have systematic divergences (data on that below) but you can see the bias in the centered model whereas the uncentered model is still biased but less obviously so. Consider the table below.

model, adapt_delta, div rate, convergence (.9<Rhat<1.1)
n      .80          12.83%    Y
n      .90          05.30%    Y
n      .95          00.35%    Y
c      .80          00.75%    N
c      .90          01.03%    Y
c      .95          00.83%    N

I don't really have a conclusion for this yet except that we should maybe investigate this further, especially if we are using these in an ensemble. Could you try reproducing these results?

[marcdotson]

Thanks for reopening this, @dgmiller. I haven't gone back to simulations, but I thought I would share some interested results in comparing parameterizations from another study using real-live data!

The object of interest in this study are the estimates of the upper-level coefficient matrix. I've run it for centered, non-centered, and "conjugate" (a random-walk Metropolis Hastings MCMC using a conjugate inverse Wishart prior).

First the marginal posteriors for the upper-level:

NON-CENTERED

CENTERED

CONJUGATE

As in simulation, they're all the same, so nothing really interesting there.

Now the traceplots (note that the conjugate model just had one chain; no fancy parallelization in the legacy code borne from my PhD).

NON-CENTERED

CENTERED

CONJUGATE

Also, not very different, though the centered parameterization has traceplots that are wider, which corroborates with what you saw in simulation. That model also complained about tree depth and maybe one other HMC-specific diagnostic.

TBH, this is the very reason that the non-centered parameterization is encouraged, isn't it? In a hierarchical setting there are strange geometries in the posterior that need to be dealt with. I know I made a stink about my initial comparison, but I think there is probably something strange in terms of the simulated data or in the fact that there was just one simulated dataset. I'll get back to testing it since I think that will be a helpful something we could share on the blog-that-is-soon-to-be. But I thought these real-world comparisons were an interesting point to consider.

As a bonus, here's how long each parameterization took:

# Centered (2000 iterations): 17990.1 seconds (Total)
# Noncentered (2000 iterations): 5634.98 seconds (Total)
# Conjugate (100000 iterations): Total Time Elapsed (in Hours):  14.89