# Convergence for a Bernoulli logit regression is weird

I tried to fit a mixed hermetical model using rstan version 2.26.22 (Stan version 2.26.1). The beta part convergence to a rational result, but the Bernoulli one seems convergence to a total irrational results. I used old stan one year ago, and at that time, the results seem to be rational.

The data generation looks like this:

 x<-matrix(NA,20,100)
z<-matrix(NA,20,100)

for (i in 1:20) {
for (j in 1:100) {
x[i,j]<-rnorm(1,7.5,1.5)
}
}

for (i in 1:20) {
for (j in 1:100) {
z[i,j]<-rnorm(1,5,2)
}
}

sigmab<-0.09164
sigmam<-0.077
b<-rnorm(20,0,sqrt(sigmab))
m<-rnorm(20,0,sqrt(sigmam))

function1 <- function(beta0,beta1){
e=matrix(NA,20,100)
mu=matrix(NA,20,100)
for (i in 1:20) {
for (j in 1:100) {
e[i,j]=exp(beta0+x[i,j]*beta1+b[i])
mu[i,j]=e[i,j]/(1+e[i,j])
}
}
mu
}
mu<-function1(1.25,0.19)

function2 <- function(alpha0,alpha1){
e=matrix(NA,20,100)
pi=matrix(NA,20,100)
for (i in 1:20) {
for (j in 1:100) {
e[i,j]=exp(alpha0 + z[i,j]*alpha1 +m[i])
pi[i,j]=e[i,j]/(1+e[i,j])
}
}
pi
}
pi<-function2(2.16,-0.7162)

w<-matrix(NA,20,100)

for (i in 1:20) {
for (j in 1:100) {
w[i,j] <- rbinom(1,1,pi[i,j])
}
}

ystar = matrix(NA,20,100)

for (i in 1:20) {
for (j in 1:100) {
ystar[i,j] <- rbeta(1,mu[i,j]*42,42*(1-mu[i,j]))
}
}

dim(ystar)

y=ifelse(w==1,w,ystar)

###############################################{SAMPLE=1/10 y}

samplex=c()
sampley=c()
samplez=c()
for (i in 1:20) {
samplei <- sample(1:100,10)
sampley = c(sampley, y[i,samplei])
samplex = c(samplex, x[i,samplei])
samplez = c(samplez, z[i,samplei])
}

areasamp  <- rep(1:20, each =10)
sampleratio <- data.frame(areasamp, samplex,sampley,samplez)
sampleratio$samplex=sampleratio$samplex-mean(sampleratio$samplex) sampleratio$samplez=sampleratio$samplez-mean(sampleratio$samplez)

databeta <- subset(sampleratio, sampley <1)
sb_list_multi <- list(N = nrow(sampleratio),
W = ifelse(sampleratio$sampley == 1, 1, 0), z = sampleratio$samplez,
County = as.numeric(as.factor(sampleratio\$areasamp)))



The Stan file looks like this:

data{
int<lower=0> N;
vector[N] z;
array[N] int<lower=0, upper=1> W;
int County[N];
}
parameters{
real<lower = 0> sigmam;
real alpha0;
real alpha1;
real u[20];
}
transformed parameters{
vector<lower=0, upper=1>[N] p;
for(i in 1:N){
p[i] =exp(alpha0 + alpha1*z[i] + u[County[i]])/(1 + exp(alpha0 + alpha1*z[i]+ u[County[i]])) 	;
}
}
model{
alpha0 ~ normal(0, 10000);
alpha1 ~ normal(0, 10000);
sigmam ~ cauchy(0, 10000);
for(i in 1:N){
W[i] ~ bernoulli(p[i]);
}

for(i in 1:20){
u[i] ~ normal(0,sigmam);
}
}



And the results look like this

         mean se_mean   sd    2.5%    25%    50%    75%  97.5% n_eff Rhat
sigmau   0.94    0.03 0.41    0.27   0.65   0.89   1.16   1.88   254 1.02
alpha0  -1.60    0.02 0.34   -2.34  -1.80  -1.56  -1.37  -1.03   427 1.01
alpha1  -0.70    0.01 0.14   -1.01  -0.79  -0.69  -0.60  -0.43   787 1.00


Where \alpha_0 should be positive and \sigma_m is way higher than setting