Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save muschellij2/70de5d3b5fc25c0ddcfc0971f50d4c71 to your computer and use it in GitHub Desktop.
Save muschellij2/70de5d3b5fc25c0ddcfc0971f50d4c71 to your computer and use it in GitHub Desktop.
## 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment