# power and sample size for a GLM with Poisson response variable

2 messages
Open this post in threaded view
|
Report Content as Inappropriate

## power and sample size for a GLM with Poisson response variable

 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