Multi-level logistic regression - probability received from the intercept is different from the original sample I am running multi-level logistic regression in order to perform case-mix adjustment (corrected estimates) for 100 clinics. I got some results but they are somehow suspicious to me. I noticed that received/modeled probability of the outcome variable (I calculated probability with the transformation from the logit: exp(x)/(1+exp(x)) is different from the orginal (sample) probability of a target variable. For example, sample probability of a target variable is 88.8%, meanwhile the received intercept of a multi-level logistic regression model is 2.20225 which is 90.0% after the transformation back to probability. This happens also on a unconditional model where covariates are not included.
I am wondering if this is something to expect? Because if I run ordinary logistic regression, both probabilities are almost identical. Also, when I run multi-level linear regression (in case of dependent variables of interval type) the sample mean is identical to the intercept of a multi-level model.
I centered all predictors (grand mean centering, weighted effect coding) but as I mentioned the situation is the same also in case of an unconditional model.
I am using R and a package glmer with a code below: 
# Unconditional Model (Intercept-only Model)
M0 <- glmer(Q15_binom ~ 1 + (1|id_hospital), data=df, family=binomial(link="logit"))

# Model with random intercept and slope and 5 predictors
M1 <- glmer(Q15_binom ~ 1 + health_c + chron_c + gender_c + Q43.wec + Q44.wec + 
            (1 + health_c + chron_c + gender_c + Q43.wec + Q44.wec|id_hospital), 
            data=df, family=binomial(link="logit"), 
            control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

I also tried to increase the nAGQ to use Gauss-Hermite approximation (for potential better accuracy) but the results are the same.
Many thanks for your feedback, help, ideas.
 A: So first let's see what happens when you run some logistic models on a toy data set. I'll modify the cbpp data set that downloads with lme4 to start. 
library(lme4)
cbpp$incidence_dich <- ifelse(cbpp$incidence==0, 0, 1)
cbpp <- cbpp[cbpp$herd!="2",]
cbpp <- cbpp[cbpp$herd!="8",] # For starters only using hers with 4 obs 

mean(cbpp$incidence_dich)

Okay so we see that the observed probability of incidence > 0 is around .6. 
What do we get out of our unconditional logistic regression? 
logit_1 <- glm(incidence_dich ~ 1, 
               data = cbpp, 
               family = binomial)

logit_to_prob <- function(x){
  prob <- exp(x)/(1 + exp(x))
  return(prob)
}

intercept <- coef(logit_1)[1]
logit_to_prob(intercept)

Alright, so far so good. As expected the uncoditional model just returns the observed probability of an incidence value being greater than 0. 
Now what happens when we deal with the nested structure of these data? Specifically, we want to specify a model that includes a random effect for herd. 
logit_2 <- glmer(incidence_dich ~ 1 + (1|herd), 
                 data = cbpp, 
                 family = binomial, 
                 control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

intercept <- fixef(logit_2)[1]
logit_to_prob(intercept)

Okay so we get a warning (boundary (singular) fit: see ?isSingular), so not great, but I don't really care about the model. For the purposes of your question, we'll ignore the warning (which is really just that the random effect variance was estimated to be 0). 
So what we see here is that we are replicating the observed probability in both modeling approaches. Now it would seem these data have a few things going on with them that guarantee that to be the case. First, the number of observations per herd was equivalent. Second, and much more importantly there was no random effect variance estimated by the model.
I'll deal with the latter in the next section. I'll simulate some data we can be confident contains a random intercept. This will also allow us to know what the "true" population value is so we can see if one approach is better than the other. 
set.seed(102619) #to ensure reproducibility
pop_prob <- .6
rand_eff <- function(x){ 
  odds <- x/(1-x)
  logit <- log(odds)
  # random effects should be normally distributed
  # but the assumption is that this is true on the link scale
  within_logit <- logit + rnorm(1, 0, 1)
  within_prob <- logit_to_prob(within_logit)
  return(within_prob)
}

# Observations for 25 groups
group <- 1:25

# 10 observations per group
n_obs <- 10

group_id <- c()
group_sim_prob <- c()
y_obs <- c()
for(i in 1:length(group)){
  group_prob <- rand_eff(pop_prob)
  for(o in 1:n_obs){
    group_id <- c(group_id, i)
    group_sim_prob <- c(group_sim_prob, group_prob)
    y <- sample(c(0,1), 
                size = 1,
                prob = c(1-group_prob, group_prob))
    y_obs <- c(y_obs, y)
  }
}

dat <- data.frame(group_id, group_sim_prob, y_obs)
head(dat)

Alright, so now we have a data set with a random intercept that varies by group. Let's see the observed probability that y_obs == 1 in these data. 
mean(dat$y_obs)

A value of .576 seems reasonable given a population value of .6. 
Now what do we get from a simple unconditional model?
logit_3 <- glm(y_obs ~ 1, 
               data = dat, 
               family = binomial)
intercept <- coef(logit_3)[1]
logit_to_prob(intercept)

Alright, matches up... so far so good. What about our fancier, generalized multilevel linear model - that accounts for nesting. 
logit_4 <- glmer(y_obs ~ 1 + (1|group_id), 
                 data = dat, 
                 family = binomial)

intercept <- fixef(logit_4)[1]
logit_to_prob(intercept)

Interesting... if you run the code you'll see that the logistic glm model replicates our observed mean perfectly (.576), but the probability estimated from the glmer model is estimated to be .5958997. Note that it may be of interest that the estimate for the nested model is much, much closer to the actual population value of .6
Maybe that is just one weird simulation right? Well what if we do the same thing 5000 times? 
n_sims <- 5000
run_name <- c()
sample_prob <- c()
pred_mlm <- c()

for(s in 1:n_sims){
  set.seed(NULL) # Because I used set.seed() earlier
  pop_prob <- .6

  # Observations for 25 groups
  group <- 1:25

  # 10 observations per group
  n_obs <- 10

  group_id <- c()
  group_sim_prob <- c()
  y_obs <- c()
  for(i in 1:length(group)){
    group_prob <- rand_eff(pop_prob)
    for(o in 1:n_obs){
      group_id <- c(group_id, i)
      group_sim_prob <- c(group_sim_prob, group_prob)
      y <- sample(c(0,1), 
                  size = 1,
                  prob = c(1-group_prob, group_prob))
      y_obs <- c(y_obs, y)
    }
  }
  dat <- data.frame(group_id, group_sim_prob, y_obs)
  logit_tmp <- glmer(y_obs ~ 1 + (1|group_id), 
                     data = dat, 
                     family = binomial, 
                     control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

  intercept <- fixef(logit_tmp)[1]
  run_name <- c(run_name, paste("Run", s))
  sample_prob <- c(sample_prob, mean(dat$y_obs))
  pred_mlm <- c(pred_mlm, logit_to_prob(intercept))
}

Let's now take a look at the results when plotted. We can see a consistent difference between the relugar glm version and the slightly more complex glmer versions of the model.  
library(bayesplot)
library(tidyverse)
sim_results<- data.frame (`Regular Logtistic Bias` = sample_prob - .6, 
                         `MLM Logistic Bias` = pred_mlm - .6)

mcmc_areas(sim_results, prob = .95)+
  geom_vline(xintercept = 0, color = 'red', lty = 'dashed')+
  theme_bw()+
  labs(title = "Difference in Model Estimates of Population Probability and Ground Truth", 
       caption = "Shaded region represents 95% interval")


What is happening here? Well in a nutshell the regular glm model appears to be ever-so-slightly undershooting the predicted population probability when random effects occur in the data. In general, I would tend to have more faith in the multilevel model when random effects are present in the data, all other things being equal. Remember the goal is to estimate the population value after all, not the sample value. 
Now there could be a number of other factors that I did not walkthrough (missing data on your x's, differences in number of observations per grouping unit, etc.). I would need to know more about your use case, or you can take this code walkthrough and modify it to more closely reflect your sample/population. 
