Hello. I am having trouble specifying a non-nested multi-level model. I want to be up front and tell everybody that this is for part of my disseration. I feel I have a good idea of what I’m trying to do theoretically, but unfortunately, my model output appears much worse than I would expect from the data. I am getting various warnings regarding Rhat, tree depth, and adapt_delta. I believe I am coding the model incorrectly. The response is a gamma distribution, because it must be positive. I have a conceptual DAG that I am using to base the code from, but I am not sure how to attach that here. Hopefully, you can find it at this link:
I can put the posterior formula here it seems (sorry, but the align commands aren’t working as expected, so it may look a little messy):
[\alpha_0, \alpha_1, \beta_1, \beta_2, \beta_3, \gamma_0, \gamma_1, \gamma_2, \sigma^2, \varsigma_m^2, \varsigma_s^2| SPM_{mse}] \propto\
\prod_{e=1}^{n} \prod_{m=1}^{M} \prod_{s=1}^{S} \text{gamma}(SPM_{mse}|g(\alpha_0, \gamma_0, \beta_1, \beta_2, \beta_3), \sigma^2)\\
\times \text{normal}(\alpha_0|\alpha_1, \varsigma_m^2)\ \text{normal}(\gamma_0|h(\gamma_1, \gamma_2), \varsigma_s^2) \\
\times \text{gamma}(\alpha_1 | 0.001, 0.001)\ \text{gamma}(\beta_1 | 0.001, 0.001)\\
\times \text{gamma}(\beta_2 | 0.001, 0.001)\ \text{normal}(\beta_3 | 0, 100)\\
\times \text{gamma}(\gamma_1 | 0.001, 0.001)\ \text{gamma}(\gamma_2 | 0.001, 0.001)\\
\times \text{inverse gamma}(\sigma^2 | 0.001, 0.001)\\
\times \text{inverse gamma}(\varsigma_m^2 | 0.001, 0.001)\\
\times \text{inverse gamma}(\varsigma_s^2 | 0.001, 0.001)
This is my model code:
data {
int<lower=0> N; //rows of data
int<lower=0> D_e; // number of event covariates
int<lower=0> D_m; // number of month covariates
int<lower=0> D_s; // number of station covariates
int<lower=0> L_m; // number of months
int<lower=0> L_s; // number of stations
matrix[N,D_e] x_e; // event design matrix (predictors)
matrix[N,D_m] x_m; // month design matrix (predictors)
matrix[N,D_s] x_s; // station design matrix (predictors)
vector[N] y; // response (outcomes)
int<lower=1, upper=L_m> m_groups[N]; // month IDs
int<lower=1, upper=L_s> s_groups[N]; // station IDs
}
parameters {
real alpha_0; // month intercept
real gamma_0; // station intercept
vector[D_e] beta; // event slope parameters
vector[D_m] alpha; // month slope parameters
vector[D_s] gamma; // station slope parameters
real<lower=0, upper=10> sigma_sq; // population variance
vector[L_m] u_m; // month random effects
vector[L_s] u_s; // station random effects
real<lower=0, upper=10> sigma_u_m_sq; // variance of station effects
real<lower=0, upper=10> sigma_u_s_sq; // variance of month effects
}
transformed parameters {
vector[N] mu_e; // individual event means
vector[N] mu_m; // individual month means
vector[N] mu_s; // individual station means
vector[L_m] alpha_0j; // varying month intercept
vector[L_s] gamma_0j; // varying station intecept
alpha_0j = alpha_0 + u_m;
gamma_0j = gamma_0 + u_s;
mu_m = alpha_0j[m_groups] + x_m * alpha;
mu_s = gamma_0j[s_groups] + x_s * gamma;
mu_e = alpha_0j[m_groups] + gamma_0j[s_groups] + x_e * beta;
}
model {
// Priors
alpha_0 ~ normal(0, 10);
gamma_0 ~ normal(0, 10);
alpha[1] ~ gamma(0.001, 0.001); // Prior for discharge
gamma[1] ~ gamma(0.001, 0.001); // Prior for distance
gamma[2] ~ gamma(0.001, 0.001); // Prior for depth
beta[1] ~ gamma(0.001, 0.001); // Prior for wind stress
beta[2] ~ gamma(0.001, 0.001); // Prior for PAR
beta[3] ~ normal(0, 10); // Prior for temperature
sigma_sq ~ inv_gamma(0.001, 0.001);
sigma_u_m_sq ~ inv_gamma(0.001, 0.001);
sigma_u_s_sq ~ inv_gamma(0.001, 0.001);
// Random effects distribution
u_m ~ normal(0, sqrt(sigma_u_m_sq));
u_s ~ normal(0, sqrt(sigma_u_s_sq));
// Data Model
y ~ normal(mu_m, sqrt(sigma_u_m_sq));
y ~ normal(mu_s, sqrt(sigma_u_s_sq));
y ~ gamma(mu_e, sqrt(sigma_sq));
}
Hopefully, the data can be downloaded here.
Does anyone see what I might be doing wrong?