Bayes Negative Binomial (or use of $MH_p)?

[John Brehm]

I am struggling with learning the bayesmh commands, attempting to write my own negative binomial codes (for count) data. Negative binomials entail an ancillary parameter, and it isn't obvious to me how I make use of the built-in variable $MH_p, which appears to allow for it. So, two alternative questions:
- Has anyone worked up a bayesmh version of a negative binomial?
- Or, does anyone have bayesmh code which makes use of an ancillary parameter, through the $MH_p or otherwise?

Thanks!

[Houssein Assaad (StataCorp)]

Dear John,

The following program, nb_ll, is the negative binomial likelihood evaluator that you can use with bayesmh to fit a Bayesian negative binomial regression.

program define nb_ll
        version 14
        args lnf xb lnalpha
        tempname m
        tempvar p mu lnfj
        scalar `m' = 1/exp(`lnalpha')
        quietly {
                gen double `mu' = exp(`xb') if $MH_touse
                gen double `p' = 1/(1 + exp(`lnalpha') * `mu') if $MH_touse
                gen double `lnfj' = ln(nbinomialp(`m',$MH_y, `p')) if $MH_touse
        }
        summarize `lnfj' if $MH_touse, meanonly
        if r(N) < $MH_n {
               scalar `lnf' = .
               exit
        }
        scalar `lnf' = r(sum)
end

Note that I followed the same notation as in section "Methods and formulas" of the nbreg command. I have also used the built-in function nbinomialp() to simplify the code (see -help density_functions-). Now, let us put it to test:

webuse rod93
set seed 123
bayesmh deaths i.cohort, lleval(nb_ll, param({lnalpha}))  ///
        prior({deaths:},normal(0,100))                    ///
        prior({lnalpha},normal(0,100)) block({lnalpha}) dots

Model summary
------------------------------------------------------------------------------
Likelihood:
  deaths ~ nb_ll(xb_deaths,{lnalpha})

Priors:
  {deaths:i.cohort _cons} ~ normal(0,100)                                  (1)
                {lnalpha} ~ normal(0,100)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_deaths.

Bayesian regression                              MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
                                                 Number of obs    =         21
                                                 Acceptance rate  =      .3299
                                                 Efficiency:  min =     .07934
                                                              avg =     .09173
Log marginal likelihood = -123.24781                          max =      .1195
 
------------------------------------------------------------------------------
             |                                                Equal-tailed
             |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
deaths       |
      cohort |
  1960-1967  |  .0610405   .3371543   .011373   .0685907  -.6008041   .7465581
  1968-1976  | -.0414399   .3409068   .012103  -.0377467  -.7467722     .60975
             |
       _cons |   4.45491   .2396205   .008461   4.438088   4.010777   4.965601
-------------+----------------------------------------------------------------
     lnalpha | -1.032418   .3328226   .009629  -1.045511  -1.638257  -.3499949
------------------------------------------------------------------------------

The results are comparable to the results from nbreg:

nbreg deaths i.cohort

Negative binomial regression                    Number of obs     =         21
                                                LR chi2(2)        =       0.14
Dispersion     = mean                           Prob > chi2       =     0.9307
Log likelihood = -108.48841                     Pseudo R2         =     0.0007

------------------------------------------------------------------------------
      deaths |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      cohort |
  1960-1967  |   .0591305   .2978419     0.20   0.843    -.5246289      .64289
  1968-1976  |  -.0538792   .2981621    -0.18   0.857    -.6382662    .5305077
             |
       _cons |   4.435906   .2107213    21.05   0.000       4.0229    4.848912
-------------+----------------------------------------------------------------
    /lnalpha |  -1.207379   .3108622                     -1.816657   -.5980999
-------------+----------------------------------------------------------------
       alpha |     .29898   .0929416                      .1625683    .5498555
------------------------------------------------------------------------------
LR test of alpha=0: chibar2(01) = 434.62               Prob >= chibar2 = 0.000

[John Brehm]

Thanks! I had missed the built-in negative binomial density function entirely. Time to scour the list of density functions more closely!