[R] Cochrans Q Test

[Marc Schwartz]

On Mon, 2006-09-18 at 14:01 +1000, peter.gremse at numbers.net.au wrote:
> Hi!
> 
> I would like to conduct a Cochran`s Q Test in R, but have not found any 
> suitable function.
> 
> My first idea was: J <- as.table(matrix(c(6,16,2,4),ncol=2, dimnames = 
> list("C" = c("Favorable","Unfavorable"),"Drug A Favorable"=c("B 
> Favorable","B Unfavorable"))))
> L <- as.table(matrix(c(2,4,6,6),ncol=2, dimnames = list("C" = 
> c("Favorable","Unfavorable"),"Drug A Unfavorable"=c("B Favorable","B 
> Unfavorable"))))
> mantelhaen.test(J,L, alternative="t")
> 
> But this is obviously the wrong function.

The CMH test does not consider the dependent nature of the measurements,
so it is indeed the wrong test, if you have dependent measures as your
data 'K' below would suggest.

Cochran's Q is a generalization of the McNemar paired chi square.

> Then I googled and found (different data):
> 
> K <- as.table(matrix(c(1,1,0,0, 1,1,0,1, 1,1,1,1, 1,1,1,1, 1,0,0,0, 
> 1,1,1,1, 1,1,1,1, 0,0,0,0, 0,1,0,1, 1,1,1,1, 0,1,0,1, 0,1,0,1),ncol=12, 
> dimnames = list("Seating type" = c("I","II","III","IV"),"Test 
> subject"=c("A","B","C","D","E","F","G","H","I","J","K","L"))))
> K
> pcochran(K,4,12)
> 
> But R said that this function does not exist.
> Can anyone help?

Here is a version of the Cochran's Q that I wrote and posted to
sci.stat.consult back in June in response to a thread there. This is
based upon Sheskin's Handbook of Parametric and Nonparametric
Statistical Procedures (2004) page 867. You might want to secure a copy
of the book and review the comments regarding the assumptions underlying
this test and the considerations for it use.

cochranq.test <- function(mat)
{
  k <- ncol(mat)

  C <- sum(colSums(mat) ^ 2)
  R <- sum(rowSums(mat) ^ 2)
  T <- sum(rowSums(mat))

  num <- (k - 1) * ((k * C) - (T ^ 2))
  den <- (k * T) - R

  Q <- num / den

  df <- k - 1
  names(df) <- "df"
  names(Q) <- "Cochran's Q"

  p.val <- pchisq(Q, df, lower = FALSE)

  QVAL <- list(statistic = Q, parameter = df, p.value = p.val,
               method = "Cochran's Q Test for Dependent Samples",
               data.name = deparse(substitute(mat)))
  class(QVAL) <- "htest"
  return(QVAL)
}


Using your matrix 'K':

> K
            Test
subject
Seating type A B C D E F G H I J K L
         I   1 1 1 1 1 1 1 0 0 1 0 0
         II  1 1 1 1 0 1 1 0 1 1 1 1
         III 0 0 1 1 0 1 1 0 0 1 0 0
         IV  0 1 1 1 0 1 1 0 1 1 1 1


> cochranq.test(K)

        Cochran's Q Test for Dependent Samples

data:  K
Cochran's Q = 23.9298, df = 11, p-value = 0.01303


BTW, a quick Google search shows that the pcochran() function is in the
'outliers' package on CRAN, which also has a cochran.test() function

HTH,

Marc Schwartz