Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: simulating random numbers from zero inflated negative binomial estimates


From   Daniel Hill-McManus <D.Hill-McManus@sheffield.ac.uk>
To   statalist@hsphsun2.harvard.edu
Subject   st: simulating random numbers from zero inflated negative binomial estimates
Date   Mon, 14 Nov 2011 16:04:17 +0000

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/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index