data {
int N; // number of observations in proportion data
int J; // number of mosquito species
int host_counts[N]; // observed host abundance data
int bites[J, N]; // bite data of length N
vector[N] flat_alpha_p; // Vector that could be informed by some other data to weight previous "observations" that form the prior
vector[N] prev_obs_p; // The previous observations that can get scaled by flat_alpha_p (for now flat_alpha_p is all 1s)
real sigma_pb_a; // first of gamma parameters for hyper-prior that establishes the variation in a mosquitoes biting preference
real sigma_pb_b; // second of gamma parameters for hyper-prior that establishes the variation in a mosquitoes biting preference
}
transformed data {
vector[N] alpha_p;
alpha_p = flat_alpha_p .* prev_obs_p; // sets up strength of prior
}
parameters {
// host proportions
simplex[N] theta_p; // vector of length = K dimensions that sums to 1
// biting preference
matrix[J, N] bite_pref; // underlying latent parameter of interest. Intrinsic preference of a mosquito on a given host
// Variation among mosquitoes in their biting preference distribution
// Goal is to have this large enough that it allows some hosts to be generalists and some hosts to be specialists
vector[J] prev_obs_bA;
vector[J] prev_obs_bB;
}
model {
// host data. host_counts is the raw data that informs theta_p (which gets mixed with the prior on host proportions) to inform the bite data
theta_p ~ dirichlet(alpha_p); // prior on host proportions
host_counts ~ multinomial(theta_p); // Actual observed abundance is some multinomial draw from the underlying proportions
// biting preference
for (j in 1:J) {
// Currently the prior on each mosquitoes preference across hosts is gamma distributed, where the parameters are drawn from the hyper
prev_obs_bA[j] ~ gamma(sigma_pb_a, sigma_pb_a);
prev_obs_bB[j] ~ gamma(sigma_pb_b, sigma_pb_b);
bite_pref[j, ] ~ gamma(prev_obs_bA[j], prev_obs_bB[j]);
}
// model fit to blood meal data
for (j in 1:J) {
bites[j, ] ~ multinomial((to_vector(theta_p) .* to_vector(bite_pref[j, ])) /
sum(to_vector(theta_p) .* to_vector(bite_pref[j, ])));
}
}