# Comparing multiple incidence rates

I asked the following question on stack overflow yesterday: Negative binomial function in R

Reading comment 2, I understand that I cannot use a negative binomial modeling approach (a Poisson model works, but I suspect the assumption of equal mean and variance is invalid-though I'm uncertain as to how I can test this with an offset) and compare the betas to a reference category. I've googled and looked through my books but cannot find any other approach to compare multiple incidence rates.

b <- data.frame(
s=c(1800,539,490,301),
pop=c(2900000,1327000,880000,268000),
reg=c("A","B","C","D")
)

summary(pois.b<-glm(s~reg,offset=log(pop),data=b,family="poisson"))


So the question is : Is there any difference between the regions with regard to incidence?

Since the question yesterday was software related and today is more statistically flavored I figured it belonged here on cross-validated.

EDIT: Aug 11:

Since there are no other covariates here and the numbers are large I guess something as simple as

pairwise.prop.test(x=b$s,n=b$pop,p.adjust.method="bonferroni")


would get me a long way.

If you only have the four data points, I think the best way to do this is with a G^2 test. You want to start by assuming the frequency is a binomial distribution (every person in the population has the condition with probability p). And your null hypothesis is that p_1=p_2=p_3=p_4.

So the overall mean is (1800+539+490+301)/(2.9m+1.327m+.88m+.268m)=0.000582.

Your expected cases in each group are 1688.7, 772.7, 512.4, and 156.1. You can calculate the G^2 statistic, but the answer I get is 192.8, which is chi-squared(3) under the null hypothesis. This is a very low p-value, so you'd reject the null and say that yes, you can be quite confident that the incidence is different between these locations.

In particular that last location is considerably higher than the other three, so that is contributing heavily to the low p-value. You can repeat this analysis for the other three and you may get something a bit different, but that is an exercise to the reader :-)

HTH

ETA: the DF is 3, not 1, as Yves pointed out in the comments.

• This is a Likelihood-Ratio test. A LR test works also for the Poisson model where the expectation of $S$ is assumed to be $\lambda \times \text{Pop}$. The null hypothesis is then $\lambda_{\text{A}} = \lambda_{\text{B}} = \lambda_{\text{C}} = \lambda_{\text{D}}$. LR could also work for overdispersed the Negative Binomial case, but more data would be required to estimate the two parameters (using a non-canonical glm link). – Yves Aug 13 '13 at 10:00
• In this LR test, I think that the number of degrees of freedom to use is 3 = 4 - 1, since the resticted 'null' model has 1 df while the unrestricted has 4 df, one by region. – Yves Aug 13 '13 at 15:24
• Would it not be easier to just use - prop.test(x=b$s,n=b$pop) - and then use -pairwise.prop.test(x=b$s,n=b$pop,p.adjust.method="bonferroni")- to determine the individual differences? – Misha Aug 18 '13 at 11:50
• That ought to work too. That test I believe uses Pearson's Chi-Squared which is asymptotically equivalent to the G^2. If you want to test that the first three are equal, you'll have to just omit the fourth one to get an overall p-value. – Mike Nute Aug 19 '13 at 16:36

Given the limited data you have to work with you may only be able to address this question by incorporating additional assumptions (or data?) regarding the process behind these incidence rates, then doing some manual modeling. Any statistical technique you use will be implicitly making such assumptions for you under the hood, so better to call those out and structure your analysis around them.

You have observations of discrete incidence counts. The forms of some common discrete distributions encode the following assumptions:

• Poisson: variance is equal to the mean
• Binomial: variance is smaller than the mean
• Negative Binomial: variance is greater than the mean

You've already started down this road by ruling out a poisson model for the underlying process, saying variance = mean is not reasonable. If the process is a contagion model then it might be very reasonable to assume variance is greater than mean, so a negative binomial distribution.

The next question is fitting the parameters of the selected model and then making your comparisons. You could approach this in a couple of ways:

1. Empirically with your four data points - calculate mean and variance, then fit to distribution with old-fashioned algebra using distribution's mean and variance formulas. (You may need to standardize your data, for same reasons you'd use an offset in the glm.) Then calculate the probabilities of all 4 data points (and perhaps different combinations of 3) using the fitted model(s); lower probability suggests the process generating the incidence rates are not equivalent.
2. Use data from existing literature/research to fit the model; then test the probability of your incidence data occurring under that model. Poor fit for one of the data points could suggest that it's incidence rate deviates from the standard process in some way (or a poor model of course if most or all do not fit well).

Such results are hardly conclusive (nothing can be with only this data imho), but just as importantly it can inform dialogue and further research into the process you're modeling.