gistfile1.txt

## Performs MCMC P-value calculation for 0-1 tables
## from Darwin's Finch data

set.seed(1)
I = 7; J = 5

testDat = matrix(sample(0:1, I * J, replace= TRUE), I, J)

A = matrix(c(
  0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0,
  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0,
  0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1,
  1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
  0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
  0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,
  0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0,
  0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
), 13, 17, byrow = TRUE
)


testStat = function(A){
  m = nrow(A)
  AAT = A %*% t(A)
  2 * sum( AAT[upper.tri(AAT)] ^ 2 ) / m / (m - 1)
}

zeroOne = function(A, test = TRUE, nosim = 1000){
  ## Error checking
  if (test) {
    ## A has to be a matrix
    stopifnot(is.matrix(A))
    ## A has to have only 0s and 1s
    stopifnot(all(A %in% (0 : 1)))
  }
  
  ## Get rid of zero rows and columns
  A = A[rowSums(A) != 0, colSums(A) != 0]
  
  
  rows = 1 : nrow(A)
  cols = 1 : ncol(A)
  
  tsobs = testStat(A)
  
  ## The first observation
  ts = NULL
  
  test1 = c(1, 0, 0, 1)
  test2 = c(0, 1, 1, 0)
  flip = sample(c(TRUE, FALSE), size = nosim, replace = TRUE)
  
  for (sim in 1 : nosim) {
    if (flip[sim]) {
      i = sample(rows, 2)
      j = sample(cols, 2)
      
      subTab = A[i, j]
      
      #print(subTab)
      
      ## Flip the coin regardless, only
      ## check if you were going to swap it anyway
      if (all(subTab == test1))
        A[i, j] = test2
      else if (all(subTab == test2)) {
        A[i, j] = test1
      }
    }
    ts = c(ts, testStat(A))
  }
  pvalue = (ts >= tsobs)
  pvalue = cumsum(pvalue) / (1 : length(pvalue))
  
  return(data.frame(statistic = ts, p.value = pvalue, 
                    stringsAsFactors = FALSE))
}

res = zeroOne(A)
plot(res$p.value)
tail(res)