SBC for compositional model with multinomial observations

Modeling
, , ,
1

I’m using sbc() in rstan on a latent compositional model with multinomial observations, and am getting very U-shaped histograms for all my parameters. My real case is a factorial experiment with overdispersion and blocking, but I get the same behaviour even in the simplest case where my compositional model has only a mean vector:

\mu \sim N(\mathbf 0, 4 \mathbf I_4)
\mathbf x = \mu
\rho = \text{ilr}^{-1}(\mathbf x)
\text{counts} \sim \text{multinomial}(n, \rho)

Here, \mu \in \mathbb R^4 is the mean vector, and \text{ilr}^{-1} is an inverse isometric logratio transformation so that \rho \in \mathbb S^4 is a parameter for a multinomial with 5 categories. Each observation is a multinomial with around 100 trials, and I have 60 observations (which I would have thought is enough for models of this kind to work). My priors are based at least to some extent on genuine prior knowledge, so I wouldn’t want to change them all that much. There’s nothing else obviously wrong (e.g. n_eff, Rhat all look OK, I’m not seeing divergent transitions, and posterior distributions of \rho are roughly centred on sample estimates even in my more complicated real case).

Any idea where things are going wrong?

Example attached. I’m using R 3.6.3, rstan 2.21.1 with stan headers 2.21.0-5, 64-bit Ubuntu 18.04.

aureliaexample.R (616 Bytes) aureliaexample.stan (1.8 KB) sbc.pdf (5.6 KB)

2

Hmm, I haven’t used rstan sbc before, so this might not be so helpful, but it looks like you’re working with a simplified model, and I think just keep simplifying until the model does work.

Like, maybe do:

for(i in 1:npanels)
  counts[i] ~ normal(mu[i], sigma);

as the likelihood. I read the sbc docs and your model and it’s not obvious to me where things are breaking.

Maybe double check that n[i] == sum(counts[i]).

Also can:

z = exp(tVinv * x);
y = z / sum(z);

be replaced by y = softmax(tVinv * x)?

3

Turns out I was just misunderstanding how to use SBC. I had my observations (counts) in the data block (I had thought that they need to be there in the same way as when fitting the model to data, but will be ignored). Instead, they need to be in the transformed data block (and named exactly as in the model block, i.e. counts rather than y in my case).

There was no problem with the model once that was fixed.

You’re right that using softmax in the ilrinv() function simplifies the code.

aureliaexample.stan (1.8 KB)

1 Like