# Hierarchical Beta Regression Working in Stan but not RStanarm

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:

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

x = getME(m0, "X")
z = as.matrix(getME(m0, "Z"))

# data
num_fix = dim(x),
num_random = dim(z),
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),
family = mgcv::betar(),
chains = 4, iter = 5000)
summary(m1, digits = 4)
``````

Should have mentioned, but the Stan code is adapted from Michael Clark

I can’t help with rstanarm (pinging @bgoodri, @Jonah), but I can suggest some cleanup for the Stan code:

``````beta_ran ~ normal(0, 5);
beta_fix ~ normal(0, 5);
phi ~ cauchy(0, 5);   // Gelman suggests a Pareto distribution here, but Cauchy is also very long-tailed
vector[N] mu = inv_logit(fix_mat * beta_fix + ran_mat * beta_ran);
y ~ beta_proportion(mu, phi);
``````

We use that beta parameterization so much we built it in as `beta_proportion`.

1 Like