Bayesian GLM where the response variable is count classes

Question

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.

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.