NB multilevel model - Finding the grand-mean posterior predictive distribution

Modeling
#1

I tag people I saw have answered similar topics before @bbbales2 , @paul.buerkner

Let’s suppose I have negative-binomial data of this kind, from three operators

image

I would like to find the grand-mean posterior predictive distribution.

Let’s start from the wrong approach I would like to improve

As for comparison, this is the wrong way to do it, ignoring the categories

brm( bright ~ 1 , data = input_dat, family = negbinomial())

Posterior predictive distribution

image

Let’s try to do things right

brm( bright ~ (1 | operator), data = input_dat,  family = negbinomial())

 Family: negbinomial 
  Links: mu = log; shape = identity 
Formula: bright ~ (1 | operator) 
   Data: input_dat (Number of observations: 900) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~operator (Number of levels: 3) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.18      0.59     0.45     2.74 1.01      675      820

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     3.89      0.54     2.82     4.99 1.01      767      976

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape    29.17      2.17    25.08    33.76 1.00     1586     1733

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

My question is how do I generate the negative binomial distribution of the grand-mean?

(pseudocode)

  1. rnbinom(exp(rnorm(Intercept, sd_Intercept )), shape =shape )
  2. rnbinom(exp(Intercept), shape =shape )

I would say (1) makes sense, but I get way to much overdispersion

image
image1239×736 13.1 KB

Thanks!

#2

What do you mean exactly with wanting to find the grand-mean posterior NB distribution?

The approach (1) that’s leads to more overdispersion of counts than observed of the data is a direct consequence of the assumption that the three operators are a sample (of N = 3) from a larger population. The estimated mean and SD of the operator-level means will therefore be very uncertain. Your approach (1) - if I understand correctly - simulates counts from that larger population, for which you have quite uncertain estimates (because N = 3).

1 Like
#3

Yes it is true

is my strategy after having thought about it. In reality, I often have more hierarchical groups than 3.

#4

You could also consider performing a mixed posterior predictive check. This means that you only resample the counts, given the group means from the posterior (i.e. not resampling the group means from the population mean).