# Comparing posteriors distributions between models

I have two models (which I'm estimating by MCMC with Stan). There are more parameters in reality, but a simplified example is:

A: y ~ (1|group)

B: y ~ X + (1|group)

I then calculate the ICC in each model from the parameter chains.

I'd like to be able to say if the ICC is meaningfully lower in model B vs model A.... that is, that does X explain some of the variance attributable to group?

Is it reasonable to:

1. Make 2 sets of draws for the ICC (one from each model).
2. Sort each list and rank them
3. Compute the difference for between draws with the same rank, and summarise this distribution (e.g. as if it were a parameter posterior)

This seems like it summarises the information from the models which is relevant to my question, but it feels hacky, and I'm worried that comparing posterior distributions in this way doesn't actually have the interpretation I'd like. Any thoughts on this much appreciated.

Comparison of different posterior distributions of parameters is tricky because the parameter spaces are of different dimensions and the interpretations change depending on what you condition on.

Let me suggest an alternate route that to answering the question "does X explain some of the variance attributable to group?" that is useful in more circumstances because you are comparing different predictive distributions from the same model:

1. Get the posterior distribution of Model B
2. Draw from the posterior predictive distribution using Model B and the observed data on X. Calculate the distribution of the variance of these predictions.
3. Draw from the posterior predictive distribution using Model B but setting X to an appropriate constant (such as the same mean) for all observations. Calculate the distribution of the variance of these predictions.
4. Compare the two distributions of variance of predictions.

In the rstanarm R package, it would look like this:

library(rstanarm)
options(mc.cores = parallel::detectCores())
model_B <- stan_lmer(y ~ X + (1 | group), data = dataset)
PPD_1 <- posterior_predict(model_B)
var_1 <- apply(PPD_1, MARGIN = 1, FUN = var)
nd <- dataset
nd$X <- mean(nd$X)
PPD_2 <- posterior_predict(model_B, newdata = nd)
var_2 <- apply(PPD_2, MARGIN = 1, FUN = var)
summary(var_1 - var_2) # for example, better to use graphs

• Ben - thank you for this - it's really smart and does exactly what I wanted so I've accepted the answer. If you did know of anywhere this approach might be mentioned in a paper or textbook that would be really cool too, but I think the logic is so clear that I can explain to reviewers without a reference. – bjw Sep 5 '17 at 20:28
• It is consistent with the Gelman and Hill book, but that was pre-Stan so it might not be explained exactly like this. – Ben Goodrich Sep 5 '17 at 20:52
• Thanks. If anyone else is interested I've also just found the Gelman et al Bayesian Data Analysis book, Chapter 6, has some more on this: stat.columbia.edu/~gelman/bayescomputation/bdachapter6.pdf See p165 – bjw Sep 6 '17 at 8:09