How to perform a power analysis for poisson distributed data

Question

I want to perform a simplified power analysis to determine how many cells (samples) to have in control and experimental conditions (assume equal n across conditions). For an alpha level of 0.05, 0.95 power, and to detect a change of 10% in counts (average count in control is 0.2 to average count in experiment is 0.22 counts per sample). For t test or normal distribution, the pwr package in R and determining cohen's d is simple, but I am not sure how to perform this when working with counts data following a poisson distribution.

In case sd is needed, use 0.3.

Answer 1

The correct calculations depend a bit on which test you want to do to establish statistical significant. However, given that the effect size is relatively small, we can reasonably assume that your final $n$ per group will be quite large. Invoking the central limit theorem allows us to assume that the sampling distribution of the mean of samples from a Poisson distribution will tend to a normal distribution. To substantiate this claim, see the empirical sampling distribution of the mean of $n = 50$ draws from a Poisson distribution with $\lambda = 0.2$, as in your description, and using the fact that the mean of a Poisson distribution is equal to it's variance.

set.seed(123)
poismeans <- replicate(10000, rpois(50, 0.2) |> mean())
plot(density(poismeans))
curve(dnorm(x, 0.2, sqrt(0.2/50)), add = TRUE, lty=2)

^{Created on 2024-06-13 with reprex v2.0.2}

Accordingly, we can employ a two-sample t-test to evaluate the difference between the groups, and we can calculate the power of the test using the power.t.test() function from the stats library.

We define the effect size to be $\delta = 0.22 - 0.20 = 0.02$ and the pooled standard deviation $\sigma_p = \sqrt{(0.22^2 + 0.2^2)/2}$. Accordingly, fill in the power as $0.95$, which yields the following code and output.

m1 <- 0.22
m2 <- 0.2
power.t.test(power = 0.95, delta = m1 - m2, sd = sqrt((m1 + m2)/2))
#> 
#>      Two-sample t test power calculation 
#> 
#>               n = 13645.41
#>           delta = 0.02
#>              sd = 0.4582576
#>       sig.level = 0.05
#>           power = 0.95
#>     alternative = two.sided
#> 
#> NOTE: n is number in *each* group

^{Created on 2024-06-13 with reprex v2.0.2}

To verify the approach, let's fill in these numbers in a small simulation, to evaluate whether we obtain similar power empirically. Here, we draw two samples from the same size ($n = 13645$) from two Poisson distributions with different mean parameters $\lambda$.

f <- function(n, lambda1, lambda2) {
  x1 <- rpois(n, lambda1)
  x2 <- rpois(n, lambda2)
  t.test(x1, x2)$p.value
}

set.seed(123)
pvals <- replicate(10000, f(13645, 0.22, 0.2))
mean(pvals < 0.05)
#> [1] 0.948

^{Created on 2024-06-13 with reprex v2.0.2}

Finally, let's quickly evaluate whether we have not artificially inflated the Type-I error rate by using a t-test for Poisson data, assuming that there is no difference between the two groups (let them have the same mean parameter $\lambda = 0.2$).

f <- function(n, lambda1, lambda2) {
  x1 <- rpois(n, lambda1)
  x2 <- rpois(n, lambda2)
  t.test(x1, x2)$p.value
}

set.seed(123)
pvals <- replicate(10000, f(13645, 0.2, 0.2))
mean(pvals < 0.05)
#> [1] 0.0492

^{Created on 2024-06-13 with reprex v2.0.2}

These results do not indicate any problem, as the results yield the nominal Type-I error rate of 0.05 (approximately).