# Bayesian chi-squared tests

I have a dataset with two groups of participants. Each participant performed a repeated measures task on which three types of errors could be made. I want to measure the difference in distributions of counts error types between groups.

As far as I understand, the conventional frequentist NHST approach here would use a $$\chi^2$$-test on the contingency table of error type counts aggregated by group.

There are a handful of Bayesian "$$\chi^2$$" tests I've found, but they typically operate in the same kind of NHST mindset (eg, Bayes factors for models that do/not encode independence assumptions). I want to report a full posterior over the differences instead.

I've developed a procedure that I think answers the question, but I'm not 100% sure of its appropriateness. I use a hierarchical model over per-subject count distributions. Loosely,

$$\mu_{population, i} \sim \mathcal{N}(0, 100)$$

$$\sigma_{population, i} \sim \mathcal{HalfCauchy}(5)$$

$$\mu_{group, i} \sim \mathcal{N}(\mu_{population, i}, \sigma_{population, i})$$

$$\sigma_{group, i} \sim \mathcal{HalfCauchy}(5)$$

$$\mu_{subject, i} \sim \mathcal{N}(\mu_{group, i}, \sigma_{group, i})$$

$$\theta_{subject, i} = \frac {e^{\mu_{subject, i}}} {\sum_{\hat{i}} e^{\mu_{subject, \hat{i}}}}$$

$$k_{subject, i} = \mathcal{Multi}(\vec{\theta}_{subject}, n_{subject})$$

Some reasoning about some of the choices in this model:

• subjects come from the same original population, so there should be a prior at the population level
• the Dirichlet distribution doesn't have well-behaved priors, so instead of using Dirichlet, I'm using a softmax function to produce a valid multinomial parameter from independent normally distributed variables, as in logistic regression. (I saw this suggested on Gelman's blog and a couple other places; seemed a reasonable approximation.)
• interestingly, when the per-group parameters are normally distributed, the sum of their rescaled squared differences is $$\chi^2$$ distributed, which makes a very natural contrast to sample the posterior over, as it mirrors the procedure and assumptions involved in the frequentist $$\chi^2$$ test:

$$\sum_j { \frac {(\mu_{0, j} - \mu_{1, j})^2} {\sigma_{0, j}^2 + \sigma_{1, j}^2} }$$

I can estimate the posterior distribution of this contrast scalar with HMC sampling (I use PyMC3).

So my questions are (1) whether this general approach is sane, (2) whether there are established similar approaches that are better, and (3) are there obvious/subtle modifications to the details of the model I describe?