From: Wassell, James T., Ph.D. <jtw2_at_CDC.GOV>

Date: Tue 28 Feb 2006 - 05:30:52 EST

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) }

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

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
*