Confidence intervals for the mean of a sample of counts

Question

I'm sure the answer to this question is obvious, but I can't seem to find it!

How do you calculate the confidence interval around the mean of a sample of counts?

I've made an example below in R, where I've bootstrapped confidence intervals. How do you calculate this using a formula? Is there a function in R?

# example sample of counts (n = 50, 'true' mean = 5)
set.seed(25); d <- rpois(50, 5)

# mean
mean(d)

# bootstrapped confidence intervals
B <- sample(d, 50 * 1000, replace = T)
B <- matrix(B, ncol = 1000)
B <- colMeans(B)
quantile(B, c(0.025, 0.975))

# mean = 4.72 (95% CI 4.14-5.34)

Answer 1

It's a bit nuanced.

You could pull out the big guns and use a poisson regression

# example sample of counts (n = 50, 'true' mean = 5)
set.seed(25)
d <- rpois(50, 5)

model = glm(d~1, family = poisson)

exp(confint(model))
>>>    2.5 %   97.5 % 
     4.143126 5.348048 

Or you could use the normal approximation since, as $\lambda$ gets big enough, the poisson starts to look a lot like a normal (though this only really happens when $\lambda$ is really really big). Be warned, if your sample is small and your $\lambda$ is close to zero, this may result in the CI covering values less than 0

mu = mean(d)
se = sd(d)/sqrt(length(d))

c(mu-1.96*se, mu+1.96*se)

>>>[1] 4.124606 5.315394

That the two estimates differ is perfectly acceptable since they make different assumptions. I would always prefer the former since the CI is computed using profile likelihood estimates.

Lastly, there is the procedure defined here though I think it adds little value over using glm.

Here are some simulations to explore the coverage properties for both methods given your data generating process:

set.seed(0)
simGlmCi<-function(){
  d <- rpois(50, 5)
  model = glm(d~1)
  ci = confint(model)
  covered = (5<ci[2])&(ci[1]<5)
  return(covered)
}

sims = purrr::rerun(1000,simGlmCi()) 
coverage = mean(unlist(sims))
coverage
>>>0.951


simNormalCi<-function(){
  d <- rpois(50, 5)
  mu = mean(d)
  se = sd(d)/sqrt(length(d))
  covered = (5<mu+1.96*se)&(mu-1.96*se)
  return(covered)
}

sims = purrr::rerun(1000,simNormalCi()) 
coverage = mean(unlist(sims))
coverage
>>>0.964

You might be better off with the profile likelihood approach since it seems the normal approach is a bit permissive.