I'm struggling with a question that I haven't seen talked about anywhere. Potentially it just illustrates some ignorance I have about how effects are estimated in mixed models.
Let's set up the problem such that we have an educational intervention treatment that is being administered to students within school settings. The hope is that this improves the student passing rate of their grade, a binary variable pass. However, the students being embedded in schools and classrooms (school, class) needs to be controlled for with random effects. Assuming the treatment effect is constant for each student we could model this at the student level with with:
glmer(pass ~ treatment + (1 | school/class), family='binomial' ...)
However, let's assume that the treatment effect is not distributed as randomly as we would like. Some classrooms administered it correctly and half their class were treated while the other half were not treated (split_treatment). Some classrooms treated all their students by accident (all_treatment), others forgot completely that they were supposed to administer this treatment (missed_treatment). The number of classes fitting each criteria are illustrated below:
split_treatment: 50 all_treatment: 500 missed_treatment: 500
In this situation, is the fixed effect for treatment essentially only estimated from the split_treatment classes that contain students who were both treated and untreated?
How would we account for or investigate potential disagreement between models? For example, imagine the specific mixed-effects model above indicates that the treatment improves student's pass rate. However, a solely fixed effects model for
glm(pass ~ treatment) shows a negative impact of the treatment, ie. we essentially see Simpson's paradox with this data. Should we just chalk it up to the selection effects of treatment administration? I'm struggling with the fact that it feels like we ignore 1000 all_treatment or missed_treatment classes while potentially only estimating based on the 50 split_treatment classes.