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]

From |
jhilbe@aol.com |

To |
statalist@hsphsun2.harvard.edu |

Subject |
Subject: st: simulating random numbers from zero inflated negative binomial estimates |

Date |
Tue, 15 Nov 2011 13:29:30 -0500 (EST) |

I hope this helps a bit. Joe Hilbe SYNTHETIC ZINB WITH LOGIT BINARY COMPONENT ===================================================

* Joseph Hilbe 15Aug2005; 10Oct2009; 5Jun2011 zinb_syn.do * LOGIT: x1= .9, x2= .1, _c= .2 * NB2 : x1=.75, n2=-1.25, _c=2, alpha=.5 clear set obs 50000 * set seed 1000 gen x1 = invnorm(runiform()) gen x2 = invnorm(runiform()) * NEGATIVE BINOMIAL- NB2 gen xb = 2 + 0.75*x1 - 1.25*x2 gen a = .5 gen ia = 1/a gen exb = exp(xb) gen xg = rgamma(ia, a) gen xbg = exb * xg gen nby = rpoisson(xbg) * BERNOULLI gen pi =1/(1+exp(-(.9*x1 + .1*x2+.2))) gen bernoulli = runiform()>pi gen zy = bernoulli*nby rename zy y * ZINB zinb y x1 x2, inf(x1 x2) nolog ==============================================

------------------------------------------------------------------------- -----

-------------+----------------------------------------------------------- ----- y |

-------------+----------------------------------------------------------- ----- inflate |

-------------+----------------------------------------------------------- -----

-------------+----------------------------------------------------------- -----

------------------------------------------------------------------------- ----- Date: Mon, 14 Nov 2011 16:04:17 +0000 From: 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

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

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

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/

- Prev by Date:
**RE: st: RE: Test of differences between geometric means** - Next by Date:
**RE: st: RE: Test of differences between geometric means** - Previous by thread:
**st: ML for logit/ologit** - Next by thread:
**st: xtpcse - "no time periods are common to all panels"** - Index(es):