Question about multilevel logistic model (mixed intercept logistic model) in Stan code

Question

I'm trying to write stan code for multilevel logistic regression. The model that I tried is a mixed intercept logistic model with two predictors. The first level is children level and the second level is mom level. When I tried to match the summary result of the code I wrote versus the one generated by function stan_glmer(), the results of fixed intercept did not match. First, the data I used as below:

library(rstanarm)
library(rstan)

data(guImmun, package = "mlmRev")
summary(guImmun)
require(dplyr)
guImmun <- guImmun %>%
  mutate(immun = ifelse(immun == "N",0,1))

Second, the stan code was written as below:

data {
    int N; // number of obs 
    int M; // number of groups 
    int K; // number of predictors
    
    int y[N]; // outcome
    row_vector[K] x[N]; // predictors
    int g[N];    // map obs to groups (kids to women)
}
parameters {
    real alpha;
    real a[M]; 
    vector[K] beta;
    real<lower=0,upper=10> sigma;  
}
model {
  alpha ~ normal(0,1);
  a ~ normal(0,sigma);
  beta ~ normal(0,1);
  for(n in 1:N) {
    y[n] ~ bernoulli(inv_logit( alpha + a[g[n]] + x[n]*beta));
  }
}

Fitting data to the model:

guI_data <- list(g=as.integer(guImmun$mom),
                y=guImmun$immun,
                x=data.frame(guImmun$kid2p, guImmun$mom25p),
                N=nrow(guImmun),
                K=2,
                M=nlevels(guImmun$mom))
ranIntFit <- stan(file = "first_model.stan", data = guI_data,
                  iter = 500, chains = 1)
summary(ranIntFit, pars = c("alpha", "beta", "a[1]", "a[2]", "a[3]", "sigma"),
                          probs = c(0.025, 0.975),
                          digits = 2)

I got the following result: results of written model However, if I use stan_glmer() function, the result would be presented as follows.

M1_stanglmer <- stan_glmer(immun ~ kid2p + mom25p + (1 | mom), 
                           family = binomial("logit"), 
                           data = guImmun,
                           iter = 500,
                           chains = 1,
                           seed = 349)
print(M1_stanglmer, digits = 2)

But the results do not match, especially the result of fixed intercept. Results generated by the stan_glmer() function

Could anyone help me figure out what's wrong with my code? Thanks!

Answer 1

So, I wouldn't expect an exact equivalence between your model in Stan and the version implemented in stan_glmer, but for models that sample well it's reasonable to expect the estimates to be similar.

However, in your case, there is yet another issue impacting your estimates:

The covariates you are using in the guI_Data$x object have values in {1,2}, where the typical implementation would use values in {0,1} to represent a binary covariate. This is what is done in stan_glmer.

This coding is apparent if you use glimpse to inspect the data structure:

> library(tidyverse)
> glimpse(guI_data)
List of 6
 $ g: int [1:2159] 1 2 3 4 5 5 6 7 7 8 ...
 $ y: num [1:2159] 1 0 0 0 0 1 1 1 1 1 ...
 $ x:'data.frame':  2159 obs. of  2 variables:
  ..$ guImmun.kid2p : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 1 2 2 ...
  ..$ guImmun.mom25p: Factor w/ 2 levels "N","Y": 1 1 1 1 2 1 1 2 2 2 ...
 $ N: int 2159
 $ K: num 2
 $ M: int 1595

This is having the biggest impact in your intercept parameter, since the intercept represents the expected linear predictor when all covariates are 0. That value will often change when covariates are transformed or added.

Actually, I would expect that the estimated coefficients from your fit and the stan_glmer model are in fact similar, once you take this transformation into consideration.

For example, consider:

• define: x_m = x + 1
• Your model (m): yhat_m = alpha_m + x_m1*beta_m1 + x_m2*beta_m2
• Stan_glmer: yhat = alpha + x_1*beta_1 + x_2*beta_2

and substitute:

yhat_m = alpha_m + (x_1 + 1)*beta_m1 + (x_2 + 1)*beta_m1

yhat_m = alpha_m + x_1*beta_m1 + beta_m1 + x_2*beta_m2 + beta_m2

yhat_m = alpha_m + beta_m1 + beta_m2 + x_1*beta_m1 + x_2*beta_m2

If we assume that yhat_m ~= yhat, beta_m1 ~= beta_1, and beta_m2 ~= beta_2... Then

alpha = alpha_m + beta_m1 + beta_m2

So, I would expect the stan_glmer alpha (-1.7) to be close to the hand-coded Stan alpha + both betas (-3.2 + 1.7 - 0.1).

Which indeed it is (-1.6).

---

If you furthermore update your Stan data to scale these covariates as {0,1} instead of {1,2}:

guI_data2 <- list(g=as.integer(guImmun$mom),
                y=guImmun$immun,
                x=data.frame(guImmun$kid2p == "Y", guImmun$mom25p == "Y"),
                N=nrow(guImmun),
                K=2,
                M=nlevels(guImmun$mom))
ranIntFit2 <- stan(file = "first_model.stan", data = guI_data2,
                  iter = 500, chains = 1)

And look at the output:

> summary(ranIntFit2, pars = c('alpha', 'beta'))
$summary
              mean     se_mean        sd       2.5%        25%        50%           75%      97.5%     n_eff      Rhat
alpha   -1.5110714 0.022982199 0.1903571 -1.8974997 -1.6318370 -1.5038593 -1.3861628334 -1.1729671  68.60488 1.0505237
beta[1]  1.5224756 0.025017739 0.1737332  1.2260666  1.4058789  1.5118314  1.6492158203  1.8673450  48.22471 1.0592955
beta[2] -0.1206084 0.009410305 0.1640406 -0.4267987 -0.2368855 -0.1267984 -0.0003187197  0.1894375 303.87510 0.9964177

You can confirm for yourself that you are in the right ballpark.

After this, the differences between your model and the stan_glmer will come down to the priors, the parameterization of hierarchical parameters, the sampling quality, etc.

Aside: there are a number of ways that categorical covariates can be coded into a model.matrix, each for a specific interpretation of the effect parameters. The models are usually equivalent, meaning that one can convert from one parameterization to another using linear transformations of effects as done above.