Hi everyone,
I’m trying to build a model that propagates uncertainty across several linked components (possibly something like a Bayesian path model / SEM). Below is a simplified version of the full model, which is behaving in some unexpected ways.
Data structure:
Site (gridref)
Site-year (gridref-year)
Transects within each site-year
I want to test how species richness predicts total abundance, separating:
-
Between-site differences in richness
-
Within-site (year-to-year) variation in richness
-
The interaction between these two.
In a later stage of the model, I want to propagate these estimates forward to predict other responses (tick abundance and Borrelia prevalence), but I’ve removed those parts here for simplicity.
First, I decompose richness into site and site-year components using random effects. I then use estimated site-level richness, site-year deviation from that richness, and their interaction as predictors of abundance, including their interaction.
The model runs without obvious diagnostics issues (no divergences, good Rhat), but there’s a posterior ridge between the intercepts for the richness and abundance components and when I add the downstream response variables (ticks and Borrelia prevalence), the geometry becomes more problematic and sampling gets ugly.
I tried to scale the site-level and site-year level richness estimates in the transformed parameters block but that also makes things worse!
What might explain this posterior ridge between the intercepts?
Is this a known identifiability issue when propagating latent predictors between model components?
Is there a better modeling strategy I should consider?
Any advice would be greatly appreciated.
data {
int n_bbs_tr_v ; // number of transect visit's
int gr_yrs ; // number of gridref year's
int grs ; // number of grid refferences (sites)
int gr_for_gr_yr[gr_yrs] ; // gridref to which each gridref year belongs
int gr_yr_for_tr_v[n_bbs_tr_v] ; // gridref year to which each transect visit belongs (only those sampled in the bbs)
int rich [n_bbs_tr_v] ; // observed species richness (transect level)
real visit[n_bbs_tr_v] ; // was the visit early or late in the season (--0.5 or 0.5)
int babnd [n_bbs_tr_v] ; // total abundance of birds
}
parameters{
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
// MODEL A - estimating gridref-year richness ////////////////////////////////
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
real rich_alpha ; // global mean richness
real rich_beta_visit ; // visit time of year coefficient
vector[gr_yrs] rich_gr_yr_raw; // raw random effect
real<lower = 0> rich_gr_yr_sigma; // gr_yr random effect scale (multiple transects within a gr_yr)
vector[grs] rich_gr_raw; // gr raw random effect
real<lower = 0> rich_gr_sigma; // gr random effect scale (multiple gr_yrs (years) within a gr (a site))
// //AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
// // MODEL B2 - estimating gridref-year bird abundance /////////////////////////
// //AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
real babnd_alpha ; // community competence global intercept
real babnd_beta_visit ;
real babnd_rich_mean ; // effect of gr (site) level richness (latent) on bird abundance
real babnd_rich ; // effect of deviation from gr (site level mean - latent) on bird abundance
real babnd_rm_int ; // an interaction between the two (latent * latent)
vector[gr_yrs] babnd_gr_yr_raw ; // the same random effect structure as above to capture variation not explained by richness
real<lower = 0> babnd_sigma_gr_yr ; // i.e site level processes other than richness could influence total abundance
vector[grs] babnd_gr_raw ; // (richness is just a proxy for some unobserved process i.e habitat degredation or something like that ...)
real<lower = 0> babnd_gr_sigma ;
}
transformed parameters{
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
// MODEL A - estimating gridref-year richness (deviation) and mean richness //
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
vector[grs] rich_gr = rich_gr_raw * rich_gr_sigma ;
vector[gr_yrs] rgr_yr =
rich_gr_yr_raw * rich_gr_yr_sigma;
vector[gr_yrs] rich_beta_gr_yr =
rgr_yr +
rich_gr[gr_for_gr_yr] ;
// extract mean gridref richness estimates
vector[gr_yrs] r_gr_mean;
for (i in 1:gr_yrs) {
int g = gr_for_gr_yr[i];
r_gr_mean[i] = rich_gr[g];
}
vector[gr_yrs] r_int = rgr_yr .* r_gr_mean ; // pre compute interaction term
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
// MODEL B2 - estimating gridref-year bird abundance /////////////////////////
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
vector[grs] babnd_gr ;
babnd_gr = babnd_gr_raw * babnd_gr_sigma ;
vector[gr_yrs] babnd_beta_gr_yr = // site-year level bird abundance is a product of ...
babnd_gr_yr_raw * babnd_sigma_gr_yr + // a site-year random effect
babnd_gr[gr_for_gr_yr] + // a site random effect
babnd_rich_mean * r_gr_mean + // site level mean richness
babnd_rich * rgr_yr + // deviance from site-level mean richness
babnd_rm_int * r_int; // and interaction between the two
}
model {
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
// MODEL A - estimating gridref-year richness ////////////////////////////////
//AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
rich_gr_yr_raw ~ normal(0, 1) ;
rich_gr_yr_sigma ~ exponential(2) ;
rich_gr_raw ~ normal(0, 1) ;
rich_gr_sigma ~ exponential(2) ;
for(i in 1:n_bbs_tr_v) {
real rich_lambda = exp(
rich_alpha +
rich_beta_gr_yr[gr_yr_for_tr_v[i]] +
rich_beta_visit * visit[i]
) ;
rich[i] ~ poisson(rich_lambda) ;
};
rich_alpha ~ normal(0, 0.5) ;
rich_beta_visit ~ normal(0, 0.5) ;
// //AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
// // MODEL B2 - estimating gridref-year bird abundance /////////////////////////
// //AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA//
babnd_alpha ~ normal(0, 0.5);
babnd_beta_visit ~ normal(0, 0.5);
babnd_rich ~ normal(0, 0.5) ;
babnd_rich_mean ~ normal(0, 0.5) ;
babnd_rm_int ~ normal(0, 0.5) ;
babnd_gr_raw ~ normal(0, 1) ;
babnd_gr_sigma ~ exponential(2) ;
babnd_gr_yr_raw ~ normal(0, 1) ;
babnd_sigma_gr_yr ~ exponential(2) ;
// Transect-level model
for(i in 1:n_bbs_tr_v) {
real babnd_mu = exp(
babnd_alpha +
babnd_beta_visit * visit[i] +
babnd_beta_gr_yr[gr_yr_for_tr_v[i]]
);
babnd[i] ~ poisson(babnd_mu);
}
}
