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]]
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
Locked
