Poisson model overestimating dispersion and inflating zeroes

[jocateme]

@martinmodrak

I would also expect that if you broke down the figure in post #5 by families or similar, you would see additional discrepancies.

You’re absolutely right; here’s how PP checks look for the four orders when no covariance matrix is added:

@martinmodrak

Also a wild guess: a Poisson hurdle model ( hurdle_poisson ) might also be a good fit as it treats zeroes as completely separate from non-zeroes, so it should allow to have low number of zeroes (as in Galliformes) and a Poisson-like tail… (unlike zero-inflated poisson which - as the name suggests - can’t have fewer zeroes than Poisson).

Just fitted a Poisson hurdle model. As you suspected, it works like a charm to fit zeroes, but the rest of the predictive distribution gets really weird, with super high values (sometimes as high as in the 1000s) predicted.

Just a brief comment on other attempted distributions:

• Zero-inflated Poisson produces almost the same estimates and PP checks as regular Poisson because of an extremely low zi (0.01; i.e., almost all of the zeroes being explained by the Poisson);
• Negative binomial also produces similar PP checks but with higher maximum values;
• I can’t seem to get the COM Poisson model converging; not sure why!

@martinmodrak

Finally, you can always fit multiple models and if the important inferences stay the same, you don’t have to worry much which of them is “most” appropriate

I believe this will be the best solution at the end of the day. Apparently, no single (easily implementable) distribution is perfect for my data, but if multiple reasonable distributions yield similar inferences, there is enough indication that they are reliable.

Thank you for all the help, @martinmodrak and @cmcd!

[martinmodrak]

Good work there. A few notes:

@jocateme

I can’t seem to get the COM Poisson model converging; not sure why!

That’s unfortunately not that surprising, it is tricky one…

I would also try the cratio ordinal model as one of the options - it can be interpreted as a process where you have separate probabilities for increase from 1 to 2, 2 to 3, etc. all depending on covariates (it requires positive response so you’d have to use trait + 1)

But otherwise I completely understand that you don’t want to invest in examining weird distributions.

Glad to have helped.

[jocateme]

So, I followed your suggestion and fitted the full model using the cratio family.

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

I’m astonished about how accurate the predictions of this model are. See for yourself:

Overall PP check

PP check by order

Plotted onto the phylogeny (inner ring = observed, outer ring = predicted)

These predictions are so accurate that I can’t help but feel that something must be wrong. I mean, is this too good to be true?!

The direction of the fixed effects are the same as in the Poisson model, which is also great:

### cratio

Group-Level Effects: 
~species (Number of levels: 762) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.79      1.05     1.39     5.48 1.00      624      773

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1]                 -0.40      5.09   -10.95    10.03 1.00     1068      617
Intercept[2]                 22.58      9.53     6.81    45.76 1.00      498      606
Intercept[3]                 44.86     17.45    20.50    89.78 1.00      543      702
continuous1                  -6.68      3.28   -14.86    -2.08 1.00      804     1030
continuous2                   7.56      3.07     3.31    15.16 1.00      836     1010
continuous1:continuous2       4.04      2.33     0.43     9.53 1.00     1244     1320

Family Specific Parameters: 
     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
disc     1.00      0.00     1.00     1.00 1.00     4000     4000


### poisson
Group-Level Effects: 
~species (Number of levels: 762) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.17      0.02     0.13     0.21 1.00     2383     2840

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                    -3.96      0.83    -5.71    -2.39 1.00     2299     2518
continuous1                  -0.71      0.23    -1.14    -0.26 1.00     3263     3043
continuous2                   0.76      0.19     0.39     1.14 1.00     3026     2594
continuous1:continuous2       0.51      0.19     0.14     0.90 1.00     3220     3325

The only thing I am wondering (and that I fear reviewers would wonder) is how valid it is to use an ordinal family when I know that a count process generated the data. But I guess that’s a discussion for another topic!

EDIT: Just found this super interesting paper. It shows that a probit ordinal regression model fits well to count data simulated under various assumptions, and that it performs well in identifying relevant covariates.

At the very least, I have another model backing up the inferences, which is fantastic. Once again, thanks for casting a light on this, Martin!

[martinmodrak]

@jocateme

These predictions are so accurate that I can’t help but feel that something must be wrong. I mean, is this too good to be true?!

To some extent, this is not very surprising: cratio has max_observed - 2 more parameters than Poisson (the additional intercepts), so it is just a much more flexible model.

@jocateme

The only thing I am wondering (and that I fear reviewers would wonder) is how valid it is to use an ordinal family when I know that a count process generated the data. But I guess that’s a discussion for another topic!

Actually cratio can be interpreted as a very flexible count model: the linear predictor (without intercept) can be seen as an overall tendency to develop spurs. The intercepts than describe how “hard” it is to move from no spurs to 1 spur, 1 spur to 2 spurs, etc. So to develop 3 spurs the species has to sequentially overcome 3 separate hurdles. You can imagine there are additional hurdles for higher spur counts up to infinity, although you can’t learn anything about them from your data (actually you could probably extract a lower bound on the intercept for “one more than max observed” from the fact that it was never overcome).

So the model is quite flexible but IMHO actually biologically interpretable :-)

You can make a similar argument for other ordinal families, but I find cratio to be most elegant for the case you describe. See https://psyarxiv.com/x8swp/ for details on the ordinal families brms supports.

[jocateme]

@martinmodrak

Actually cratio can be interpreted as a very flexible count model: the linear predictor (without intercept) can be seen as an overall tendency to develop spurs. The intercepts than describe how “hard” it is to move from no spurs to 1 spur, 1 spur to 2 spurs, etc. So to develop 3 spurs the species has to sequentially overcome 3 separate hurdles. You can imagine there are additional hurdles for higher spur counts up to infinity, although you can’t learn anything about them from your data (actually you could probably extract a lower bound on the intercept for “one more than max observed” from the fact that it was never overcome).

Loved this interpretation! Plus, I just realized there is a biological reason why using spur count as an ordinal variable makes sense in my case. What I’m really interested in is how energetically costly spurs are to bear in locomotion. This depends on the total mass of the spurs that an individual carries. But individual spurs have variable mass across species, meaning that, among species that possess, say, 2 spurs, the actual total spur mass is variable across species (while probably higher than species with 1 spur and lower than ones with 3 spurs). So, it makes total sense to think of total spur mass as a continuous latent variable, and spur count as the rank of total spur mass.

@martinmodrak

You can make a similar argument for other ordinal families, but I find cratio to be most elegant for the case you describe. See https://psyarxiv.com/x8swp/ for details on the ordinal families brms supports.

Fantastic recommendation; learned a lot reading this tutorial by @paul.buerkner ! In case anyone is interested, the peer-reviewed version can be found here.