gistfile1.txt

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

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, nosim = 1000){
  ## Just useful for sampling 
  rows = 1 : nrow(A)
  cols = 1 : ncol(A)
  
  ## The observed test statsitic
  tsobs = testStat(A)
  
  ## The test condition
  test1 = c(1, 0, 0, 1)
  test2 = c(0, 1, 1, 0)

  tsCurrent = tsobs
  
  pvalue = 1 / nosim
  
  for (i in 1 : nosim){
      i = sample(rows, 2)
      j = sample(cols, 2)
      ## the tetrad subtable
      subTab = A[i, j]
      ## check if it's one of the two types
      if (all(subTab == test1)) {
        A[i, j] = test2
        tsCurrent = testStat(A)
      }
      else if (all(subTab == test2)) {
        A[i, j] = test1
        tsCurrent = testStat(A)
      }
      #else no change
    pvalue = pvalue + (tsCurrent >= tsobs) / nosim
  }  
  return(pvalue)
}