functions { /** inverse ilr transformation of a vector x, using the inverse of the transpose of the V matrix of the ilr (tVinv) */ vector ilrinv(matrix tVinv, vector x, int ntaxa) { vector[ntaxa] z; vector[ntaxa] y; z = exp(tVinv * x); y = z / sum(z); return y; } } data { int ntaxa; //number of taxa int npanels; //number of panels matrix[ntaxa, ntaxa - 1] tVinv; //back-transformation matrix for ilr transformation int counts[npanels, ntaxa]; //observed counts int n[npanels]; //total number of points on each panel } transformed data { int s = ntaxa - 1; vector[s] mu_; vector[s] x_[npanels]; //prior predicted logratio coordinates vector[ntaxa] rho_[npanels]; //prior predicted relative abundances int y[npanels, ntaxa]; //prior predicted counts for(i in 1:s){ mu_[i] = normal_rng(0, 2); } for(i in 1:npanels) { x_[i] = mu_; rho_[i] = ilrinv(tVinv, x_[i], ntaxa); y[i] = multinomial_rng(rho_[i], n[i]); // realizations from prior predictive distribution } } parameters { vector[s] mu; //intercept } transformed parameters { vector[s] x[npanels]; //predicted logratio coordinates vector[ntaxa] rho[npanels]; //predicted relative abundances for(i in 1:npanels){ x[i] = mu; rho[i] = ilrinv(tVinv, x[i], ntaxa); } } model { for(i in 1:npanels) { counts[i] ~ multinomial(rho[i]); // observation model } mu ~ normal(0, 2); } generated quantities{ vector [npanels] log_lik; vector[s] pars_; int y_[npanels, ntaxa] = y; int ranks_[s]; for(i in 1:s){ ranks_[i] = mu[i] > mu_[i]; pars_[i] = mu_[i]; } for (j in 1:npanels) log_lik[j] = multinomial_lpmf(counts[j] | rho[j]); //log likelihood for WAIC }