Hierarchical Model for ragged/ unbalanced data (in STAN)

Question

(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?