If I’m reading this correctly, you have a cohort of N objects under observation of T time steps that belong to one group for each of 4 categories, where there are (5, 2, 13, 4) groups in these categories.
You then pre-compute the (N x 520) indicator matrix, reshape it into a (Nx1) matrix, with external indices representing the group information, and replicate it to a (TxN) matrix. Thus all([x in vars12_indexes for x in vars123_indexes]) should be true. So your hierarchical model looks like
\beta_{1a\gamma} \sim N(\beta_{1a}, \sigma_\gamma^2)
with
[\beta_1, \dots ,\beta_5](t) \sim \mathrm{RW}(\Sigma)(t)
If I got this correct, then the model looks well-specified. In the absence of other covariates on which to condition e.g. \beta_2(t), you can’t do much better than a random walk (though you might look into letting \beta_2(t) be a function of global principal components, which would reduce the number of parameters from 10+5T to 5 n_{PC}).
At any rate, this is a very complicated model with well over 1000 parameters, well into the territory where practical sampling requires a specifically-tailored MCMC. General samplers will certainly improve in the future, so I suspect this will not always be the case, but for now NUTS can’t work magic.
I can see two approaches for tractability:
- Use ADVI or low-rank ADVI (Householder flow)
- Marginalize out the pieces you don’t want; using N(0, t\Sigma) in place of the random walk; or choosing conjugate priors for categories you’re not actively analyzing and marginalize them out.
Good luck!