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