R/poissonGoF.R

In helophilus/ColsTools: A variety of convenience tools and short-cuts

#' goodness-of-fit test for a Poisson distribution
#'
#' \code{poissonGoF(classes, freq, perm = 1000)}
#'
#' @param classes is an object giving the consecutive count classes, e.g., 0:5
#' @param freq are the frequencies of the classes in an object of the same length
#' as 'classes'.
#' @param is the number of randomizations to use, default is 1000.
#'
#' @return \code{p.value} is the probability dervived from the randomizations
#' @return \code{lambda} is the mean of the original data
#' @return \code{Chisq} is the calculated value of chi-squared from the actual data
#' @return \code{expected.values} are the expected values based upon a Poisson
#' distribution with the mean of the data for the given number of observations.
#'
#' @details The actual data are used to produce a Poisson distributed data set with
#' the same mean and number of values. This is compared to the actual data (a goodness
#' of fit test) by calculating a value of chi-squared. Randomized data sets are then
#' produced by sampling in the same proportions as dictated by the Poisson distribution.
#' A chi-squared value is calculated for the goodness-of-fit of each of these
#' 'randomized' distributions to the Poisson distribution, and the proportion of these
#' at least as big as the reference chi-squared value is the probability.

poissonGoF <-
  function (Classes, freq, perm = 1000)
  {
    options(warn = -1)
    classes = min(Classes):max(Classes)
    Mn = mean(rep(Classes, freq))
    props = Mn^classes/(factorial(classes) * exp(Mn))
    L = length(props)
    props[L] = 1 - sum(props[-L])
    EXP = props * sum(freq)
    S = sum(freq)
    TD = table(factor(rep(Classes, freq), levels = classes))[]
    chisqO = rep(NA, length(EXP))
    for (x in 1:length(EXP)) {
      chisqO[x] = (TD[x] - EXP[x])^2/EXP[x]
    }
    Chiref = round(sum(chisqO), 3)
    CHI2 = rep(NA, perm)
    for (i in 1:perm) {
      SP = sample(classes, prob = props, size = S, replace = T)
      TSP = table(factor(SP, levels = min(classes):max(SP)))[]
      LD = length(EXP) - length(TSP)
      if (LD > 0) {
        for (a in 1:LD) (TSP = c(TSP, 0))
      }
      chisq = rep(NA, length(EXP))
      for (b in 1:length(EXP)) {
        chisq[b] = (TSP[b] - EXP[b])^2/EXP[b]
      }
      CHI2[i] = sum(chisq)
    }
    options(warn = 0)
    PR = mean(CHI2 >= Chiref)
    cat("\n", "Poisson G-o-F Randomization test", "\n", "\n",
        "original chi-squared =", Chiref, "\n", "mean of data (lambda)=",
        Mn, "\n", "randomisation p-value =", PR, "\n", "\n",
        "number of random samples =", perm, "\n")
    cat("\n", "expected values: ", round(EXP, 2), "\n", "\n")
    OUT = list(p.value = PR, lambda = Mn, Chisq = Chiref, expected.values = EXP)
    return(invisible(OUT))
  }