I'm trying to correctly set up Bayesian parameter estimation for a mixed-design study with one 2-level between-groups independent variable and one 2-level within-subjects independent variable. The dependent variable is a binary outcome (ie Bernoulli distributed), and I'm interested in differences both between groups and between within-subjects variables.
For one variable it's straightforward to set up the problem using a hierarchical model. (see Doing Bayesian Data Analysis, BEST paper, 2010 TiCS review, etc):
y ~ Bernoulli(y | theta_subject) theta_subject ~ Beta(theta_subject | mu_group, sigma_group) mu_group ~ Beta(mu_group | mu_pop, sigma_pop) sigma_group ~ Gamma(sigma_group | scale_pop, rate_pop) etc - truncating priors for brevity
One can then estimate group differences by sampling differences between
mu_group for each between-subjects level.
I could use separate hierarchical models for both between-groups and within-subjects variables, but this would essentially ignore any interaction between the two variables. And because Bayesian ANOVA typically assumes a linear model, I can't apply it directly here.
What I can do is use an ANOVA-like hierarchical setup, but use a logistic model instead of a linear model. However, this would be a pretty substantive change in model configuration wrt the single-dependent-variable model. In the hierarchical single-variable setup, the group means are parameters for beta distribution (having reparameterized beta distribution in terms of mean and variance) drawn from a beta distribution. But in the logistic regression setup, the group means would be parameters for normal distributions for regression coefficients. This would require a different interpretation of the model.
So my question is a bit open-ended: what is the best way to set up a Bayesian hierarchical model for differences estimation in a 2x2 between/within experimental design?