N-1 Pearson's Chi-Square in R

Question

I have two categorical variables and was looking into doing a chi-square test. I then noticed I had some low frequencies in my contingency table and thought Fisher's Exact Test may be useful. I've now come full circle after doing some reading and want to use Pearson's Chi Squared with n-1 correction. Is there a way in R to run chisq.test with the n-1 correction (discussed here: Given the power of computers these days, is there ever a reason to do a chi-squared test rather than Fisher's exact test?)?

If not, how would I apply the correction to the output of Pearson's chi-squared?

Presuming a sample size of 80:

(80-1)/80 = 0.9875

Do I simply multiply the Chi-Squared statistic by 0.9875 and then use this value to derive the p value?

2.9687 * 0.9875 = 2.931591

1-pchisq(2.931591,4)

p = 0.569338

Answer 1

According to this page the N-1 correction is very simple; just multiply $\chi^2$ by (N-1)/N. You could then use the pchisq function in R to get the right p value (the exact code would be, I believe, something like

newchisq = ((N-1)/N) * oldchisq
newp <- 1 - pchisq(newchisq, df)

Answer 2

Following the above discussions I wrote an R function implementing recommended policy from cited here Campbell (2007) paper which is:

1. Where all expected numbers are at least 1, analyse by the ‘N − 1’ chi-squared test (the K. Pearson chi-squared test but with N replaced by N − 1).
2. Otherwise, analyse by the Fisher–Irwin test, with two-sided tests carried out by Irwin’s rule (taking tables from either tail as likely, or less, as that observed).

campbell2x2.test <- function(t) {
  min_exp_val <- min(colSums(t))*min(rowSums(t))/sum(t)

  if (min_exp_val < 1) {
    # in Campbell's naming: Fisher–Irwin test by Irwin’s rule
    result <- fisher.test(t, alternative = "two.sided")
    result$method <- paste("Optimal 2x2 test according to Campbell(2007) recommendation\n\n",
                       paste("Minimal expected cell count: ", round(min_exp_val, 3), "\n\n", sep = ""),
                       paste("Performing", result$method)
    )
    return(result)
  } else {
    #  'N − 1' Pearson's Chi-squared test
    n1chisq.test <- function(t) {
      chisqtst <- chisq.test(t, correct = FALSE)
      N <- sum(chisqtst$observed)
  chisqtst$statistic = ((N-1)/N) * chisqtst$statistic
  chisqtst$p.value <- 1 - pchisq(chisqtst$statistic, chisqtst$parameter)
      chisqtst$method <- paste("'N-1'", chisqtst$method)
      return(chisqtst)
    }
    result <- n1chisq.test(t)
    result$method <- paste("Optimal 2x2 test according to Campbell(2007) recommendation\n\n",
                       paste("Minimal expected cell count: ", round(min_exp_val, 3), "\n\n", sep = ""),
                       paste("Performing", result$method)
    )
    return(result)
  }
}

Example usage:

campbell2x2.test(matrix(c(1, 5, 3, 2), nrow = 2))

Answer 3

Sorry for showing up so late to this party. I just happened across this discussion and saw that some of my notes were being discussed. I gather that this is the section of the notes that was causing some confusion for some folks:

Where Campbell describes replacing N with N -1, he is referring to this formula for Pearson's chi-square:

        chi-square = N(ad-bc)^2 / (mnrs)

where:

N is the total number of observations
a, b, c, and d are the observed counts in the 4 cells
^2 means "squared"
m, n, r, s are the 4 marginal totals

If one has the regular Pearson chi-square (e.g., in the output from statistical software), it can be converted to the 'N - 1' chi-square as follows:

           'N -1' chi-square = Pearson chi-square x (N -1) / N

If you want N(ad-bc)^2 / (mnrs) to become (N-1)(ad-bc)^2 / (mnrs), surely you must divide the whole expression by N and then multiply by (N-1).

I hope this clarifies things.

Answer 4

Busing, Weaver and Dubois (2015) (http://onlinelibrary.wiley.com/doi/10.1002/sim.6808/full) give good advice on how to implement the N-1 chi-squared test in different software packages.

One possibility is to use the uncorrected Mantel-Haenszel chi-squared statistic which is equivalent to using N-1 chi-squared when you are analyzing a single 2x2 table.

In R, the uncorrected Mantel-Haenszel chi-squared can be obtained as follows:

stratum <- rep(1, N)
mantelhaen.test(variable1, variable2, stratum, correct = FALSE)

For the code to work as intended, the elements in the stratum vector must all have the same value.