Implements Mantel's chi-square test for ordinal categorical data

test.mantel.R

test.mantel <- function (counts, rowScores = NULL, colScores = NULL) {
  ## This implements an ordinal chi-square test following Mantel, 1963,
  ## as described in Agresti, 2007, section 2.5.1 "Linear Trend Alternative to
  ## Independence"
  
  METHOD <- "Mantel's Ordinal Chi-Square Test"
  DNAME = deparse(substitute(counts))
  
  n = sum(counts) # total obs
  if (n < 30)
    warning("Chi-squared approximation may be incorrect (n < 30)")
  
  ## marginals
  wt_c = colSums(counts)
  wt_r = rowSums(counts)
  
  ## default to uniform spacing
  if (is.null(rowScores))
    rowScores <- seq(nrow(counts))
  if (is.null(colScores))
    colScores <- seq(ncol(counts))
  
  PARAMETERS <- c(rowScores, colScores)
  names(PARAMETERS) <- c(paste0("rowScore_", rowScores), 
                         paste0("colScore_", colScores))
  
  ## score means
  mu_r = sum(wt_r*rowScores)/n
  mu_c = sum(wt_c*colScores)/n
    
  ## mean-centered scores
  rowScores.norm = rowScores - mu_r
  colScores.norm = colScores - mu_c
  
  ## weighted covariance
  wt_cov = sum(outer(rowScores.norm, colScores.norm) * counts)
  
  ## weighted variances
  wt_rowVar = sum(rowScores.norm * rowScores.norm * wt_r)
  wt_colVar = sum(colScores.norm * colScores.norm * wt_c)
  
  ## weighted correlation
  r = wt_cov / sqrt(wt_rowVar*wt_colVar)
  
  ## test statistic ~ ChiSquare(df=1)
  STATISTIC = (n - 1) * r**2
  names(STATISTIC) <- "M-square ~ Chi-square(df=1)"
  
  PVAL = pchisq(STATISTIC, df = 1, lower.tail = FALSE)
  
  structure(list(statistic = STATISTIC, p.value = PVAL, method = METHOD,
                 observed = counts, parameters = PARAMETERS, data.name = DNAME),
            class = "htest")
}

## Comparison with Pearson's chi-square test
test.mantel(matrix(c(100,82,82,100), nrow = 2, ncol = 2), c(1,2), c(1,2))
chisq.test(matrix(c(100,82,82,100), nrow = 2, ncol = 2), correct = FALSE)

## 6 x 4 Example with uneven spacing
cts = cbind(c(10,15,16,17,15,10),
            c(0,0,0,0,0,0),
            c(20,15,10,5,5,6),
            c(5,5,5,5,5,5))
rscores = c(1,2,3,4,5,6)
cscores = c(5, 10, 80, 90)

test.mantel(cts, rscores, cscores)

## Contrast Two Groups
test.mantel(cts, c(0,0,0,1,1,1), cscores)

## Note that the test statistic is invariant to linear transformations
## of the scores. On the other hand, non-linear transformations can 
## significantly impact the result.