0
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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... $\endgroup$ – if_the_correlations_are_zero Jan 16 at 18:20
  • $\begingroup$ Ah, I see what you're saying. I was interpreting the Exp(B1) in G*Power as the rate in the intervention group $\endgroup$ – TPM Jan 16 at 18:22
  • 1
    $\begingroup$ It's the rate, relative to the control group. When you change that do you get consistent answers? $\endgroup$ – if_the_correlations_are_zero Jan 16 at 18:23
1
$\begingroup$

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

$\endgroup$

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.