[R] How to do poisson distribution test like this?

[(Ted Harding)]

On 29-Jul-09 03:28:24, glen_b wrote:
> Corrected links (the originals somehow aquired an extra space, sorry)
> 
> Paper: http://www.plantphysiol.org/cgi/content/full/148/3/1189
> Table I: http://www.plantphysiol.org/cgi/content/full/148/3/1189/TBL1

I think I've cracked it (using the informationm in Table 1).

What they seem to have done is:
1: Sum the oberved numbers N[i] of F-Box Genes: total = N
2: Sum the lengths L[i] of the chromosomes: total = L
3: Calculate expected number Exp[i] for chromosome i as
   Exp[i] = N*L[i]/L
4: Use the Poisson distribution as:
     ppois(N[i],Exp[i])    if N[i] < Exp[i]
     1-ppois(N[i],Exp[i]   if N[i] > Exp[i] (equality does not occur)

I have implemented this in R, using the Table 1 data, as follows.
The variables Exp0 and P0 are the values of Exp and P from the Table;
the variables Exp1 and P1 are as calculated according to the above.

Len  <- c(32.16,23.44,17.45,15.08,17.04,17.68,11.90,
          15.43,12.41,19.21,13.17,13.03,11.50,13.68,
          10.19,12.82,5.44,12.44,10.23)
Obs  <- c(36,25,13,10,19,19,12,15,12,15,17,13,13,13,12,9,2,12,2)
Exp0 <- c(30,22,17,14,16,17,11,15,12,18,13,12,11,13,10,12,5,12,10)
P0   <- c(0.137,0.235,0.235,0.159,0.196,0.242,0.340,
          0.391,0.395,0.273,0.082,0.353,0.207,0.421,
          0.176,0.231,0.112,0.398,0.004)

N <- sum(Obs) ; L <- sum(Len)
Exp1 <- N*Len/L
Low <- (Obs < Exp1)
P1 <- Low*ppois(Obs,Exp1) + (1-Low)*(1-ppois(Obs,Exp1))

cbind(Len,Obs,Exp0,Exp1,P0,P1)
         Len  Obs  Exp0       Exp1     P0           P1
#  [1,] 32.16  36    30  30.429265  0.137  0.136529344
#  [2,] 23.44  25    22  22.178544  0.235  0.234714001
#  [3,] 17.45  13    17  16.510904  0.235  0.234942347
#  [4,] 15.08  10    14  14.268449  0.159  0.158563376
#  [5,] 17.04  19    16  16.122969  0.196  0.196445135
#  [6,] 17.68  19    17  16.728526  0.242  0.241956811
#  [7,] 11.90  12    11  11.259585  0.340  0.340015730
#  [8,] 15.43  15    15  14.599613  0.391  0.390970369
#  [9,] 12.41  12    12  11.742139  0.395  0.394571188
# [10,] 19.21  15    18  18.176187  0.273  0.273013405
# [11,] 13.17  17    13  12.461238  0.082  0.082372362
# [12,] 13.03  13    12  12.328772  0.353  0.353596077
# [13,] 11.50  13    11  10.881112  0.207  0.207821286
# [14,] 13.68  13    13  12.943792  0.421  0.420776167
# [15,] 10.19  12    10   9.641611  0.176  0.175748117
# [16,] 12.82   9    12  12.130074  0.231  0.231213063
# [17,]  5.44   2     5   5.147239  0.112  0.112786229
# [18,] 12.44  12    12  11.770524  0.398  0.397809438
# [19,] 10.23   2    10   9.679458  0.004  0.003598524

There is now agreement between P0 (from the article) and P1
(as calculated above) -- subject to the rounded third decimal
place of P0 being 1 out in a few cases ([12], [13], [17]).

The puzzle clearly qrose in the first place becaue the authors
reported their Expected Numbers as integers, not fractions!

Well, well, well!
Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 29-Jul-09                                       Time: 11:59:11
------------------------------ XFMail ------------------------------