Hello,

I am trying to run a hierarchical beta regression using the CarTask dataset found in the betareg package. I have never ran a hierarchical beta model before, so I was surprised that it sampled correctly on the first run. To make sure I was running the model correctly, I decided to check it using rstanarm, but to my surprise, I kept getting warning errors on the number of divergent transitions I was having. I attempted to solve the problem by using 10,000 iterations, but that did not fix the problem. The first question I have is if my Stan code is correct for the model below. The second is why is rstanarm having so many issues? I am not complaining; I like having an excuse to use pure Stan rather than rstanarm, but I would have figured this model would have converged fairly easily in rstanarm.

Thanks.

Model:

probability ~ NFCCscale + (1|task)

Stan:

```
stan_beta_hier <- "
data {
int<lower=1> N; // sample size
vector<lower=0,upper=1>[N] y; // response
int <lower = 0> num_fix; //number of fixed effects including intercept
matrix[N,num_fix] fix_mat;
int <lower = 0> num_random; //number of subjects
matrix[N,num_random] ran_mat;
}
parameters {
vector[num_fix] beta_fix; // coefficients for beta_fix
vector[num_random] beta_ran; // coefficients for beta_ran
real<lower = 0> phi; // dispersion
}
model {
// model calculations
vector[N] LP; // linear predictor
vector[N] mu; // transformed linear predictor
vector[N] A; // parameter for beta distn
vector[N] B; // parameter for beta distn
LP = fix_mat*beta_fix + ran_mat*beta_ran;
for (i in 1:N) {
mu[i] = inv_logit(LP[i]);
}
A = mu * phi;
B = (1.0 - mu) * phi;
// priors
beta_ran ~ normal(0, 5);
beta_fix ~ normal(0, 5);
phi ~ cauchy(0, 5); // different options for phi
// likelihood
y ~ beta(A, B);
}
"
# being lazy with my design matricies
m0 <- lmer(probability ~ NFCCscale + (1|task), data = CarTask)
x = getME(m0, "X")
z = as.matrix(getME(m0, "Z"))
# data
dat = list(N = length(CarTask$task),
y = CarTask$probability,
num_fix = dim(x)[2],
num_random = dim(z)[2],
fix_mat = x,
ran_mat = z)
# model/sampling
mod <- stan_model(model_code = stan_beta_hier)
fit <- sampling(mod,
data = dat,
chains = 4, iter = 2000, verbose = FALSE)
```

rstanarm:

```
library(rstanarm)
m1 <- stan_glmer(probability ~ NFCCscale + (1|task),
data = CarTask,
family = mgcv::betar(),
chains = 4, iter = 5000)
summary(m1, digits = 4)
```