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