Nick Cutler <s0455078 <at> sms.ed.ac.uk> writes: > > I would like to be able to randomise presence-absence (i.e. binary) > matrices whilst keeping both the row and column totals constant. Is > there a function in R that would allow me to do this? > > I'm working with vegetation presence-absence matrices based on field > observations. The matrices are formatted to have sites as rows and > species as columns. The presence of a species on a site is indicated > with a 1 (absence is obviously indicated with a 0). > > I would like to randomise the matrices many times in order to construct > null models. However, I cannot identify a function in R to do this, and > the programming looks tricky for someone of my limited skills. > > Can anybody help me out?
Nick, For a 1001 x 1001 matrix, this method takes less than 2 seconds on my 2 year old Windows PC. ronetab( marg1, marg2 ) returns a table of 0's and 1's according to the marginal contraints. ck.ronetab( marg1, marg2 ) checks that all the constraints were honored. msample <- function(x,marg) { ## Purpose: sample at most one each from each cell of a margin ## ---------------------------------------------------------------------- ## Arguments: x - number to sample, marg - a vector of integers ## ---------------------------------------------------------------------- ## Author: Charles C. Berry, Date: 28 Apr 2007, 08:17 ## GPL 2.0 or better if (!(x<=sum(marg) && all(marg>=0))) browser() wm <- which(marg!=0) if (length(wm)==1) { wm } else { sample( seq(along=marg), x, prob=marg ) } } ronetab <- function(m1,m2,debug=F) { ## Purpose: sample from a table with fixed margins and {0,1} cell values ## ---------------------------------------------------------------------- ## Arguments: m1, m2 - two margins ## ---------------------------------------------------------------------- ## Author: Charles C. Berry, Date: 28 Apr 2007, 08:21 ## GPL 2.0 or better stopifnot( sum(m1)==sum(m2)|| max(m1)>length(m2) || max(m2)>length(m1) ) i.list <- j.list <- list() k <- 0 while( sum(m1)>0 ){ k <- k+1 if ( sum(m1!=0) > sum(m2!=0) ){ i <- which.max( m1) j <- msample( m1[i], m2 ) i.list[[ k ]] <- rep( i, m1[i] ) j.list[[ k ]] <- j m1[i] <- 0 m2[ j ] <- m2[ j ] - 1 } else { j <- which.max( m2 ) i <- msample( m2[j], m1 ) i.list[[ k ]] <- i j.list[[ k ]] <- rep( j, m2[j] ) m2[j] <- 0 m1[ i ] <- m1[ i ] - 1 } } res <- array(0, c(length(m1), length(m2) ) ) res[ cbind( unlist(i.list), unlist(j.list) ) ] <- 1 res } ck.ronetab <- function(m1,m2){ tab <- ronetab(m1,m2) m1.ck <- all(m1==rowSums(tab)) m2.ck <- all(m2==colSums(tab)) cell.ck <- all( tab %in% 0:1 ) res <- m1.ck & m2.ck & cell.ck if (!res) attr(res,"tab") <- tab res } I'll warn you that I have not worked through what looks to be a tedious proof that this randomly samples matrices under the constraints. The heuristics seem right, and a few simulation spot checks look reasonable. If you do not want to trust it, you can still use it to generate a starting value for an MCMC run. HTH, Chuck ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.