[R] power and sample size for a GLM with Poisson response variable

From: Wassell, James T., Ph.D. <jtw2_at_CDC.GOV>
Date: Tue 28 Feb 2006 - 05:30:52 EST


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]]



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Feb 28 05:43:09 2006

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:42:46 EST