2
$\begingroup$

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) 

ADD #1

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:

enter image description here

a snippet of this later reference is:

enter image description here

ADD #2

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:

enter image description here

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)))
$\endgroup$
1
$\begingroup$

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

$\endgroup$
13
  • $\begingroup$ Greg, I added "Add #1" to the question, is this more clear? $\endgroup$
    – B_Miner
    Jan 11 '17 at 19:03
  • $\begingroup$ 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. $\endgroup$
    – B_Miner
    Jan 11 '17 at 23:58
  • 1
    $\begingroup$ @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. $\endgroup$
    – Greg Snow
    Jan 13 '17 at 16:18
  • 1
    $\begingroup$ @B_Miner, it is looking good. I would still suggest using log(Y0) so that they are more on the same scale. $\endgroup$
    – Greg Snow
    Jan 17 '17 at 17:17
  • 1
    $\begingroup$ @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. $\endgroup$
    – Greg Snow
    Jan 18 '17 at 20:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.