Bayesian random parameter poisson regression using bayesmh lleval option

[Behram Wali]

Dear all,

I am trying to investigate an empirical problem that requires estimation of Bayesian Poisson regression with random parameters.

I am aware that Poisson distribution is one of the many built-in likelihoods that comes built-in Stata's bayesmh command.

As such, I tried the following:

sysuse auto                                            //Example 
gen id = _n                                            //id created just for the sake of demonstration
//Bayesian Random Coefficients Poisson using bayesmh "likelihood" option
set seed 123
    bayesmh trunk i.id#c.weight i.id#c.price, likelihood(poisson)  ///
        prior({trunk:i.id#c.weight},normal({trunk:weight},{Sigma1})) 
        prior({trunk:i.id#c.price},normal({trunk:price},{Sigma2}))    
        prior({trunk: _cons weight price}, normal (0,100))                                    ///
        prior({Sigma1}, igamma(0.001, 0.001))                                                ///
        prior({Sigma2}, igamma(0.001, 0.001))                                                ///                                                         
         block({Sigma1}, gibbs)                                                        ///
        block({Sigma2}, gibbs)                                                        ///
        block({trunk:i.id#c.weight})                                                        ///
        block({trunk:i.id#c.price})                                                         ///                                                    ///
        exclude({trunk:i.id#c.weight}{trunk:i.id#c.price}) dots                                         
        bayesstats ic

My understanding is that bayesmh likelihood option can be slow for complex problems. Am i right? Perhaps, that's why Stata has taken almost two hours before even starting the burn-in iterations.

With this said, i thought that using bayesmh "llevaluator" option can be faster. Below is the user-defined Poisson likelihood which is then called by bayesmh llevaluator option.

//User Defined Poisson Likelihood
program drop wal_poiss
program define wal_poiss
  version 14.2
  args todo b lnf        
  tempvar theta1         
  mleval `theta1' = `b', eq(1)   
  local y "$ML_y1"    
  tempvar lnyfact mu
  quietly gen double `lnyfact' = lnfactorial(`y')
  quietly gen double `mu' = exp(`theta1')
  mlsum `lnf' = -`mu' + `y'*ln(`mu') - `lnyfact'
end
ml check

//ml model trunk weight price
//ml maximize

//Bayesian Random Coefficients Poisson using bayesmh "lleval" option                                                                                
set seed 123
    bayesmh trunk i.id#c.weight i.id#c.price, lleval(wal_poiss)  ///         I suspect am missing important attribute of the log-likelihood evaluator in the lleval option
        prior({trunk:i.id#c.weight},normal({trunk:weight},{Sigma1}))                         ///
        prior({trunk:i.id#c.price},normal({trunk:price},{Sigma2}))                             ///   
        prior({trunk: _cons weight price}, normal (0,100))                                    ///
        prior({Sigma1}, igamma(0.001, 0.001))                                                ///
        prior({Sigma2}, igamma(0.001, 0.001))                                                ///                                                         
         block({Sigma1}, gibbs)                                                                ///
        block({Sigma2}, gibbs)                                                                ///
        block({trunk:i.id#c.weight})                                                        ///
        block({trunk:i.id#c.price})                                                         ///                                                    ///
        exclude({trunk:i.id#c.weight}{trunk:i.id#c.price}) dots                                         
        bayesstats ic

When i run this, Stata returns the following error:

"ml model not found" and "wal_poiss could not be evaluated."

I am certain that i am missing an important attribute of the log-likelihood evaluator in the lleval option in above code.

Can someone please guide me in this regard?

Any help is highly appreciated!
-Behram

[Nikolay Balov (StataCorp)]

Hi Behram,

The maximum likelihood evaluators used by ml model are not compatible with bayesmh's likelihood evaluators. Even if you could make it work however, you would not see any speed up in simulation, due to the inefficient use of factor variable specification of random-effects in your bayesmh model.

I would suggest a much more efficient specification of the same model using the redefine option

set seed 123
bayesmh trunk = ({trunk__cons} + {weight:}*weight + {price:}*price), ///
    likelihood(poisson)                                              ///
    redefine(weight:i.id) redefine(price:i.id)                       ///
    prior({weight:i.id}, normal({trunk:weight},{Sigma1}))            ///
    prior({price:i.id},  normal({trunk:price}, {Sigma2}))            ///
    prior({trunk: _cons weight price}, normal (0,100))               ///
    prior({Sigma1}, igamma(0.001, 0.001))                            ///
    prior({Sigma2}, igamma(0.001, 0.001))                            ///                                                         
    block({Sigma1}, gibbs)                                           ///
    block({Sigma2}, gibbs)                                           ///
    exclude({weight:i.id}{price:i.id}}) burnin(5000) dots

Now, the simulation efficiency of this model is not high because of model misfit, but applied to appropriate data should work fine.

Nikolay