Confidence interval for population mean when sample is a series of counts?

Question

I have count data for each of a sample of individuals (it's the number of times each independent individual performed a certain behaviour during a standardised observation of that individual). How can I calculate a confidence interval for the population mean? All over the internet there are Poisson methods for calculating the CI for a single count, but I can't find a method for a sample of separate counts. R code would be appreciated!

Edit: information about my data is useful. I have ten cells of data: Gender, by a five level ordinal variable, with small n in each cell. I am not certain the data is Poisson. I have done a possible dodgy Poisson regression which shows a significant effect of Var and an interaction with Gender. I would like to plot means with CIs to see if I believe the model.

Gender  Var  Counts
M       1    0 2 1 1 13 6 11 5 1
M       2    4 6 0 2 0 6 7
M       3    0 0
M       4    0 1
M       5    3 3 2 0 4 8
F       1    3 1 4 7 5 0 13 1
F       2    4 6
F       3    0 6
F       4    1 0 0
F       5    0 0 0 0 1 2

Answer 1

Here'a a quick display of the data:

If you're prepared to assume that within each Gender/Var combination the counts are Poisson with the same mean, then this is reasonably straightforward.

The variance of a sum of independent Poissons is Poisson, so you can form a confidence interval for the expected sum within each Var/Gender category using methods you say you're already familiar with. Since the mean is the sum divided by a constant, you simply divide the endpoints of your confidence interval for the expected sum by that constant. Done!

However, there's a fairly clear suggestion of more variation within groups than the Poisson model would indicate ("overdispersion").

Some possibilities - you could look at fitting an overdispersed Poisson model, or you could look at a negative binomial model, and then in either case obtain a confidence interval for each subgroup mean from that.

For example, here's the fitted means $\pm 1$ standard error of the mean for a main effects overdispersed Poisson model with log link (done in R):

I'd expect similar results from a negative binomial model. (You would possibly prefer a model with an interaction term though the results would look fairly similar to these.)

These standard errors would be about 1.8 times as big as you'd get by assuming a straight Poisson model.

Code:

a=data.frame(
        Gender=as.factor(c(rep("M",26),rep("F",21))),
        Var=rep(c(1:5,1:5),c(9,7,2,2,6,8,2,2,3,6)),
        Counts=c(0,2,1,1,13,6,11,5,1,4,6,0,2,0,6,7,0,0,0,1,3,3,2,
                0,4,8,3,1,4,7,5,0,13,1,4,6,0,6,1,0,0,0,0,0,0,1,2)
       )

VarGen=with(a,Var+(as.numeric(Gender)-1.5)/5)
plot(Counts~jitter(VarGen),pch=as.character(Gender), data=a,col=Gender)

VarGenfit=glm(Counts~as.factor(Var)+Gender,data=a,family=quasipoisson)
summary(VarGenfit)
aa=predict(VarGenfit,type="response",se.fit=TRUE)

plot(Counts~jitter(VarGen),pch=as.character(Gender), 
   data=a,col=c("grey","pink2")[a$Gender])
points(aa$fit~VarGen,pch="-",cex=1.8,col=a$Gender)
points((aa$fit+aa$se.fit)~VarGen,pch="-",cex=1.4,col=a$Gender)
points((aa$fit-aa$se.fit)~VarGen,pch="-",cex=1.4,col=a$Gender)
segments(VarGen,aa$fit-aa$se.fit,y1=aa$fit+aa$se.fit,col=a$Gender,lwd=2)