# Cochran's Q for dichotomous variates
# Author: M.V.Wickerhauser
# Date: 28 February 2007
#
# Input: array of 0s and 1s
# Example:
# > data <- c(1,1,1,0, 0,1,0,1, 0,0,1,0)
# > table <- array(data, c(4,3))
# > table
# [,1] [,2] [,3]
# [1,] 1 0 0
# [2,] 1 1 0
# [3,] 1 0 1
# [4,] 0 1 0
#
# [Note that as "data" is read from left to right,
# the first, or row, index in "table" changes fastest. ]
#
# Then call "cochranQ(table)" to get
# Q, for blocks {1,1,1,0}, {0,1,0,1}, and {0,0,1,0}
# Qt, for blocks {1,0,0}, {1,1,0}, {1,0,1}, and {1,0,0}
# and their degrees of freedom and p-values (from Chi squared density).
#
cochranQ <- function ( table ) {
a <- length(table[1,]) # number of blocks
A <- array(0, a)
for (i in 1:a) {
A[i] <- sum(table[,i])
}
b <- length(table[,1]) # number of groups
B <- array(0, b)
for (j in 1:b) {
B[j] <- sum(table[j,])
}
Q <- (b-1)*(sum(B^2) - sum(B)^2/b)/(sum(A) - sum(A^2)/b);
Qt <- (a-1)*(sum(A^2) - sum(A)^2/a)/(sum(B) - sum(B^2)/a);
cat("\nCochran's Q\n")
cat("Q =\t",Q ,"\t\tDF =",b-1,"\t\tp =",1-pchisq(Q,b-1),"\n")
cat("Qt =\t",Qt,"\t\tDF =",a-1,"\t\tp =",1-pchisq(Q,a-1),"\n")
}