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: RE: Power calculation for ZINB and Poisson Models


From   "MacLennan, Graeme" <g.maclennan@abdn.ac.uk>
To   "statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu>
Subject   st: RE: Power calculation for ZINB and Poisson Models
Date   Thu, 2 May 2013 21:00:27 +0100

Christina, a baby example. I dont have access to Stata at the moment, but this should still work

* quick power calculation for detecting
* half a standard deviation with a sample size
* of 50 in each group (two sample t-test)
sampsi 0 0.5, sd(1) n(50)


* replicate using a simple Monte Carlo simulation
program drop _all
program define ttestsim, rclass
       version 11.2
        syntax [, obs(integer 1)]
        drop _all
        set obs `obs'
        gen group = mod(_n, 2)
        gen outcome = rnormal()
        replace outcome = outcome + 0.5 if group == 1
        ttest outcome, by(group)
        return scalar p = r(p)
end

* create data set of 1000 simulated datasets returning p
* values from ttests
quietly simulate p = r(p), reps(1000): ttestsim, obs(100)
* the proportion of p values less than 0.05 is  power
gen power = p <= 0.05
tab power



You'll see that the simulation power is very similar to the out of the can number.
HTH,
Graeme,


________________________________________
From: owner-statalist@hsphsun2.harvard.edu [owner-statalist@hsphsun2.harvard.edu] On Behalf Of Chris Ansen [lakridstina@gmail.com]
Sent: Wednesday, May 01, 2013 9:37 AM
To: statalist@hsphsun2.harvard.edu
Subject: Re: st: RE: Power calculation for ZINB and Poisson Models

Thank you Graeme.

I can now do a simulation fitted to my data. However the link you
refer to does not explain how to convert the simulated model (or
several simulated models) into a power estimation.

Anybody who knows how to do that??? I would really appreciate your input.

Kind regards
Christina

On Tue, Apr 30, 2013 at 11:05 AM, MacLennan, Graeme
<g.maclennan@abdn.ac.uk> wrote:
> Christina, I will not comment on code below, but instead suggest this post which covered ZINB simulation.
>
> http://www.stata.com/statalist/archive/2011-11/msg00748.html
>
> HTH. Graeme.
>
> -----Original Message-----
> From: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Chris Ansen
> Sent: 30 April 2013 08:16
> To: statalist@hsphsun2.harvard.edu
> Subject: st: Power calculation for ZINB and Poisson Models
>
> Dear all
>
> I am trying to do a power calculation with simulated data. Originally I wanted to do one for a zero-inflated negative binomial model, testing three null hypothesis as suggested by Williamson et al. "Power Calculations for ZIP and ZINB Models", Journal of Data Science 2007.
> It is found easily by typing the title to a Google browser. Here three null-hypotheses are tested:
>
> Three null-hypotheses:
> H0: lamda1 equals 0                      nbreg-part of the model
> H0: beta1 equals 0                         logistic part of the model
> H0: [beta1, lamda1] equals [0,0]      joined logistic and nbreg
>
> However I started with simulation of the simpler Poisson inspired by:
> http://www.stata-journal.com/sjpdf.html?articlenum=st0010
> and
> http://www.stata.com/support/faqs/statistics/power-by-simulation/
>
> My problem is that even though I copy the two different approaches and execute them - both gets stuck by the command "generate x =`d1´ in 1 with the error message "in not found" r(111);
>
> The commands looks like this:
> args N r d1 d2 d3 b1
> drop _all
> set obs 3
> generate x =`d1´ in 1
> and continues exactly as explained in the link
>
> The other rather long approach (see link under Poisson with a loop within a loop) - does not work either! I am using STATA SE 11.2
>
> Both examples are copied from the two mentioned links, the only thing I have changed is the `var´, because they copy differently. I have also changed them for this mail not to be bounced from statalist.
>
> I have tried the long approach with version 7.0, 11.0 and 11.2 - so this is not the case.
>
> I will welcome any advises and solutions - also a more simple approach is very welcome.
>
> Thank you for your precious time.
>
> Kind regards
> Christina
>
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/faqs/resources/statalist-faq/
> *   http://www.ats.ucla.edu/stat/stata/
>
>
> The University of Aberdeen is a charity registered in Scotland, No SC013683.
>
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/faqs/resources/statalist-faq/
> *   http://www.ats.ucla.edu/stat/stata/

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


The University of Aberdeen is a charity registered in Scotland, No SC013683.

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


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