Description of data

I have to analyze some data where the response variable is the counts of number of insects observed feeding on a bait at many time points. The treatments are three different types of bait and there are several replicates per treatment. However it was not possible to count the number of insects exactly at each time point. Instead, each time point was assigned a count class. The designations and ranges of values for each class are as follows:

  • 0
  • 1
  • 2
  • 5 (3-6)
  • 10 (7-14)
  • 20 (15-34)
  • 50 (35-74)
  • 100 (75-154)
  • 200 (155-249)
  • 300 (>249)

As you can see the number used to designate the class is roughly but not always exactly the midpoint of the range.

Model fitting problem

I want to fit a model to help determine whether the abundance of insects differs between treatments.

My first thought was to just treat the class designations as if they were true counts and use a zero-inflated poisson response distribution. That seems wrong because it ignores the fact that a value of 200 may not actually represent 200 insects; the true underlying value could be anywhere from 155-249.

I also was considering fitting an ordered multinomial and treating the count classes as classes. That seems more appropriate but it also bothers me, because it is throwing away information. The count classes do convey at least some information about the relative number of ants but if you convert it to ordinal classes, you lose all that information except for class 1 < class 2 < class 3, and so on.

Snippets of data and code

The first few rows of the data look like this, after converting the abundance class to an ordered factor using the ordered() function in R.

Rep    Trtmt Time_posttrt Abund_class Abund_class_ordinal
<dbl> <chr>        <dbl>       <dbl> <ord>              
 1     2 Ctrl             5          20 20                 
 2     2 Ctrl            10          10 10                 
 3     2 Ctrl            15          50 50                 
 4     2 Ctrl            20         100 100                
 5     2 Ctrl            25         100 100                
 6     2 Ctrl            30          50 50                 
 7     2 Ctrl            35         100 100                
 8     2 Ctrl            40         100 100                
 9     2 Ctrl            45         100 100                
10     2 Ctrl            50         100 100                

My brms formulas in R that I have tried look like this (not including priors etc.):

  • Zero-inflated Poisson: brm(Abund_class ~ Trtmt + (1 + Time_posttrt | Rep), family = zero_inflated_poisson)

  • Ordered multinomial: brm(Abund_class_ordinal ~ Trtmt + (1 + Time_posttrt | Rep), family = cumulative(link = 'logit'))

I would appreciate any general advice on fitting a model to this type of data. I can provide reproducible data and code if it would be helpful, but I was hoping for more general advice on model fitting with this type of data.

  • 1
    $\begingroup$ A first idea that comes to mind is to try to model this using a latent variable model. The full model would be $(Z_i, Y_i)$, where $Z_i$ is the number of insects (latent) and $Y_i$ is an indicator of the range interval as you described (observed). If you assume $Z_i$ follows some count distribution data such as the Poisson or NegBin, maybe you could use EM, or its bayesian version, to estimate the parameters. However, this is only "a first idea" since I do not know which software you could use to fit such a model, and maybe you would have to work out the details yourself. $\endgroup$ Feb 16, 2022 at 21:26
  • 1
    $\begingroup$ Did you try an ordinal regression model? statsmodels.org/dev/examples/notebooks/generated/… $\endgroup$ Feb 17, 2022 at 3:32
  • 1
    $\begingroup$ @RajivSambasivan I'm working on ordinal regression at the moment for this problem, it seems to be working well so far. $\endgroup$
    – qdread
    Feb 17, 2022 at 14:44

1 Answer 1


I don't think this is going to be super easy, but none the less let's try.

Let $d(x;{\theta})$ be the distribution of counts for a group (txt or control) which is parameterized by $\theta$ (so if the distribution is Poisson, then $\theta$ would be $\lambda$ for example). Let $D(x; \theta)$ be the cumulative distribution of $d$.

Much of your data is of the form $D(x+z;\theta) - D(x;\theta)$. For example, the proportion of your sample in bucket 5 the treatment group is an estimate of $D(6; \theta) - D(3;\theta)$.

You can consider your data in txt/control $\mathbf{y}$ as multinomial with multinomial parameter $\psi$

$$ y \sim \mbox{Multinomial}(\psi, N)$$

Where the $i^{th}$ element of $\psi$ can be expressed as the difference in CDF intervals (as I've shown above).

To my knowledge, there is no brms equivalent for this model. I did however write a similar model here (Stan code included) so if you are so inclined you could probably change my Stan code to fit your case.

  • $\begingroup$ This post and the model you link to are extremely helpful. The data you are fitting your model to are very similar to what I'm working with. I hope to code up an extension of your model that includes a time series component. Unfortunately that's a big job and I'm not sure if I will have time to do it in the near future. I will report back here if and when I get the chance to write it!!! $\endgroup$
    – qdread
    Feb 17, 2022 at 17:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.