Counts of events, based on the Poisson distribution, is a frequently encountered model in medical research. Examples of this are number of falls, asthma attacks, number of cells, and so on. The Poisson parameter Lambda (λ) is the total number of events (k) divided by the number of observation units (n) in the data (λ = k/n). The unit forms the basis or denominator for calculation of the average, and need not be individual cases or research subjects. For example, the number of asthma attacks may be based on the number of child months, or the number of pregnancies based on the number of women years in using a particular contraceptive.

Poisson is different to the binomial parameter of proportion or risk where proportion is the number of individuals classified as positive (p) divided by the total number of individuals in the data (r = p/n). Proportion or risk must always be a number between 0 and 1, while λ may be any positive number.

For examples, if we have 100 people, and only 90 of them go shopping in a week then the binomial risk of shopping is 90/100 = 0.9. However, some of the people will go shopping more than once in the week, and the total number of shopping trips between the 100 people may be 160, and the Poisson Lambda is 160/100 = 1.6 per person week

Large Lambda (λ=k/n) values, say over 200, assumes an approximately normal or geometric distribution, and the count (or sqrt(count)) can be used as a Parametric measurement. If the events occur very few times per individual, so that individuals can be classified as positive or negative cases, then the binomial distribution can be assumed and statistics related to proportions used. In between, or when events are infrequent, the Poisson distribution is used.

Some clarification of nomenclature may be useful.

• Counts of events (e.g. number of asthma attacks recorded) are represented by k. This count must be in terms of how many events over a defined period or environment (e.g. in 100 attacks in 300 children over 6 months, or 10 cells seen in 5 microlitres of fluid),
• The mean count, or count rate (k/n) is represented by λ. e.g. 100 attacks(k) in 1800 children months (n) produces λ=100/1800 = 0.06 attacks per child month (λ)

• The most commonly used method, initially described by Przyborowski and Wilenski (see reference), is known as the Conditional Test (the C Test). The test is based on the null hypothesis that the ratio of the two count rates (λ₂ / λ₁) is equal to 1.
• More recently, Krishnamoorthy and Thomson (see reference) proposed an improvement on the C Test, where the null hypothesis is that the difference between the two count rates (λ₂ - λ₁) is equal to 0. Althought computation for this test is more complex, the advantages are that it is more robust, and the results have greater power
• Whitehead (see reference), in his text book on unpaired sequential analysis, provided algorithms to determine sample sizes for non-sequential methods, and a method for comparing two counts at the end of the sequence. This test depends on a transformation of the difference into a normally distributed mean and compares the means against the null hypothesis of 0.
• Comparing the 3 tests: In most cases, the results from the 3 tests on the same set of data are approximately the same. There is only a need to choose when the ample size is small or when the difference between the two groups are minor, so that the statistical significance is ambiguous.
  • The advantage of the Whitehead algorithm is the speed of computation, as there is no need for repeated estimation of factorial numbers that are necessary in the other two tests. Users may prefer this test if the sample sizes are large or there are numerous calculations to be done
  • The C Test has been used the longest, and quoted by most text books.
  • The E Test is the most robust, least likely to result in a false statistical significance, and so is the safest one to use. With large numbers however, it is also the most computational intensive and time consuming
  In most cases, therefore, the C test should be used. If the sample size is too large and the calculations take too long, then the Whitehead algorithm can be used. When the statistical significance of the results is ambiguous (e.g. close to p<0.05) then the E Test is preferred.

Data Input

The columns are

• Col 1: k1, the total count in group 1
• Col 2: n1, the sample size in group 1
• Col 3: k2, the total count in group 2
• Col 4: n2, the sample size in group 2

More Complex Models

References

Poisson Distribution by Wikipedia

Steel RGD, Torrie JH Dickey DA (1997) Principles and Procedures of Statistics. A Biomedical Approach. The McGraw-Hill Companies, Inc New York. p. 558

Przyborowski J and Wilenski H (1940) Homogeneity of results in testing samples from Poisson series. Biometrika 31:313-323.

Krishnamoorthy, K and Thomson, J. (2004). A more powerful test for comparing two Poisson means. Journal of Statistical Planning and Inference, 119, 249-267. Program adapted from FORTRAN program by same author, downloaded from https://userweb.ucs.louisiana.edu/~kxk4695/statcalc/pois2pval.for

Whitehead John (1992). The Design and Analysis of Sequential Clinical Trials (Revised 2nd. Edition) . John Wiley & Sons Ltd., Chichester, ISBN 0 47197550 8. p. 48-50

• Subroutines that creates and access a vector of log(factorial) numbers, as these are repeatedly used in C and E Tests
• The input data and setting up of log(factorial) vectors as these are common to all 3 algorithms
• Three (3) algorithms. The whitehead algorithm, the C Test, and the E Test
• Results

# log(factorial) vector for repeated calculations of factorial and binomial coefficients
logFactVector = vector()

MakeLogFactVector <- function(n) # make log(factorial vector
{
  x = 0
  for(i in 1:n)
  {
    x = x + log(i)
    logFactVector <<- append(logFactVector,x)
  }
}

LogFact <- function(i)
{
  if(i==0) i=1
  return (logFactVector[i])
}

# binomial coefficient
BinomCoeff <- function (n,k)
{
  if(k==0 | k==n) return (1)
  return (exp(LogFact(n) - LogFact(k) - LogFact(n-k)))
}

# input data
myDat = ("
k1   n1   k2   n2
13   10   8   10
10   30   10  50
12   100  4    110
         ") 
myDat
df <- read.table(textConnection(myDat),header=TRUE)  # conversion to data frame    
df      # display input data (count(k) and sample size (n) for groups 1 and 2)

# Create log(factorial) vector
maxN = max(df)
MakeLogFactVector(maxN)

#   Pgm 1 Whitehead's algorithm
Pw <- vector()     # probability of Type I Error (Whitehead) one tail
for(i in 1:nrow(df))
{
  k1 = df$k1[i]
  n1 = df$n1[i]
  k2 = df$k2[i]
  n2 = df$n2[i]
  n = n1 + n2;               
  zi = (n1*k2 - n2*k1)/n
  vi = n1*n2*(k1+k2)/(n*n);
  z = abs(zi / sqrt(vi))
  Pw <- append(Pw,1 - pnorm(z))   # Whitehead alpha one tail
}
df$Pw <- Pw   # add Whitehead's algorithm results to data frame 

# Program 2: C Test
Pc <- vector()     # probability of Type I Error (C Test) one tail
for(i in 1:nrow(df))
{
  k1 = df$k1[i]
  n1 = df$n1[i]
  k2 = df$k2[i]
  n2 = df$n2[i]
  x = 1         # for ratio = 1
  if(k1/n1!=k2/n2)
  {
    x1=0
    x2=0
    k = k1 + k2
    p = 0;
    if(k1>k2)
    {
      p = (n1/n2)/(1.0 + (n1/n2))
    }
    else 
    {
      k1 = k2
      p = (n2/n1) / (1.0 + (n2/n1))
    }
    for(j in 0:1)
    {
      x = 0;
      if(j==0)
      {
        for(ii in k1:k) x = x + BinomCoeff(k,ii) * p^ii * (1-p)^(k-ii) # choose is binomial coefficient
      }
      else
      {
        for(ii in 0:k) x = x + BinomCoeff(k,ii) * p^ii * (1-p)^(k-ii)
      }
      if(j==0) { x1 = x  } else { x2 = x }
    }
    if(x1<x2) { x = x1; } else { x = x2; }
    if(x>1.0)x = 1.0;
  }
  Pc <- append(Pc,x / 2)
}
df$Pc <- Pc   # add C Test results to data frame 

#Program 3: E Test
#   Result vectors
Pe <- vector()     # probability of Type I Error (E Test) one tail
for(i in 1:nrow(df))
{
  k1 = df$k1[i]
  n1 = df$n1[i]
  k2 = df$k2[i]
  n2 = df$n2[i]
  diff = 0
  x = 1
  if(k1 / n1 != k2 / n2)
  {
    Vk = k1 / (n1*n1) + k2 / (n2*n2)
    Tk1k2 = abs((k1 / n1 - k2 / n2 - diff) / sqrt(Vk))
    Lambda2k = (k1+k2) / (n1+n2) - diff * n1 / (n1 + n2)
    SumT = 0
    SumF = 1
    x1 = 0
    while (x1<1000 && SumF>1e-10)
    {
      x2 = 0
      SumF = 0
      f = 1;
      while (x2<1000 & f>1e-10)
      {
        Vx = x1 / (n1*n1) + x2 / (n2*n2)
        if(abs(x1 / n1 - x2 / n2)<=diff) { Tx1x2 = 0 } 
        else { Tx1x2 = abs(x1 / n1 - x2 / n2 - diff) / sqrt(Vx); }
        if (Tx1x2>=Tk1k2)
        {
          v1 = exp(-n1 * (Lambda2k + diff))
          v2 = (n1 * (Lambda2k + diff))^x1
          v3 = exp(-n2 * Lambda2k) * (n2 * Lambda2k)^x2
          f = exp(log((v1*v2*v3)) - (LogFact(round(x1)) + LogFact(round(x2))))
          SumF = SumF + f
        }
        x2 = x2 + 1;
      }
      SumT = SumT + SumF
      x1 = x1 + 1
    }
    x = SumT / 2
  }
  Pe <- append(Pe,x)
}
df$Pe <- Pe   # add E Test results to data frame 

df # display input data and the results (Type I Error, one tail), Pw=Whitehead, Pc=C Test, Pe= E Test

> df      # display input data (count(k) and sample size (n) for groups 1 and 2)
  k1  n1 k2  n2
1 13  10  8  10
2 20  30 10  50
3 12 100  4 110
> 
  k1  n1 k2  n2         Pw         Pc         Pe
1 13  10  8  10 0.13761676 0.09582758 0.14336618
2 10  20 10  50 0.01694743 0.49414271 0.04229829
3 12 100  4 110 0.01415499 0.01249091 0.01611372
>