# Power simulation for Poisson Regression

I am interested in a power analysis for a Poisson regression of the form:

$Y_{i,1} = \beta_{0}+\beta_{1}Y_{i,0}+\beta_{2}Z_{i}$

where...

$Y_{i,1}$=The count from the 'post period' for case $i$

$Y_{i,0}$=The count from the 'pre period' for case $i$

$Z_{i}$=An indicator (0 or 1) of treatment for case $i$

There are $n_{0}$ and $n_{1}$ samples from each the control and treatment cases respectively.

I am wondering how to best simulate for a power analysis? What I know:

• We need to declare a minimally important effect size in terms of a multiplier for the treatment, the $\beta_{2}$. Lets say this is a six percent increase on average so effect is log(1.06)
• I do have data without a treatment and ran the model above with the treatment indicator and have the parameter estimates $\beta_{0}$ = 5.083739 and $\beta_{1}$ = 0.002056
• Since I have the values of $Y_{i,0}$ from the time without a treatment, I could use fitdist() in R for example to estimate the parameters of this data in say a poisson

Question: What uncertainty needs to be captured in a power analysis like this?

Is it sufficient to do the following (repeated nsim times) or is this not capturing enough uncertainty by using the static values of $\beta_{0}$ and $\beta_{1}$ ? If not how?

#loop this nsims times
n0<-30
n1<-30
B0<-5.083739
B1<-0.002056
effect<-1.06

#treatment indicator
treat<-c(rep(0,n0),rep(1,n1))

#sample of values for the "pre-period" based on prior data this is how the pre-period is distributed
Y0<-as.integer(rpois(n0+n1,505))

lin_pred<-B0+B1*Y0+log(effect)*treat
Y1<-as.integer(exp(lin_pred))

mod<-glm(Y1~Y0+treat,family=poisson)
#extract the p-value from here
summary(mod)


To Greg's questions, the data is paired in the sense that the same units (geographical units such as markets) are measured twice. $Y_{i,1}$ is the response variable after the intervention period (only some markets are treated of course) and $Y_{i,0}$ is the value before the intervention period....a "Pre-test" sometimes called.

This methodology comes from Google's Data Science team here and Andrew Gelman (Here)

a snippet of the former is: a snippet of this later reference is: Here is another take at this which I believe is closer to Greg's point about $Y_{1}$ not having any randomness except as inherited from $Y_{0}$

poisson.power.1cov<-function(alpha=0.05,nsims=NULL,n0=NULL,n1=NULL, covariate=NULL, B0=NULL, B1=NULL, possible.effects=NULL )
{
#hold power estimates for each effect
powers <- rep(NA, length(possible.effects))

#treatment indicator
treat<-c(rep(0,n0),rep(1,n1))

for (i in 1: length(possible.effects)) #loop over effects sizes
{
hold_power<-c()

for (j in 1:nsims)
{
#linear predictor using assumed B0, B1 and the simulated covariate Y0
lin_pred<-B0+B1*Y0+log(possible.effects[i])*treat

#exponentiate
mu<-exp(lin_pred)

#simulate Y1 based on mu (so here we are creating n0+n1 Y1 values, each of which has a different mu)
Y1<-rpois(n=length(mu),lambda = mu)

#fit model
mod<-glm(as.integer(Y1)~as.integer(Y0)+treat,family=poisson)

#extract p-value
sum_mod <-summary(mod)
hold_power<-rbind(hold_power,sum_mod$coefficients[3,4]) } #loop sims powers[i]<-mean(hold_power < alpha) } #loop effects return(powers) } #effects of interest effects<-seq(from=1.01, to=1.20, by=0.01) #determined from aprior data Y0<-as.integer(rpois(30,505)) #run power analysis with fixed B0 and B1 powers<-poisson.power.1cov(alpha = 0.05,n0 = 15,n1 = 15,nsims = 1000,covariate = Y0, B0 = 5.083739, B1 = 0.002056, possible.effects = effects) #plot plot(effects,powers) lines(effects,rep(0.8,length(effects)))  Which yields: ADD #3 Adds in simulation of Y0 at each sim, versus a static set poisson.power.1cov<-function(alpha=0.05,nsims=NULL,n0=NULL,n1=NULL, covariate_mu=NULL, B0=NULL, B1=NULL, possible.effects=NULL ) { #hold power estimates for each effect powers <- rep(NA, length(possible.effects)) #treatment indicator treat<-c(rep(0,n0),rep(1,n1)) for (i in 1: length(possible.effects)) #loop over effects sizes { hold_power<-c() for (j in 1:nsims) { #simulate Y0 Y0<-rpois((n0+n1),lambda = covariate_mu) #linear predictor using assumed B0, B1 and the simulated covariate Y0 lin_pred<-B0+B1*Y0+log(possible.effects[i])*treat #exponentiate mu<-exp(lin_pred) #simulate Y1 based on mu (so here we are creating n0+n1 Y1 values, each of which has a different mu) Y1<-rpois(n=length(mu),lambda = mu) #fit model mod<-glm(as.integer(Y1)~as.integer(Y0)+treat,family=poisson) #extract p-value sum_mod <-summary(mod) hold_power<-rbind(hold_power,sum_mod$coefficients[3,4])

} #loop sims
powers[i]<-mean(hold_power < alpha)

} #loop effects

return(powers)

}

#effects of interest
effects<-seq(from=1.01, to=1.20, by=0.01)
#determined from aprior data
covariate_mu <-505
#run power analysis with fixed B0 and B1
powers<-poisson.power.1cov(alpha = 0.05,n0 = 15,n1 = 15,nsims = 1000,covariate_mu =covariate_mu , B0 = 5.083739, B1 = 0.002056, possible.effects = effects)

#plot
plot(effects,powers)
lines(effects,rep(0.8,length(effects)))


Your $Y1$ variable does not have any randomness to it other than what it inherits from $Y0$.

Is this meant to be paired data? i.e. 2 measurements on each subject, one in each time period or with each treatment? If so, then you should probably be analyzing using a glmm model (lme4 package) instead of a glm model.

Also it is not completely clear what the relationship/difference is between your time periods and your treatment.

In general I would create the lin_pred variable based on your estimated coefficients, then generate $Y0$ and $Y1$ from a Poisson distribution with mean equal to exp(lin_pred).

• Greg, I added "Add #1" to the question, is this more clear? Jan 11 '17 at 19:03
• Does that change your recommendation? I am not sure I completely follow simulating both from exp(lin_pred) as it is based off of Y0. Jan 11 '17 at 23:58
• @B_Miner, yes that is much more along the lines that I was thinking. Though you may want to use log(Y0) for the lin_pred. And it is better to pre-specify the size of hold_power rather than expanding using rbind. Jan 13 '17 at 16:18
• @B_Miner, it is looking good. I would still suggest using log(Y0) so that they are more on the same scale. Jan 17 '17 at 17:17
• @B_Miner, your method is still a variant on the bootstrap, just not the usual sample with replacement version (your ecdf is not the traditional step shaped, but linear between the points). Still a good approach. Jan 18 '17 at 20:03