I’m relatively new to Stan: I’ve done a few basic tutorials and have used Staninterfacing packages in R
, such as brms
, etc. A while ago, I asked a question about a betabinomialbased model for ant preferences on stackexchange and (perhaps understandably) didn’t get any answers. I tried my best to do as much as possible with what I could find online, and after that, came here. I’ve refined my questions, so here goes.
General goal: I want to know whether ants of a certain species have individual preferences in terms of which types of food they bring back to the nest, or if they just get want the colony needs.
Setup: Suppose there are N ant colonies of a certain species, where colony i (where i \in \big\{1,2,..., N\big\} ) has N_i worker ants. Ant j_i (where j_i \in \big\{1,2,..., N_i\big\} ) goes on N_{ji} trips, where on each trip it decides to bring back either Apples or Bananas (a binary choice) to the colony.
Given a record of individual ants’ trips (across different colonies), I want to see how likely it is that ants of a certain species have specialized preferences (i.e., some ants get Apples a lot, while others go for Bananas mostly) or whether they all just get whatever their colony wants in general.
My approach: I’m assuming that the number of times Apples are brought back by an ant j_i is sampled from Binomial(N_{ji}, p_{ji}), where p_{ji} is ant j_i's ‘preference’ toward Apples as opposed to Bananas. I’m assuming (for one of my hypotheses) that p_{ji} \sim Beta(\mu_i, t_i), (parameterized with \mu and \nu) where \mu_i is the colony’s ‘preference’ for Apples, and t_i is the total number of trips taken by ants in colony i.
In order to make a model that can use data from multiple colonies simultaneously, I connect the colonies’ preferences (the \mu_i\mathrm{s}) via a hyperparameter distribution Beta(\mu_{species}, \mathrm{var}) (parameterized with \mu and variance) and where \mu_{species} is the overall species preference for Apples and \mathrm{var} is the variance.
Alternative hypothesis: The alternative hypothesis—i.e., if individual ants don’t have preferences—would be that every time an ant goes on a trip, it basically flips a coin weighted according to the colony’s preference \mu_i. I.e., the number of times ant j_i brings back Apples is sampled from Binom(N_{ji},\mu_i) instead.
Additional (optional) goal: If possible, I would like to be able to assign a particular value to a species to indicate how “specialized” its ants are, something like the “Polarization Index” mentioned in this article by Rungie, Brown, Laurent & Rudrapatna (2005). Essentially, since the highest variance of beta distribution is bounded, the Polarization Index is “a standardized form of the variance: the variance divided by its upper bound”, where 0 is when all the decision makers (i.e., ants) are making decisions at random, and 1 is when they all are completely constant.
MY PROBLEMS:

Does my approach make sense? I’m particularly worried about how my model incorporates data from multiple colonies—I want to test hypotheses about the species, and I’m not sure how my model can do this.

Is my Stan model (below) written correctly? I’ve included my best attempt at a Stan model of the generative process I’ve just described, as well as sample data that might go in. Any advice on making it better would be appreciated!

How do I hypothesis test with this model? I’m pretty much in the dark here. My instinct would be to do a Bayes factor with the probability that the current model is true given the data vs. the alternative hypothesis model is true given the data. Is that kosher? How do I do that?

(Optional) Can I get a ‘specialization’ metric from the data? The ideal scenario for me would be to be able to take in ant trip data and output a metric of how “specialized” an ant species’ workers are. Preferably, this would be something like the “Polarization Index” I mentioned earlier, but I’m not sure if that fully captures the intuition I’m going for and I’m not sure how I would aggregate data from multiple colonies into a single value.
Any help on anything would be appreciated—I’m at the very limits of my understanding for this project!
Below are what I’ve written for my model and some toy data to show how the data would be structured:

Here’s some toy data in the ‘expected’ format:
[3, // number_of_colonies
[2, 1, 3], // ants_per_colony
6, // number_of_ants (technically, redundant)
[90, 17, 4, 1, 100, 33], // trips_per_ant
[80, 0, 2, 1, 25, 33] // apples_per_ant
]
The ants that made 90
and 17
trips belong to the first colony, the one that made 4
trips to the second, and the 1
, 100
, and 33
trip ants belong to the last colony (as indicated by ants_per_colony
).
Here’s my Stan code:
data {
int<lower=1> number_of_colonies; // # total number of ant colonies
int<lower=1> ants_per_colony[number_of_colonies]; // # total number of ants per individual colony
int<lower=1> number_of_ants; // # total number of ants
int<lower=1> trips_per_ant[number_of_ants]; // # number of trips per ant, in order of colony
int<lower=0> apples_per_ant[number_of_ants]; // # total number of apples returned per ant, in order of colony
}
// # This is just a way of assigning a colony to each ant via ants_per_colony and trips_per_colony:
// # Essentially, it just gets the start and stop indices in the ant list of each colony
transformed data {
int<lower=1> ant_index_start[number_of_colonies];
int<lower=1> ant_index_end[number_of_colonies];
vector<lower=1>[number_of_colonies] trips_per_colony;
if (number_of_colonies == 1) {
ant_index_start[1] = 1;
ant_index_end[1] = ants_per_colony[1];
} else {
int c;
c = 1;
for (x in 1:number_of_colonies) {
ant_index_start[x] = c;
c = c + ants_per_colony[x];
ant_index_end[x] = c  1;
}
}
for (x in 1:number_of_colonies) {
trips_per_colony[x] = sum(trips_per_ant[ant_index_start[x]:ant_index_end[x]]);
}
}
parameters {
real<lower=0,upper=1> hyper_mean;
real<lower=0> hyper_variance;
vector<lower=0,upper=1>[number_of_colonies] colony_mean;
vector<lower=0,upper=1>[number_of_ants] ant_preferences;
}
// # Reparametrizing parameters from things I care about to alphas and betas
transformed parameters {
vector<lower=0>[number_of_colonies] colony_alpha;
vector<lower=0>[number_of_colonies] colony_beta;
real<lower=0> colony_hyper_alpha;
real<lower=0> colony_hyper_beta;
colony_alpha = colony_mean .* trips_per_colony;
colony_beta = (1  colony_mean) .* trips_per_colony;
colony_hyper_alpha = hyper_mean * (hyper_mean * (1  hyper_mean) / hyper_variance  1);
colony_hyper_beta = (1  hyper_mean) * (hyper_mean * (1  hyper_mean) / hyper_variance  1);
}
model {
// Colonies' mean preferences for apples are connected via hyperparameters
colony_mean ~ beta(colony_hyper_alpha, colony_hyper_beta);
for (i in 1:number_of_colonies) {
for (j in ant_index_start[i]:ant_index_end[i]) {
ant_preferences[j] ~ beta(colony_alpha[i], colony_beta[i]);
}
}
apples_per_ant ~ binomial(trips_per_ant, ant_preferences);
}