Hierarchical Beta Regression Working in Stan but not RStanarm

Modeling
,
1

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

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

3

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