st: simulating random numbers from zero inflated negative binomial estimates

[Daniel Hill-McManus <D.Hill-McManus@sheffield.ac.uk>]

Hi,
I came across this code recently that Paul wrote in response to a query 
(see below).

I've also found it useful, thank you. But working through it I cannot 
understand why the linear predictor p2 is not exponentiated before being 
used in rgamma(1/alph, alph*p2).

I'd be grateful if someone could point out to me why this is.

Dan

On 3/06/2011 2:56 a.m., E. Paul Wileyto wrote:

    I'm not sure whether anyone has answered this yet.


    First, read the help on zinb post-estimation commands. There are 
many flavors of "predict" listed there. You will need three of them 
before you start generating random numbers. The first one you will need is:

    predict p1, pr


    That will generate a new variable, p1, which will be the predicted 
probability of an inflated zero. All the work is done for you.

    The second predicted quantity you will need is:

    predict lp, xb


    That will generate the linear combination of predictor variables 
weighted by coefficients for the negative binomial part of the model. 
Finally, you will need:

    predict alpha, xb eq(#3)


    which will generate a variable containing the overdispersion 
parameter for the negative binomial. With those three bits, you can get 
on to simulating.

    Here's my script:

    zinb cignums drug  week, inf(drug  week)
    predict p1 , pr
    predict p2 , xb
    predict lnalpha , xb eq(#3)
    gen alph=exp(lnalpha)
    gen xg=rgamma(1/alph, alph*p2)
    gen pg=rpoisson(xg)
    gen zi=runiform()>p1
    gen newcigs=zi*pg

    zinb newcigs drug  week, inf(drug  week)


    Paul


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/