I’m evaluating a moderately large network of clinical programs, which is composed of about a hundred programs organized under eight hubs.

I have a race variable that can take one of six values: White, Black, Hispanic, AsianAPI, Multiracial, or Other.

I’m trying to estimate whether different racial groups have different treatment completion rates. I’d like to estimate overall variability by race, but I also want to see whether that racial variability is itself variable by program and by hub. The hope is to assess whether some hubs/programs are seeing different outcomes than others for certain racial groups.

Intuitively, it seems like this could be modeled as follows:

stan_glmer(completed ~
(1 | race) + # varying completion rates by race
(1 | program) + # varying completion rates by program
(1 | hub) + # varying completion rates by hub
(1 | program:race) + # varying completion rates by race by program?
(1 | hub:race), # varying completion rates by race by hub?
data = data,
family=binomial)

However, this seems somehow ridiculous? Is this a “reasonable” model specification?

I just can’t seem to wrap my head around the appropriate modeling strategy for this, and am hoping for help. Sorry if this is a dumb question.

It seems reasonable to me. You’re specifying a big interaction model with shrinkage. You may not care about the main effects so much as some balanced prediction conditional on the parameter draws (for example, treatment completion by race, assuming the same distributions of programs and hub, being sure to average across the actual completion rates, not the logits of completion).

You might take a look at Si et AL 2020 for an example of more structured priors for interactions. But the idea is the same as what you already have: allow the completion rates to vary by cell and shrink them towards the mean. It just differs in how the shrinkage works.

Thank you for the response, and that’s encouraging to see a similarly specified model in that paper.

I see what you mean about using this for prediction rather than worrying about fixed effects. In fact it would be nice to get reasonable estimates of the simple effects of race, too, rather than pure rates. I think I can get my race effect estimates by looking at the intercepts for (1|race). For example I guess I can subtract the logits for the black intercept from the white intercepts, exponentiate and take the mean for an odds ratio between them.

I agree that looks reasonable. One wrinkle in these situations is that at times it makes no sense to consider one effect without the other in that the main effect always shows up with one of the interactions also present. So sometimes you get a better fit with the interaction effect only. Worth checking, but what you have should work fine in most cases. Worst case, the main effect collects very little variance and the interaction gets the vast majority of it.

It’s not obvious to me the interaction should capture more variance than the main effect, but I’ve tried the model both ways and it does seem as though the effect of hub is more muted after the interaction is added. Is there an intuition for why that would be the case: that interactions somehow absorb more of the variance than the main effect itself? That seems backwards to me.

I’m not exactly sure; it’s been one of those “learn by doing” things for me. Perhaps when some interactions are more dominant in their relative effect? There could also be some latent effect currently without any variable representation that is not equally represented among the interactions.