I am having trouble getting the chains to mix properly in this hierarchal model. The model can be described mathematically as:

\lambda_k \sim Normal(\mu, \sigma)

\delta_k \sim Normal(\nu, \rho)

\tau_{k,i} \sim Uniform(0,T_k)

\alpha_{k,i} \sim Poisson(m_i N_{k,i} \lambda_k \tau_{k,i})

\beta_{k,i} \sim Poisson(n_i N_{k,i} \lambda_k (T_k - \delta_k - \tau_{k,i}))

With appropriate priors on \mu, \sigma, \nu, \rho.

Briefly, I am trying to estimate the age (\tau) of a patient when copy number aberrations - e.g. going from 2 (m_i) chromosomes to 4 (n_i) - occur based on the number of mutations that occurred before (\alpha_{k,i}) and after (\beta_{k,i}) the CNA. I assume that, for each patient, mutations occur at a constant rate per base pair (\lambda_k) on a region of size N_{k,i}. One complication is that we can only detect mutations up until a given time in the past, which I have characterised by including a patient specific lag-time \delta. If I use a non-hierarchal model (i.e \lambda_k and \delta_k are the same for each patient) and apply it to simulated data, than the chains converge and I can recover the input parameters. However, introducing a hierarchal model causes the chains to fail to converge. Have I misspecified my model, or do I need to use more informative priors?

```
data {
int<lower=1> N; // # of chromosomal events
int<lower=0> K; // # of groups
int<lower=0> alpha[K, N];
int<lower=0> beta[K, N];
vector<lower=0>[K] age;
vector<lower=0>[N] region_size[K];
vector<lower=1>[N] num_strands_before[K];
vector<lower=1>[N] num_strands_after[K];
}
parameters {
vector<lower=0, upper=1>[K] t_offset_raw;
vector<lower=0, upper=1>[N] t_raw[K];
vector<lower=0>[K] rate;
real<lower=0> rate_mu;
real<lower=0> rate_sigma;
real<lower=0> t_offset_mu;
real<lower=0> t_offset_sigma;
}
transformed parameters{
vector[K] t_offset;
vector[N] t[K];
vector[N] lam1[K];
vector[N] lam2[K];
t_offset = age .* t_offset_raw;
for (k in 1:K){
t[k] = (age[k] - t_offset[k]) * t_raw[k];
lam1[k] = rate[k] * num_strands_before[k] .* region_size[k] .* t[k];
lam2[k] = rate[k] * num_strands_before[k] .* region_size[k] .* (age[k] - t[k] - t_offset[k]);
}
}
model {
t_offset_mu ~ normal(5, 3);
t_offset_sigma ~ normal(2, 3);
rate_mu ~ normal(0, 0.1);
rate_sigma ~ normal(0, 0.05);
for (k in 1:K){
alpha[k, :] ~ poisson(lam1[k]);
beta[k, :] ~ poisson(lam2[k]);
}
}
```