# G*Power and simulation produce very different sample size estimates for a Poisson regression

I've been attempting to conduct a sample size calculation for a Poisson regression. G*Power produced a sample size of 472.

Parameters for G*Power

Tails = 2

Exp(B1) = 1.233

alpha = 0.05

Power = 0.8

Base Rate Exp(B0) = 1.37

Mean exposure = 1

R2 from other variables = 0

X Distribution = Binomial

X parm π = 0.5

However, when I try to simulate a sample size calculation in R, a sample size of 472 tends to result in about 25% power

N_per_grp = 236
significant_count = 0

for (i in 1:10000) {

grp_1 = data.frame(rep(1,N_per_grp), rpois(N_per_grp, 1.233))
names(grp_1) = c("grp", "outcome")

grp_0 = data.frame(rep(0,N_per_grp), rpois(N_per_grp, 1.37))
names(grp_0) = c("grp", "outcome")

df = rbind(grp_1, grp_0)

fit = summary(glm(outcome ~ grp, data = df, family = poisson(link = "log")))

p = fit$coefficients[2,4] if (p <= .05) { significant_count = significant_count + 1 } }  Does anyone know why this is? • Shouldn't the rate in the treatment group be 1.233 * 1.37? The way you're simulating it now, there's a base rate of 1.37, and a RR of about 1.233/1.37 = 0.9... – if_the_correlations_are_zero Jan 16 '19 at 18:20 • Ah, I see what you're saying. I was interpreting the Exp(B1) in G*Power as the rate in the intervention group – TPM Jan 16 '19 at 18:22 • It's the rate, relative to the control group. When you change that do you get consistent answers? – if_the_correlations_are_zero Jan 16 '19 at 18:23 ## 1 Answer You're not programming the relative rate of 1.233 correctly. The rate in group 1 needs to 1.233 times larger than the rate in group 0...Not just 1.233 set.seed(1) N_per_grp = 236 significant_count = 0 for (i in 1:10000) { grp_1 = data.frame(rep(1,N_per_grp), rpois(N_per_grp, 1.233*1.37)) names(grp_1) = c("grp", "outcome") grp_0 = data.frame(rep(0,N_per_grp), rpois(N_per_grp, 1.37)) names(grp_0) = c("grp", "outcome") df = rbind(grp_1, grp_0) fit = summary(glm(outcome ~ grp, data = df, family = poisson(link = "log"))) p = fit$coefficients[2,4]

if (p <= .05) {
significant_count = significant_count + 1
}
}
significant_count/10000
[1] 0.7937


Looks about like 80% power to me....