Let's say I run a single-level logit model of "being depressed", as a function of age (in years) and opiate use (binary). The model looks like this
$Pr(Y_i=1)=\frac{exp(\beta_0+\beta_1Age_i+\beta_2OpiateUse_i)}{1+exp(\beta_0+\beta_1Age_i+\beta_2OpiateUse_i)}$
The model gives me values for the various coefficients (and the constant term), which I can then plug back into the formula to estimate the probability that a particular type of person will be depressed. So by plugging in age=27 and opiate use=0 I can get "the predicted probability of depression for a 27 year old who does not use opiates."
Now let's say that my respondents (indexed by i) are nested within school districts (indexed by j). To account for this I include a random effect $\zeta_j$ for the intercept, which I assume is distributed normally with a mean of zero and some estimated variance.
Now the model looks like this
$Pr(Y_{ij}=1)=\frac{exp(\beta_0+\beta_1Age_{ij}+\beta_2OpiateUse_{ij}+\zeta_j)}{1+exp(\beta_0+\beta_1Age_{ij}+\beta_2OpiateUse_{ij}+\zeta_j)}$
I want to generate the same kind of predicted probabilities as before, to know the predicted probability of depression for (say) a 27 year old who doesn't use opiates. This raises the question of what to do about the random effect $\zeta_j$. I'm aware that in earlier years when you asked a program (like Stata) to estimate predicted probabilities after multi-level logit models it could only do so by setting the random effects at zero. However, it seems like this approach is frowned upon and now the default approach (at least for Stata - when you run "margins, at(age=27 opiate=0)" after "melogit") is to "integrate over the random effects" or estimate probabilities "marginal on the random effect" (For example https://www.stata.com/stata14/marginal-margins/).
My problem is that I'm not entirely sure what this means or how (if at all) it affects the interpretation of the results.
So my first question is: what is the formula being used to evaluate a predicted probability from a multi-level logit model when "integrating over random effects"? I know that it's NOT simply replacing $\zeta_j$ with a particular value (e.g. the area under the normal curve describing the distribution of random effects), but what IS going on here?
My second question is: How do we interpret the results? If I substitute "27" for age and "0" for opiate use but "integrate over random effects" can the result still be interpreted as "the predicted probability of depression for a 27 year old who does not use opiates?" Or do I need to modify this statement in some way to reflect what is going on with the random effects?
