R help - Feb 2006 - power and sample size for a GLM with poisson response variable

Hi all,

I would like to estimate power and necessary sample size for a GLM with 
a response variable that has a poisson distribution.  Do you have any 
suggestions for how I can do this in R?  Thank you for your help.

Sincerely,
Craig

-- 
Craig A. Faulhaber
Department of Forest, Range, and Wildlife Sciences
Utah State University
5230 Old Main Hill
Logan, UT 84322
(435)797-3892

Craig A Faulhaber wrote:
> Hi all,
> 
> I would like to estimate power and necessary sample size for a GLM with 
> a response variable that has a poisson distribution.  Do you have any 
> suggestions for how I can do this in R?  Thank you for your help.
> 
> Sincerely,
> Craig
> 
package asypow (on CRAN) or simulation.

Kjetil

Craig,   I found the package asypow difficult to use and it did not
yield results in the ballpark of other approaches.   (could not figure
out the "constraints" vector).   

I wrote some simple functions, one asy.pwr uses the non-central
chi-square distn.

asy.pwr<-function(counts=c(7,1),py=c(1081,3180))
{ # a two group Poisson regression power computation 
#  py is person-years or person-time however measured 
group<-gl(2,1)
ncp<-summary(glm(counts~group+offset(log(py)),family=poisson))$null.devi
ance
q.tile<-qchisq(.95,1)  #  actually just the X2 critical value of
3.841459
list(power=round(1-pchisq(q.tile,df=1,ncp),2))}

The second function, sim.pwr,  estimates power using simulated Poisson
random variates.  The most time consuming step is the call to rpois.
(Maybe someone knows a more efficient way to accomplish this?).  The
"for" loop is rather quick in comparison.   I hope you may find this
helpful, or if you have solved your problem some other way, please pass
along your approach.    Note, that for this problem,  very small values
of lambda, the two approaches give  much different power estimates  (96%
vs. 55% or so).   My problem may be better addressed as binomial
logistic regression, maybe then the simulation and the asymptotic
estimates my agree better.    

sim.pwr<-function(means=c(0.0065,0.0003),ptime=c(1081,3180),nsim=1000)
{ # a two group poisson regression power computation
# based simulating lots of Poisson r.v.'s 
# input rates followed by a vector of the corresponding person times
# the most time consuming part is the r.v. generation. 
# power is determined by counting the how often p-values are <= 0.05
group<-as.factor(rep(c(1,2),ptime))
rej<-vector(length=nsim)
y<-rpois(ptime[1]*nsim,means[1])
y<-c(y,rpois(ptime[2]*nsim,means[2]))
y<-matrix(y,nrow=nsim)
cat(sum(ptime)*nsim,"Poisson random variates
simulated","\n")
for(i in 1:nsim){res<-glm(y[i,]~group,family=poisson())
rej[i]<-summary(res)$coeff[2,4]<=0.05}
list(power=sum(rej)/nsim) }





	[[alternative HTML version deleted]]

Thanks for the functions, James.  I, too, did not understand the 
"constraints" vector in asypow.  Can anyone on the listserve explain
this?

Craig

Wassell, James T., Ph.D. wrote:
> Craig,   I found the package asypow difficult to use and it did not 
> yield results in the ballpark of other approaches.   (could not figure 
> out the "constraints" vector).  
>
> I wrote some simple functions, one asy.pwr uses the non-central 
> chi-square distn.
>
>
-- 
Craig A. Faulhaber
Department of Forest, Range, and Wildlife Sciences
Utah State University
5230 Old Main Hill
Logan, UT 84322
(435)797-3892