I am an evolutionary biologist trying to fit a model (`brm.count`

) in which:

- the response variable is the count of a phenotypic structure (e.g., toes) that a species possesses (
`trait.count`

); - the population-level effects are two continuous variables and their interaction (
`continuous1`

and`continuous2`

); - intercepts can vary over the group-level effect
`species`

, and species are correlated according to the`phylogeny`

covariance matrix.

Since the response variable is count data, I opted to choose a Poisson distribution. There is slight overdispersion in the data (index of dispersion, \sigma^2/\mu = 1.1) and 78% of the observations are zeroes. This is how the distribution of counts (`trait.count`

) looks like:

The code that I used to run the model is as follows:

```
brm_count <- brm(trait.count ~ continuous1*continuous2 + (1|gr(species, cov = A)),
data = data,
family = poisson(),
data2 = list(A = phylogeny),
chains = 2,
cores = 2,
iter = 4000)
```

When I run a posterior predictive check on the output of this model, hereâs what I get:

```
pp_check(brm_count, type="bars", nsamples = 100)
```

Besides underestimating the number of ones, the model seems to be overestimating the number of zeroes. This is further evidenced by checking the proportion of zeroes in samples from the posterior predictive distribution, compared to the observed proportion.

```
prop_zero <- function(x) {sum(x == 0)/length(x)}
ppc_stat(y = data$trait.count, yrep = posterior_predict(brm_count, draws = 1000), stat="prop_zero")
```

But what I found the weirdest is that the model also overestimates dispersion!

```
dispersion <- function(x) {var(x)/mean(x)}
ppc_stat(y = data$trait.count, yrep = posterior_predict(brm_count, draws = 1000), stat="dispersion")
```

My question is: how come the dispersion of my data (1.1) is inflated to ~2 by a model that assumes no overdispersion (index = 1)? Shouldnât I expect the opposite: for dispersion to be underestimated? Am I missing something?

(By the way, I also fitted this model with zero-inflated Poisson, negative binomial and Poisson with observation-level random effects `(1|obs)`

, and the issue persists.)

- Operating System: macOS 10.12.6
- brms Version: 2.13.0

Thank you so much for any help!

UPDATE: I ran `ppc_stat()`

with variance and mean to check the source of the overdispersion. Turns out the mean is pretty accurately estimated, while the variance is overestimated (0.3 in the observed data, ~0.6 in the posterior predictive samples).