What is the best way to perform a goodness-of-fit test of data to a continuous distribution with the chi-square method in R?

Question

I'm trying to write a piece of code in R that

• finds the most-fitting distribution to a set of data, by
• performing goodness-of-fit tests to a list of distributions, and then
• finding the most fitting one.
• This program should be able to run in real-time, so analysis should be very light on computational load. What I mean by this is that it should be able to process, say, a fit every second or few seconds at the most, so the simpler the program, the better.

For instance, I've already written the following code:

for(i in 1:numfit) {
if(distrib[[i]] == "negative binomial"){
  gf_shape = "negative binomial"
  fd_nb <- tryCatch((fitdistr(data, "negative binomial", start=list(size=1,prob=0.5))),
    error = function(fd_nb) fd_nb <- fitdistr(data, "negative binomial"))
  est_size = fd_nb$estimate[[1]]
      est_prob = fd_nb$estimate[[2]]
  gfn = goodfit(data,type="nbinomial",method="MinChisq",par = list(size = est_size))
  tidied = tidy(summary(gfn))
  results[i,] = c(gf_shape, est_lambda, "NA", tidied$X.2, tidied$P...X.2.)
}

else if(distrib[[i]] == "poisson"){
  gf_shape = "poisson"
  fd_p <- fitdistr(data, "poisson")
  est_lambda = fd_p$estimate[[1]]
  gf = goodfit(data,type="poisson",method="MinChisq",par = list(lambda = est_lambda))
  tidied = tidy(summary(gf))
  results[i,] = c(gf_shape, est_lambda, "NA", tidied[1,1], tidied[1,3])
}
results = rbind(c("distribution", "parameter 1", "parameter 2", "chi-squared test statistic", "P > X2"),   results)
return(results)
}

that performs a chi-square goodness-of-fit test of my data to a Poisson and negative binomial distribution, from which I can then find the distribution with the lowest chi-squared test statistic and infer the most suitable distribution to the data from there.

My question is how to do this with continuous data/continuous distributions. Using the goodfit package in R worked really well for me, but it only works with discrete distributions. I do, however, need to use the chi-square goodness-of-fit method (project requirement), so I'm not sure which package or method to turn to and how to implement this.

Can anyone help me with an easy idea of how to implement something similar for, say, a normal distribution? At the moment I think I'm just going to write a piece of code to bin the data into categories, the way you would do the chi-square test by hand, but any way to optimize this process would help. Help would be much appreciated.

Answer 1

It's a pity that a chi-square goodness-of-fit test is a "project requirement," because in general there is much to be lost by binning continuous variables. If possible, try to convince those in charge to allow methods more appropriate for continuous variables in this context of distribution fitting.

If you are stuck with binning you might have a problem. Essentially, if you define the bin boundaries based on parameters of a distribution estimated from the data, then the chi-square test itself may no longer be valid. This is beyond the simple loss of degrees of freedom in chi-square due to estimating parameters; the chi-square statistic as usually calculated may be no longer distributed as chi-square. See the "Level 3" section of the answer by @cardinal on that linked page; that answer includes some recommendations for how to proceed (which are beyond my particular expertise). As you proceed, you do have to consider what you are trying to accomplish by fitting different types of distributions to the same data, and to recognize that results might be hard to generalize even if you find a reliable way to define useful bin boundaries.

Answer 2

One thing is to perform a goodness-of-fit test and another to use it to select an appropriate distribution for the data. For the former, I wrote the following function:

#' Chi-squared Goodness-of-fit Test for Continuous Data
#'
#' Performs the chi-squared goodness-of-fit test for a continuous distribution
#' by grouping data into *bins* (also called *intervals*, *classes* or *cells*)
#' with equal probabilities under the null hypothesis.
#' @param x numeric vector containing the observed values.
#' @param distribution character string naming a continuous distribution,
#' such as "norm" or "unif" (the cumulative distribution `p<distribution>()` and
#' quantile `q<distribution>()` functions must exist).
#' @param nclasses number of bins.
#' @param output logical; if `TRUE` an histogram is plotted and a table with
#' the results for each class is printed.
#' @param nestpar number of estimated parameters (composite null hypothesis).
#' @param ...  parameters of the distribution (specified by `distribution`).
#' @return A list with class `"htest"` containing the following components:
#' \item{statistic}{the value the chi-squared test statistic.}
#' \item{parameter}{the degrees of freedom of the approximate chi-squared
#' distribution of the test statistic.}
#' \item{p.value}{the p-value for the test.}
#' \item{method}{a character string indicating the type of test performed.}
#' \item{data.name}{a character string with the actual `x` argument name.}
#' \item{classes}{a character vector with the class labels.}
#' \item{observed}{the observed counts.}
#' \item{expected}{the expected counts under the null hypothesis.}
#' \item{residuals}{the Pearson residuals, `(observed - expected) / sqrt(expected)`.}
#' @examples
#' nx <- 30
#' x <- rnorm(nx)
#' chisq.cont.test(x, distribution = "norm", nestpar = 2,
#'                 mean = mean(x), sd = sqrt((nx - 1) / nx) * sd(x))
#' @export
chisq.cont.test <- function(x, distribution = "norm",
                            nclasses = floor(length(x)/5), output = TRUE, nestpar = 0, ...) {
  # Funciones distribución
  q.distrib <- eval(parse(text = paste("q", distribution, sep = "")))
  d.distrib <- eval(parse(text = paste("d", distribution, sep = "")))
  # Puntos de corte
  q <- q.distrib((1:(nclasses - 1))/nclasses, ...)
  tol <- sqrt(.Machine$double.eps)
  xbreaks <- c(min(x) - tol, q, max(x) + tol)
  # Gráficos y frecuencias
  if (output) {
xhist <- hist(x, breaks = xbreaks, freq = FALSE,
              lty = 2, border = "grey50")
curve(d.distrib(x, ...), add = TRUE)
  } else {
xhist <- hist(x, breaks = xbreaks, plot = FALSE)
  }
  # Cálculo estadístico y p-valor
  O <- xhist$counts  # Equivalente a table(cut(x, xbreaks)) pero más eficiente
  E <- length(x)/nclasses
  DNAME <- deparse(substitute(x))
  METHOD <- "Pearson's Chi-squared test"
  STATISTIC <- sum((O - E)^2/E)
  names(STATISTIC) <- "X-squared"
  PARAMETER <- nclasses - nestpar - 1
  names(PARAMETER) <- "df"
  PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
  # Preparar resultados
  classes <- format(xbreaks)
  classes <- paste("(", classes[-(nclasses + 1)], ",", classes[-1], "]",
                   sep = "")
  RESULTS <- list(classes = classes, observed = O, expected = E,
                  residuals = (O - E)/sqrt(E))
  if (output) {
    cat("\nPearson's Chi-squared test table\n")
    print(as.data.frame(RESULTS))
  }
  if (any(E < 5))
    warning("Chi-squared approximation may be incorrect")
  structure(c(list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL,
                   method = METHOD, data.name = DNAME), RESULTS), class = "htest")
}

which could be used as a basis for the latter...