I noticed when trying your model, that the logit_p calculation might have some type instability, because when I do @code_warntype, it appears red.
I encountered this issue multiple times, and always found it really annoying, not having a good way to deal with this multiple array indexing type thing that is quite frequent whenever you have multiple categorical indexes co-occuring.
But then I encountered CartesianIndex, it seems to good to be true though, but benchmarking it (outside the model at least) seemed significantly better in every respect, min/median/max/mean runtime, memory estimate and allocs estimate.
It really does seem to good to be true, but it also does not appear as type unstable in the model.
The way I did it then is this:
logit_p = g[treatment] .+ α[CartesianIndex.(treatment, actor)] .+ β[CartesianIndex.(treatment, block)]
I just feel like there must be a catch, can anyone confirm?