# Special case of comparing poisson-distributed counts

I have count data (microbiology lab work) in which I have a serial 10-fold dilution of cells, and then grow the dilutions in 2 different conditions. For each dilution, I have 24 wells that could 'grow' or 'not grow' after some time (same amount of time for both conditions). The count data at the "extinction dilutions" would tell me how many cells were initially present in undiluted sample. The count data might look something like this:

Dilution 1:10 1:100 1:10^3 1:10^4 1:10^5 1:10^6
number of wells grown in condition 1 24 24 24 24 24 8
number of wells grown in condition 2 24 24 14 1 0 0

In this case, doing a statistical test is largely unnecessary to tell me the counts from conditions 1 and 2 are different due to something other than random chance (a biological phenomenon such as growth inhibition). In the case of my experiments, it tells me that only a small proportion of the cells are capable of growing in condition 2.

However, if the differences in counts are much closer, as in:

Dilution 1:10 1:100 1:10^3 1:10^4 1:10^5 1:10^6
number of wells grown in condition 1 24 24 24 24 24 8
number of wells grown in condition 2 24 24 24 24 19 4

...running some kind of statistical test becomes necessary to reject the null. I understand that counts of cells in this case fall from a poisson distribution. If I was comparing counts from the same dilution, I could use poisson.test() in R. However, I am not sure how to compare counts from across different dilutions as in the above case. How would I go about comparing 19/24 wells grown in 1:10^5 dilution in condition 2 and 8/24 wells grown in 1:10^6 dilution in condition 1?

Bonus points if you can help me understand how to combine counts from across dilutions to make use of all the observed counts (19 in 1:10^5 and 4 in 1:10^6 vs. 8 in 1:10^6). Extra bonus points if you can help me construct a confidence interval for the ratio of the cell counts from each condition.

• Sometimes it's convenient to think of this as having a Poisson distribution but it's really binomial.
– DWin
Aug 17, 2022 at 22:50
• "I understand that counts of cells in this case fall from a poisson distribution." -- I doubt this. It looks like each of 24 wells has some probability of producing a "1" / success (otherwise a "0"; your counts are then the sum of those). If the probability is constant within those categories (within the set of 24) and the events are independent, then you'd have a binomial count, not a Poisson. If the probabilities are not constant then it would be Poisson-binomial instead (assuming independence still pertains). Aug 18, 2022 at 22:47
• If there's some known model for the general form of the way probability should change with dilution (due to physics, chemistry, etc), that information should inform the model for the counts. Aug 18, 2022 at 22:52

In other types of microbiology assays, for example those evaluating the number of bacterial colonies found in a Petri dish, Poisson models are appropriate. The situation you describe, as noted in comments, is better handled with a binomial model, as each of the 24 wells is rated as either a success or a failure.

You can evaluate the probability of a "successful" well, $$p$$, as a function of dilution and treatment. I show that here for a simple model treating the dilution and treatment as categorical predictors, but the principle could be applied to incorporate dilution factors on a continuous scale or concentration-response type handling of treatment conditions.

In R, set up a data frame with the number of successes (grow) and failures (no-grow) and the corresponding dilutions and treatments. Omit the conditions in which all wells were successes:

dilDat
#   success fail  dil cond
# 1      24    0 five  one
# 2      19    5 five  two
# 3       8   16  six  one
# 4       4   20  six  two


A binomial logistic regression estimates the log-odds ($$\log[p/(1-p)]$$) of success as a function of dilution (dil) and condition (cond):

summary(glm(cbind(success,fail)~ dil + cond,
data=dilDat, family=binomial))

# (some lines omitted)
#
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept)   3.0862     0.6810   4.532 5.84e-06
# dilsix       -3.5884     0.6591  -5.445 5.19e-08
# condtwo      -1.4622     0.6271  -2.332   0.0197
#
# (other lines omitted)


I omitted some further details from the model summary that are useful for evaluating the quality of the model but not needed to make the basic point here.

The (Intercept) is the estimated log-odds of success in the reference situation: dilution of 1:10^5 ("five" in the data frame) and treatment condition "one." The coefficient for dilsix is the difference in log-odds between that situation and a dilution of 1:10^6 ("six" in the data frame) under treatment condition "one"; the coefficient for condtwo is the difference in log-odds from the reference situation for treatment condition "two."

This model assumes that the associations of dil and cold with outcome are additive on the log-odds scale. The model gets the correct success numbers within about $$\pm 1$$:

 24*predict(glm(cbind(success,fail)~dil+cond,data=dilDat,family=binomial),type="response")
#         1         2         3         4
# 22.951682 20.048318  9.048318  2.951682


You can use the standard errors of the coefficients to get confidence intervals for any differences in log-odds, which can then be translated to probabilities. The 95% CI for each coefficient in the table (difference from reference condition) is 1.96 times the corresponding SE value around the point estimate in the log-odds scale.

For comparisons other than against the reference condition, you need to take the coefficient covariance matrix (not shown, but present in the model) into account. Many find the emmeans package helpful for such calculations, as it handles many types of models and can transform results from log-odds into probabilities.