I collected cells from 28 mice that received 6 different treatments (3-5 per group).

The number of cells I collected is highly variable across treatments; one group only has a couple hundred per animal, and some have a few million.

I’m trying to build a multi-level model with a random effect for each animal.

So, something like:

Cells \sim Poisson (\lambda)

log(\lambda) = mu0 + bT[Treatment] + bM[Mice]

The thing is, I cannot get the model to converge.

My model looks like

```
data{
int Cells[28];
int Treatment[28];
int Mice[28];
}
parameters{
real<lower=0> mu0;
vector[28] bM;
vector[6] bT;
}
model{
vector[28] mu;
bT ~ normal( 0 , 6 );
bM ~ normal( 0 , 2 );
mu0 ~ gamma( 10/1 , 1/1 );
for ( i in 1:28 ) {
mu[i] = mu0 + bM[Mice[i]] + bT[Treatment[i]];
mu[i] = exp(mu[i]);
}
Cells ~ poisson( mu );
}
```

Data from dput

```
structure(list(Cells = c(6227606, 5844779, 1338474, 10365562,
6349543, 9831545, 10968084, 15513791, 2564861, 6725653, 1609213,
716699, 991468, 822666, 521778, 11940, 345215, 42238, 7612, 10746,
448, 448, 149, 597, 149, 14477, 1492, 7910), Mice = 1:28, Treatment = c(1,
1, 1, 1, 1, 2, 2, 2, 2, 2, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 4, 4,
4, 4, 4, 3, 3, 3)), row.names = c(NA, -28L), class = c("tbl_df",
"tbl", "data.frame"))
```

I think my priors are justifiable.

We typically expect to collect 0.5-20 million cells per mouse, so I think `mu0~gamma(10, 1)`

is reasonable.

I think I am having a hard time parameterizing `bT`

. The distribution is clearly not normal. Some treatment does nothing while some treatment kills almost every cell. I’ve tried to use `bT~cauchy(0, 2)`

for heavier tails but to no avail.

Looking for suggestions to help me (1) reparameterize my model completely and/or (2) pick better priors.

I am a biologist by training, so please also let me know if I’m doing anything that is obviously incorrect.

Many thanks in advance.

Run everything with

```
"
data{
int Cells[28];
int Treatment[28];
int Mice[28];
}
parameters{
real<lower=0> mu0;
vector[28] bM;
vector[6] bT;
}
model{
vector[28] mu;
bT ~ normal( 0 , 6 );
bM ~ normal( 0 , 2 );
mu0 ~ gamma( 10/1 , 1/1 );
for ( i in 1:28 ) {
mu[i] = mu0 + bM[Mice[i]] + bT[Treatment[i]];
mu[i] = exp(mu[i]);
}
Cells ~ poisson( mu );
}" -> stan_code
structure(list(Cells = c(6227606, 5844779, 1338474, 10365562,
6349543, 9831545, 10968084, 15513791, 2564861, 6725653, 1609213,
716699, 991468, 822666, 521778, 11940, 345215, 42238, 7612, 10746,
448, 448, 149, 597, 149, 14477, 1492, 7910), Mice = 1:28, Treatment = c(1,
1, 1, 1, 1, 2, 2, 2, 2, 2, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 4, 4,
4, 4, 4, 3, 3, 3)), row.names = c(NA, -28L), class = c("tbl_df",
"tbl", "data.frame")) -> dd
rstan::stan(model_code = stan_code, data = list(Cells = dd$Cells, Mice = dd$Mice, Treatment = dd$Treatment))
```