(I'm fairly new to Bayesian modelling please forgive me any minor accidents in my questions)
I'm trying to model a data set in STAN, but don't understand why I get large no. divergent transitions.
The data sets consists of data from $k=3$ different field experiments (conducted in different years) in which ($N_k = 13, 14, 16$) samples were taken. In order to be actually analyzed, each sample had to be "measured", i.e. genetically sequenced. Now, some of the samples (in total 6) were sequenced more than once (2 or 4 times), i.e. replicated and the replications vary within those samples. Due to financial reasons genetic sequencing cannot be done to get a balanced data set. At last, the data are transformed via principle component analysis and projected onto the resulting first component (which accounts for about 70% of the total variance).
I want to do Bayesian inference about differences in field experiments while taking into account the variance induced by the measurement process/ genetic sequencing.
Model and idea:
There are 3 levels: Experiment -> Measurement -> Replication, i.e.
$$ \begin{aligned} y_{ijk} & \sim Normal\left(\mu_{jk}, \tau_{k}\right) & \text{Replications} \\ \mu_{jk} & \sim Normal\left(\mu_k, \sigma_k\right) & \text{Measurements} \\ \mu_k & \sim Normal\left(\mu, \sigma\right) & \text{Experiments} \\ \mu & \sim Normal\left(\alpha, \sigma_\alpha\right) & \text{Prior} \\ \sigma, \sigma_k & \sim Normal(0, 1) & \text{Prior} \\ \alpha, \sigma_\alpha & \sim \sigma_\alpha^{-1} & \text{Hyper-Prior} \\ \tau_k & \sim \text{Inv-} \chi ^2\left(\nu_n, s^2\right) & \text{Posterior Variance} \\ \end{aligned} $$
(1) because the data is transformed via PCA and the variance maybe due to measurement error I assume it's normally distributed and constant within each field experiment ($\tau_k$).
(2) To estimate $\tau_k$ for each experiment analytically using conjugate scaled inverse chi-quare, sample directly from it's posterior and plug-in the results on the data level $y_{ijk}$. ($s_{jk}^2$ as the replication sample variance and $s_k^2$ is $max(s^2_{jk})_{j}$).
The corresponding STAN code is
int<lower=1> N[3];
real y1[N[1]];
real y2[N[2]];
real y3[N[3]];
int<lower=1> K;
int<lower=1> M; // no. experiments
int<lower=1> S[M]; // lengths of replication data for each experiment
real r[K]; //replications, ragged array
real<lower=0> sigma0; //prior replication variance
int<lower=0> nu0; //prior replication degrees of freedom
}
transformed data {
real v[M]; // mean
real s[M]; // empirical variance
real<lower=0> nu[M]; // posterior degrees of freedom
int pos = 1;
for(i in 1:M) {
// mean of each set of replications
v[i] = mean(segment(r, pos, S[i]));
// variance of each set of replications
s[i] = variance(segment(r, pos, S[i]));
// posterior degrees of freedom (inv-chi-square) for each set of replications
nu[i] = nu0 + S[i];
pos = pos + S[i];
}
}
parameters {
real<lower=0> sigma_mu1; //sigma_1
real<lower=0> sigma_mu2; //sigma_2
real<lower=0> sigma_mu3; //sigma_3
real<lower=0> sigma_theta; //sigma
real<lower=0> sigma_alpha; //prior variance for mu
real mu1a; //mu_1
real mu2a; //mu_2
real mu3a; //mu_3
real thetaa; // mu
real alpha; // prior mean for mu (mean of experiments)
real mura1; //mu_j1
real mura2; //mu_j2
real mura3; //mu_j3
real<lower=0> sigmar[M]; //tau_k
}
transformed parameters {
real theta = alpha + thetaa * sigma_alpha;
real mu1 = theta + mu1a * sigma_theta;
real mu2 = theta + mu2a * sigma_theta;
real mu3 = theta + mu3a * sigma_theta;
real mur1 = mu1 * mura1 * sigma_mu1;
real mur2 = mu2 * mura2 * sigma_mu2;
real mur3 = mu3 * mura3 * sigma_mu3;
}
model {
// replications (every data point is a replication)
y1 ~ normal(mur1, sigmar[1]);
y2 ~ normal(mur2, sigmar[2]);
y3 ~ normal(mur3, sigmar[3]);
// sample replication variances
sigmar ~ scaled_inv_chi_square(nu, s);
// measurements
// mur1 ~ normal(mu1, sigma_mu1);
// mur2 ~ normal(mu2, sigma_mu2);
// mur3 ~ normal(mu3, sigma_mu3);
mura1 ~ std_normal();
mura2 ~ std_normal();
mura3 ~ std_normal();
// experiments
mu1a ~ std_normal();
mu2a ~ std_normal();
mu3a ~ std_normal();
// prior
thetaa ~ std_normal();
// scale parameters
sigma_mu1 ~ normal(0, 1);
sigma_mu2 ~ normal(0, 1);
sigma_mu3 ~ normal(0, 1);
sigma_theta ~ normal(0, 1);
sigma_alpha ~ normal(0, 1);
}
Unfortunately, this doesn't work because it results in large number of divergent transitions. When I model it without multiple replications, then the model runs just fine (0 divergent transitions). My guess is that it's hard for STAN to sample for single measurements in an additional Replication-level, but I don't know what to do about it.
Anyone any idea or can point me in the right direction?