Hi all,

I am having issues using beta_binomial. I started with trying to fit m12.1 in Statistical Rethinking (beta_binomial fit to classic UC Berkeley admissions data). The ulam generated stancode follows.

```
data{
int N[12];
int A[12];
int gid[12];
}
parameters{
vector[2] a;
real<lower=0> phi;
}
transformed parameters{
real theta;
theta = phi + 2;
}
model{
vector[12] pbar;
phi ~ exponential( 1 );
a ~ normal( 0 , 1.5 );
for ( i in 1:12 ) {
pbar[i] = a[gid[i]];
pbar[i] = inv_logit(pbar[i]);
}
A ~ beta_binomial( N , pbar*theta , (1-pbar)*theta );
}
```

Warnings such as maximum treedepth, large R-hat, low ESS etc. are seen. I tried adjusting the priors with no luck.

Separating the beta_binomial into hierarchical pbar ~ beta(alpha,beta) and admits ~ binomial(apps, pbar) worked (code below). Moreover, the posterior samples (of a, phi, theta) had very similar summary statistics to the ones in Statistical Rethinking (R code 12.3 on pg. 372).

```
data{
int N;
int admits[N];
int apps[N];
int gid[N];
}
parameters{
vector[2] a;
real<lower=0> phi;
vector<lower=0, upper=1>[N] pbar;
}
transformed parameters{
vector[N] mu;
vector[N] alpha;
vector[N] beta;
real<lower=2> theta;
theta = phi + 2;
for (i in 1:N){
mu[i] = inv_logit(a[gid[i]]);
}
alpha = mu * theta;
beta = (1-mu) * theta;
}
model{
pbar ~ beta(alpha, beta);
admits ~ binomial(apps, pbar);
a ~ normal(0, 1.5);
phi ~ exponential(1);
}
```

Should the separation into beta then binomial not be equivalent to using beta_binomial in one step? (some naming and other slight changes to practice my stan codingā¦). I also tried more or less exactly the above code but recombining into beta_binomial, which makes the beforeseen errors resurface.

```
data{
int N;
int admits[N];
int apps[N];
int gid[N];
}
parameters{
vector[2] a;
real<lower=0> phi;
vector<lower=0, upper=1>[N] pbar;
}
transformed parameters{
vector[N] mu;
vector[N] alpha;
vector[N] beta;
real<lower=2> theta;
theta = phi + 2;
for (i in 1:N){
mu[i] = inv_logit(a[gid[i]]);
}
alpha = mu * theta;
beta = (1-mu) * theta;
}
model{
admits ~ beta_binomial(apps, alpha, beta);
a ~ normal(0, 1.5);
phi ~ exponential(1);
}
```

Would appreciate any illumination into how to use beta_binomial better, or if my assumption of the equivalent model specifications is wrong. Thank you!