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)
```

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.

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