Power calculation by simulation for poisson regression

[Madu Abuchi]

Hello all,

I found this code for calculating power and sample size for a poisson regression model at https://www.stata.com/support/faqs/s...by-simulation/ but unfortunately the code comes up with error :

generate x=`d1' in 1
in not found
r(111);

Can any one help suggest what could have been the problem. The code is

/* poispow.do - does power estimation for Poisson regression model

model: log(mu) = b0 + b1*x
       y ~ Poisson(mu)

Specifically, power is estimated for testing the hypothesis b1 = 0, against a
user-specified alternative for a user-specified sample size. Without loss of
generality, we can assume the true value of b0 is 0.

In the `args' statement below:

        N is number of iterations in the simulation
        r is the number of rats receiving the jth dose (j=1,3)
        d1, d2 and d3 are the fixed dosages
        b1 is the "true" value of b1 (the alternative hypothesis).
*/

version 15  
args N r d1 d2 d3 b1
drop _all
set obs 3
generate x=`d1' in 1
replace x=`d2' in 2
replace x=`d3' in 3
expand `r'
sort x
generate mu=exp(0 + `b1'*x)   /* (the "zero" is to show b0 = 0) */
save tempx, replace

/* Note: Here, I generated my "x" and mu-values and stored them in a dataset
-tempx- so that the same values could be used throughout the simulation  */

/* set up a dataset with N observations which will eventually contain N
"p"-values obtained by testing b1=0.  */
drop _all
set obs `N'
generate pv = .

/* Loop through the simulations */
local i 0
while `i' < `N' {
     local i=`i'+1
     preserve

     use tempx,clear        /* get the n = 3*r observations of x
                               and the mean mu into  memory  */
     gen xp = rpoisson(mu)             /* generate n obs. of a Poisson(mu)random
                               variable in variable -xp- */
     quietly poisson xp x           /* do the Poisson regression */
     matrix V=e(V)          /* get the standard-error matrix */
     matrix b=e(b)          /* get the vector of estimated coefficients */
     scalar tv=b[1,1]/sqrt(V[1,1])          /* the "t"-ratio */
     scalar pv = 2*(1-normal(abs(tv)))    /* the p-value  */
     restore /* get back original dataset with N observations */
     quietly replace pv=scalar(pv) in `i'           /* set pv to the p-value for
                                               the ith simulation */
      _dots `i' 0
}

/*The dataset in memory now contains N simulated p-values. To get an
estimate of the power, say for alpha=.05, just calculate the proportion
of pv's that are less than 0.05:  */

count if pv<.05
scalar power=r(N)/`N'
scalar n=3*`r'
noisily display "Power is " %8.4f scalar(power) " for a sample size of " /*
    */ scalar(n) " and alpha = .05"
exit


******************************************************
/*We now run this program for various values of r (10, 20, 30, 40, 50, 60, and 70), w
    ith 100 simulations per run and with α = 0.05 hard coded in (of course this
    could easily be changed or made a user input)*/.

set seed 1234

run poispow 100 10 .2 .5 1.0 0.64
*Power is   0.3100 for a sample size of 30 and alpha = .05

run poispow 100 20 .2 .5 1.0 0.64
*Power is   0.5700 for a sample size of 60 and alpha = .05

run poispow 100 30 .2 .5 1.0 0.64
*Power is   0.7300 for a sample size of 90 and alpha = .05

[Maarten Buis]

The example is correct, the problem is the way you executed the example. On the FAQ you see that what you put in one .do file is actually two .do files. That is important. In essence the poispow.do file has become a program. What you need to do is store the .do file poispow.do and than you can from another .do from another do file that runs the simulations.

[Madu Abuchi]

Hi Maarten, thanks for your comment, but sorry, I don't seem to understand. Even when I save the first set of the .do file poispow.do, it still doesn't run. The error occurs at the beginning of the program where you have

version 15  
args N r d1 d2 d3 b1
drop _all
set obs 3
generate x=`d1' in 1

[Jean-Claude Arbaut]

Works for me.

You must store everything from

/* poispow.do - does power estimation for Poisson regression model

to

exit

in a file called poispow.do, and you must not execute that file or its contents right now.

Then, you run the lines that follow exit. They will call poispow.do with arguments.

The error comes from the fact that you ran the code directly, hence there are no arguments, hence the d1 local macro is empty, and the line generate x=`d1' in 1 is interpeted as generate x=in 1, and of course in does not exist (it's not a variable nor a scalar).


Note (read only if you understand the preceding, as it may be a bit confusing): as a matter of fact, it's harmless if you store the lines that run poispow.do in the same file (after exit), because the exit command ensures those lines will not be executed when the poispow.do file is called by run. You can just select those lines in Stata editor and run them, and this will also work.

Hope this helps,

Jean-Claude Arbaut

[Madu Abuchi]

Thanks Jean-Claude, very helpful. It works now
Cheers,
Madu