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, 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 <>
Subject   st: simulating random numbers from zero inflated negative binomial estimates
Date   Mon, 14 Nov 2011 16:04:17 +0000

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.


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)


*   For searches and help try:

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