Hi all,
This is my first time to use brms for Bayesian hierarchical model, please excuse me if my questions are stupid.
I have a data frame with the outcome of participants (the test result of if a participant gets infected or not, say n participants in total), age group (8 age groups in total), and regions (9 regions in total). 0 and 1 indicate the test result is negative or positive. Now I want to fit a Bayesian hierarchical model on the probability of being positive in each age group in each region. So my model has three levels, 9 regions, each region has 8 age groups, and each age group has a certain number of participants. Note that the number of participants in each age group in each region is not exactly the same, but very similar.
The mathematical expression of my model would be:
y_{ijk} \sim Bernoulli (1, p_{ij})
logit(p_{ij}) = \alpha + \alpha_{age[i,j]} + \alpha_{region[j]}
\alpha_{age_{[i,j]}} \sim Normal(0, \sigma_{age})
\alpha_{region_{[i]}} \sim Normal(0, \sigma_{region})
\alpha \sim Normal(0,1)
\sigma_{age} \sim HalfCauchy(0,1)
\sigma_{region} \sim HalfCauchy(0,1)
y_{ijk}: the test result of the kth individual in the ith age group in the jth region, where i = 1, …, 8, j = 1 and k = 1, … , n
\alpha_{age[i,j]}: ith age group in jth region
\alpha_{region[i]}: jth region
In brms, ignore priors for now, I write the fomula as:
brm (outcome ~ 1 + (1  region / age)).
Based on what I descbribed above, I have two following questions:

Did I build correct model mathematically? If so, is there any terms can be added in the model to improve the predictive perfomance?

Does the brm formula I used correctly reflect the mathematical model? If not, how should I change the formula?
Any comments and suggestions would be much appreciated!