In 1975, Diamond  proposed rules which govern the combination of several bird communities on islands in New Guinea, such as “Some pairs of species never coexist, either by themselves or as a part of a larger combination”. However, Connor & Simberloff  disputed these rules, suggesting that the combination of species on the islands was not exceptionally different from random. Aside from the ensuing feud within the ecological community [9, 2, 26, 20, 25, 12, 13], Connor & Simberloff sparked interest in a novel form of “null model analysis” [16, 14], a general approach which compares an observed matrix statistic to the distribution of that statistic under a null model. The null model of choice in this case was the set of all random combinations of species on islands satisfying the following three constraints:
The number of species on each island matches the number observed.
The number of islands inhabited by each species matches the number observed.
Each species inhabits islands of comparable size to the islands on which it is observed.
Connor & Simberloff reason that if the observed species combination is not substantially different from the random combinations with respect to some statistic of interest, then it becomes difficult to argue that the observed pattern is anything noteworthy.
To generate combinations which adhere to the above constraints, they represent the species/ island combinations as a matrix, with rows representing species and columns representing islands. The element in a particular row and column is set to 1 if the corresponding species inhabits the corresponding island, and 0 otherwise. In the matrix representation, the first two constraints imply fixing the marginal sums of the matrix, while the third constraint limits which elements are permitted to be 1. Thus, the aim in the “null model” analysis is to generate a collection of matrices which possess the desired marginal sums and exclude the prohibited 1’s. Connor & Simberloff propose two methods, one which begins with a 0-filled matrix and sequentially fills in each row randomly, being careful to adhere to the constraints, and another which begins with the observed matrix, and repeatedly swaps 1’s such that the constraints are maintained. Subsequent authors have proposed additional methods for generating the desired matrices, which are generally either a “fill” [28, 26, 25, 13] or “swap” [2, 24, 26, 20, 29] type method.
Early debates surrounding the use of this particular null model focused on the appropriateness of imposing Connor & Simberloff’s constraints, the choice of statistic upon which to base comparison, and the interpretability of null model tests [9, 24].
However, focus eventually shifted from whether to use this approach to which method of producing matrices is best .
This accompanied theoretical work in binary matrices with fixed marginal sums  and Markov Chain Monte Carlo (MCMC) sampling
and Markov Chain Monte Carlo (MCMC) sampling, paving the way for swap and fill methods with proven statistical properties [19, 5].
Generally speaking, sampling from a null model can be useful for quantifying the degree of “non-randomness” in an observed matrix (although care should be exercised in the selection of a test statistic and mechanistic interpretation of results).
This approach has been used for testing for significance of species co-occurrence , and in the analysis of contingency tables
Generally speaking, sampling from a null model can be useful for quantifying the degree of “non-randomness” in an observed matrix (although care should be exercised in the selection of a test statistic and mechanistic interpretation of results). This approach has been used for testing for significance of species co-occurrence, evaluating goodness of fit for the Rasch item-response model [2, 5]
, and in the analysis of contingency tables.
In this paper, we extend the capability of null model analysis by allowing for non-uniform sampling. This permits the construction of more sophisticated “null” models corresponding to non-uniform probability distributions.
1.1 Review of Sampling Methods
Here we briefly review the variety of existing swap and fill methods that have been employed for sampling binary matrices with fixed margins. To highlight the contribution of this work, we point out that most efforts at sampling these matrices thus far have focused on uniform sampling, where each matrix has an equal probability of being included in the sample. Early work gave little attention to actually proving uniformity and few authors address situations in which some matrix elements must be 0, so called structural zeros, with exceptions noted below.
As mentioned, these sampling methods can be categorized into fill and swap methods. Fill methods are typically faster than swap methods, but risk hitting “dead ends”, where the process must restart or backtrack in order to avoid violating row and column margin constraints. Swap methods, on the other hand, are arguably simpler to understand and implement, but rely on local changes, so sampled matrices tend to be highly correlated. For both categories, initial methods lacked much theoretical treatment, but more recent work has considered the statistical properties of the algorithms.
1.1.1 Fill Methods
Fill methods begin with a matrix of all 0’s, and randomly fill in 1’s in such a way that adheres to the margin (row and column) constraints.
Connor & Simberloff’s  first placed the row sums in decreasing order and column sums in increasing order (“canonical form”), then proceeded row by row, randomly selecting columns to fill with 1’s, until the desired row sum is reached. If the column is already 1, or it has met its column sum, another random selection is made. If the procedure reaches a point where no 1 can be filled without violating some constraint (a dead end), the process restarts with an empty matrix. No mention is made of the probability of generating a particular matrix. In their “Milne method”, Stone & Roberts  adapt the fill method of Connor & Simberloff to check whether the proposed fill will lead to a dead end.
Sanderson et al  propose a “knight’s tour” algorithm which randomly samples a row and column index for filling, rather than filling one row at a time, and backtracks rather than restarts if the algorithm reaches a dead end. They claim that since row and columns are sampled randomly, then each valid matrix occurs with equal probability. However, Gotelli & Entsminger show that the “knight’s tour” algorithm does not sample uniformly, and propose a “random knight’s tour” which only samples a subset of available elements before backtracking, and backtracks randomly rather than sequentially. They provide no proof of uniform sampling, instead showing on a small example that a sample generated by the random knight’s tour does not reject a lack of fit test.
Chen et al  propose a sequential importance sampling approach which randomly fills in each row, according to a conditional Poisson proposal distribution that depends on the previously filled rows.
The matrices in the sample are then re-weighted to approximate a uniform distribution.
This method also approximates the size of the space of valid matrices.
propose a sequential importance sampling approach which randomly fills in each row, according to a conditional Poisson proposal distribution that depends on the previously filled rows. The matrices in the sample are then re-weighted to approximate a uniform distribution. This method also approximates the size of the space of valid matrices.
Taking a different approach, Miller & Harrison  devise a recursive sampling algorithm based on the Gale-Ryser Theorem to generate matrices from the uniform distribution exactly, and compute the size of the space exactly. The recursion is over row and column sums, so this method is limited to matrices of moderate size.
Finally, Harrison & Miller  adopt the sequential importance sampling approach of Chen et al for non-uniform sampling. This is, to our knowledge, the only paper dedicated towards deliberate non-uniform sampling. They propose a multiplicative weights model that forms the basis of our weights model below, and permit specifying structural zeros by setting a corresponding weight to 0.
1.1.2 Swap Methods
Swap methods have developed in parallel to fill methods, beginning with the original paper of Connor & Simberloff . These methods involve identifying “checkerboard” patterns of 0’s and 1’s in an existing valid matrix, and swapping the 1’s for the 0’s to produce a mew matrix which still adheres to the constraints (see Figure 1). Repeated swaps result in a Markov chain Monte Carlo algorithm that produces a dependent sample of valid matrices. Of particular importance is the method used to identify potential swaps, and what happens if a proposed swap is invalid, which have computational as well as theoretical implications. Initially, there was some question about whether swapping methods could produce all possible valid matrices by starting at the observed matrix; in other words, if the space of valid matrices is connected under checkerboard swaps. However, Brualdi  showed this to be true using elementary circuits. This paved the way for more development in swap methods [26, 20].
As with fill methods, early swap methods gave no proof that matrices are sampled uniformly, until Kannan et al  proved this for the following procedure: Select two rows and two columns uniformly at random. If they form a checkerboard, then perform a swap, and if not, keep the previous matrix. This result influenced the development of subsequent swapping algorithms. Gotelli & Entsminger’s  version of the swapping algorithm importantly picks a swap randomly from the set of eligible swaps with equal probability, not as prescribed by Kannan et al, and consequently demonstrate that the resulting distribution which this samples from is not uniform. This approach is also taken by Zaman & Simberloff , who correct for the imbalance by reweighting each matrix when computing summary statistics. These approaches require identifying all possible swaps at each step in order to sample uniformly from them, which is prohibitive for large matrices.
Finally, in the “curveball” algorithm proposed by Strona et al  multiple swaps are performed at once by first selecting two rows (or columns), and permuting the elements in their symmetric difference. Following the proof technique of Kannan et al., Carstens  proves that curveball samples uniformly and demonstrate that multiple swaps can lead to faster mixing compared to Kannan et al’s algorithm.
1.2 Non-Uniform Sampling
Most existing null matrix sampling algorithms focus on the uniform distribution, the goal being to compare the observed matrix to other matrices with the same marginal properties, without preference for particular matrices over others. We consider, however, that for some data analyses we may wish to lend more importance to some matrices compared to others, when testing for a specific effect. For instance, Diamond  mentions several factors which may give rise to certain species combinations, such as resource overlap, dispersal ability, proneness to local extinction, and distributional strategy, and criticised Connor & Simberloff  for not accounting for such “significant structuring forces”. We argue that it may be possible to control for such external effects in a modified null model, which does not give equal probability to all possible matrices.
Here we consider the task of non-uniform sampling of matrices in a general setting. Following the formulation of Harrison & Miller , we assume the existence of a weight matrix indicating the relative likelihood of a 1 in each element and define the probability of a matrix as a function of these weights. Careful attention must be paid to the consequences of weights equal to 0, which further restrict the set of matrices with positive probability, and affect the validity of the proposed sampling methods.
In the rest of this paper, we describe a weight-based probability model, then present a novel MCMC checkerboard swap algorithm which incorporates the weights and samples according to the given model. The weights introduce the notion of a structural zero, which can impact the efficiency and even correctness of sampling. We then present a novel curveball MCMC algorithm, which in most circumstances is a faster alternative to the weighted checkerboard swap algorithm. We investigate the correctness, efficiency under different weighting and structural zero schemes via simulation. Finally, we revisit the bird species-island dataset first discussed by Diamond  and Connor & Simberloff  and consider incorporating a weight matrix in the null distribution to show the impact that non-uniform sampling can have on scientific conclusions.
2 Non-Uniform Sampling from binary Matrices
2.1 Probability Model
Let be the set of all binary matrices satisfying row sums and column sums . Furthermore let be a non-negative weight matrix representing the relative likelihood of a 1 in element . Following Harrison & Miller , define the probability of observed matrix as
We call the set of elements structural zeros, because any matrix with positive probability must have if . Let be the set of matrices with positive probability. We interpret the weights by considering the relative probability between two matrices :
which is governed by the weights of the elements not shared between them. Furthermore, we can write the conditional probability:
This forms the basis for the swapping probability defined in the next section.
2.2 Weighted Checkerboard Swap Algorithm
0 1 1 0
We propose a swapping algorithm similar to that of Brualdi , which identifies submatrices with diagonal 1’s and 0’s, as in Figure 1, where the rows and columns are not necessarily consecutive. Let and be the indices of the 1’s in the submatrix, so that and are 0’s. The positions of the 1’s and 0’s can be swapped without altering the row or column sums. This swap operation produces a matrix which differs from (see Figure 1). In the uniform sampling case, all matrices have equal probability, so swaps are performed with a fixed probability. In our version, the swap is performed with probability
This modification implies matrices with zero probability will be visited with probability zero.
In Algorithm 1, we present an MCMC procedure based on weighted checkerboard swap. As with any MCMC algorithm, this algorithm can be modified to incorporate burn-in and thinning to reduce autocorrelation between saved samples.
If all weights are positive, then we have the following result:
Given binary matrix and weight matrix with , then Algorithm 1 generates a Markov chain with stationary distribution given by .
2.3 Structural Zeros
As mentioned, structural zeros have received little attention since originally introduced by Connor & Simberloff . We show here that structural zeros can be problematic for swapping algorithms, first addressing sampling efficiency then correctness.
Structural zeros can result in some checkerboard swaps having zero probability, which removes transitions from the Markov chain. This can increase the number of swaps necessary to transition between two states, in turn reducing the sampling efficiency, and is reflected in a larger diameter of the Markov state space. Here we show that Brualdi’s  upper bound on the diameter of the Markov space no longer holds with structural zeros.
Given matrices the upper bound on the maximum number of swaps required to transition from one matrix to the other is given by , where is the Hamming distance between and , and is the largest number of non-overlapping “elementary circuits” in , which we briefly explain here. Consider as a weighted adjacency matrix of a bipartite graph, where the rows and columns form separate partitions of a vertex set. Suppose an entry of indicates an edge from the row vertex to the column vertex, and an element of indicates an edge from column vertex to the row vertex. An elementary circuit is then defined as a cycle of edges in this graph.
To show how the bound can be violated, consider the following matrices, both with unit row and column sums where structural zeros are represented with grey fill:
1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0
The Hamming distance between them is and there is a single elementary circuit, giving a bound of . However, seven swaps are required to transition from to , which can be described as a sequence of column pairs, since there is only one 1 per column: .
Additionally, the following example from Rao et al  illustrates how structural zeros can render it impossible to transition between two matrices via checkerboard swaps:
0 1 0 0 0 1 1 0 0 0 0 1 1 0 0 0 1 0 .
Because no sequence of swaps transforms into , the Markov chain described in Algorithm 1 is reducible and hence there is no unique stationary distribution.
In many cases, the presence of a few structural zeros will not result in a reducible chain, but in the general case structural zeros can be fatal to checkerboard swap based MCMC sampling algorithms. However, let us consider a special case where they are not problematic, specifically where structural zeros satisfy the following:
The structural zeros of a matrix are monotonic if there exists a permutation of the rows and columns of , such that implies for .
This situation can arise when the rows and columns reflect different moments in time.
For example, let
This situation can arise when the rows and columns reflect different moments in time. For example, letrepresent a citation network among academic papers, where a 1 indicates that the row article cites the column article. Suppose the rows and columns are arranged chronologically with increasing row/column index. Clearly, articles cannot cite into the future, so we prohibit some 1’s using structural zeros. In this example, the upper triangle (including the diagonal) of are structural zeros. Furthermore these structural zeros are monotonic as will satisfy the definition if the order of the columns of are reversed.
When structural zeros are monotonic, Algorithm 1 produces a valid sampler:
Given binary matrix and weight matrix where any structural zeros are monotonic, then Algorithm 1 generates a Markov chain with stationary distribution given by .
2.4 Weighted Curveball Algorithm
Strona et al.  proposed a faster version of checkerboard swapping called the curveball, which imagines two rows trading their 1’s like baseball cards. Two rows are randomly selected, and the 1’s in their symmetric difference become candidates for the trade, where several 1’s may move between rows in a trade. In the uniform setting, the proposed trade is always performed. In the non-uniform setting, we modify the algorithm in the same way as for weighted checkerboard swapping, and introduce a probabilistic trade, with probability dependent on the weights of the elements involved in the trade. The entire procedure is presented in Algorithm 2.
Carstens  proves that the unweighted curveball algorithm can be used to sample uniformly, and we contribute the following theorem for the weighted case:
Given binary matrix and weight matrix where any structural zeros are monotonic, then Algorithm 2 generates a Markov chain with stationary distribution given by .
We performed simulations to investigate the behavior of the weighted swap and weighted curveball algorithms. In the first simulation, we verify the weighted checkerboard swap algorithm samples correctly from a non-uniform distribution for a small example. In the second simulation, we compare the mixing performance of both algorithms for different weighting and structural zero schemes. Finally, we investigate the effect of different weighting schemes on the resultant sampling distribution of a summary statistic.
For the second and third simulations, we will make use of a global statistic termed diagonal divergence, which quantifies how far the 1’s of a square matrix are from its diagonal:
where is the number of 1’s in and is the indicator function.
3.1 Small Example
We first consider a small example, such that every Markov state can be enumerated and exact probabilities can be computed. Consider the following binary matrix , and weight matrix :
Including , there are six possible states with row and column sums equal to 1 (see Figure 2). Under a uniform distribution, all matrices would have probability . induces a non-uniform probability distribution among the matrices (see Table 1), where matrices and have lower probability because their non-zero elements correspond to 1’s in the weight matrix. We compare the true distribution to the distribution of a sample produced by weighted swapping. To sample, we used a burn-in of swaps before retaining every th matrix, to generate a sample of matrices.
Table 1 shows agreement between the true and empirical probability distributions. The KL-Divergence between the empirical and true distributions is .
|0.056||0.051 (+/- 0.007)|
|0.222||0.237 (+/- 0.013)|
|0.222||0.217 (+/- 0.013)|
|0.056||0.058 (+/- 0.008)|
|0.222||0.222 (+/- 0.013)|
|0.222||0.215 (+/- 0.014)|
MCMC standard errors.
3.2 Mixing Performance
In MCMC sampling, efficient mixing is desirable in order to maximize inferential capability and minimize computation time. Due to autocorrelation in sequential observations, the precision of a statistic from MCMC sampling is less than it would be for independent, identically distributed (iid) sampling. The difference in precision can be measured by effective sample size (ESS), which gives the number of iid samples required to match the precision of the MCMC sample . We simulated matrices using weighted checkerboard swapping and weighted curveball under different weight schemes and structural zero configurations in order to assess their impact on mixing, as measured by ESS. ESS was computed for the diagonal divergence statistic using the coda package .
To assess the impact of weights, we constructed weight matrices by sampling each element independently from one of the distributions: Exponential(1), Uniform(0,1), and Uniform(0.5, 1), and also a weight matrix of all 1’s, corresponding to uniform sampling. After constructing the weight matrices, we introduced structural zeros with one of two methods. The first method was to randomly set elements equal to 0 with fixed probability , where was either 0.1, 0.25, or 0.5. The resulting matrices will likely correspond with irreducible Markov chains, given the relatively small probabilities. The second method was to designate a lower triangular region of the weight matrix as all zeros, so that they would be monotonic. Triangles covered either 10%, 25%, or 50% of the weight matrix. We generated a binary starting matrix by setting each element to 1 independently with probability 0.25, resulting in a density of 23.5% before enforcing structural zeros.
For each combination of weights and structural zeros, we ran each sampling algorithm for 10,000 iterations (with no burn-in or thinning). Figure 3 shows ESS for the different conditions, where each line is the average over 10 individual chains with different random seeds. For illustration, Figure 4 shows ESS for all chains for the condition in the bottom-right panel of Figure 3.
Comparing the weighted swap algorithm to the weighted curveball algorithm, we see that the curveball algorithm tends to mix more rapidly. This reflects the fact that weighted curveball can perform multiple trades simultaneously, whereas weighted swap only swaps two 1’s at a time, increasing the autocorrelation between sequential observations. This difference is less pronounced for the Exponential(1) and Uniform(0,1) weighting schemes compared to the Uniform(0.5, 1) and “All 1’s” weighting schemes. This is because the former weight schemes contain weights near 0, making 1 in the corresponding position of possible but unlikely. Since curveball attempts to make several trades simultaneously, this increases the chances of attempting to trade a 1 into an unlikely position. This makes the trade itself unlikely, so fewer trades take place overall, removing the advantage that weighted curveball has over weighted swap. The latter two weight schemes do not produce weights near 0 (ignoring structural zeros), so more curveball trades are successful, restoring the ESS gap between the algorithms. One instance where weighted swap may have an advantage over weighted curveball is in sparse matrices with weights near 0. The sparsity would limit the ability of weighted curveball to perform multiple trades, and the low weights would favor the “smaller steps” taken by weighted swap.
Comparing random structural zeros to monotonic structural zeros, we see little difference for the top two weighting schemes. As mentioned before, these two schemes already tend to produce matrices with weights near 0, so setting some elements equal to zero does not change the weight matrices much, whether done randomly or monotonically. However, the bottom two weight schemes produce weights away from 0, so setting some as structural zeros is a bigger change. In these cases, mixing is better when structural zeros are monotonic compared to random, with the difference being more pronounced when there are more structural zeros. This is unsurprising given the discussion in Section 2.3 showing that structural zeros may give rise to reducible chains, whereas for monotonic structural zeros irreducibility is preserved, and the space has bounded diameter (see Appendix C).
Finally, we acknowledge that “well mixing” may be a questionable label for any of the scenarios shown, since ESS is relatively paltry compared to the number of iterations. In practice, this is ameliorated by careful algorithm implementation (which we do not discuss here) and aggressive thinning.
3.3 Effect of Heterogeneous Weights
The matrix adds flexibility to the sampling algorithms by allowing for non-uniform stationary distributions, however it is unclear how much influence has over the resultant sample given the fixed row and column margins. To investigate this issue, we compare the sampling distribution for a statistic under different “strength” weight matrices, which are constructed by taking the elements of a base weight matrix to a positive or negative power. In this case, we let have dimension , and define the base weight matrix as:
which has value on the diagonal and value near in the corners (see Figure 5). Sampling proceeded with a burn-in period of proposed swaps, before retaining every th matrix, to generate a sample of retained matrices.
Figure 6 shows the approximate sampling distributions of the diagonal divergence for selected powers of the weight matrix. Since favors matrices with 1’s near the diagonal, larger powers of favor lower values for , and negative powers of favor higher values for . Larger values of emphasize having 1’s close to the diagonal ( positive) or packed in the corners ( negative), where there are relatively few distinct matrices, and there is a corresponding decrease in the variance of the samples in those cases. This figure illustrates how vast the space
negative), where there are relatively few distinct matrices, and there is a corresponding decrease in the variance of the samples in those cases. This figure illustrates how vast the spaceis as there is virtually no overlap in the sampling distributions under different weighting schemes.
4 Example: New Hebrides Bird Species
We revisit the original Vanuatu (formerly New Hebrides) bird species data provided by Diamond in 1975 , to illustrate the impact weighting can have on the analysis of real data. Of interest is whether or not there exists competitive exclusion between different bird species, or groups of species, on these islands. That is, do certain species or groups exclude one another from islands due to competition over resources? The data are represented in a matrix , where rows are species, and columns are islands (see Figure 6(a)). The elements of the matrix are 1 if the corresponding species is found on the corresponding island, and 0 otherwise. If a checkerboard pattern exists between two species on two islands, this can indicate competitive exclusion (see ). Accordingly, Stone & Roberts  propose the C-score to measure competitive exclusion by counting how many checkerboards exist between species pairs in the matrix:
Here, is the number of shared islands between species and , is the number of islands inhabited by species , and is the number of species. A high C-score is interpreted as evidence of competitive exclusion. We will analyze the data using a null model analysis twice, once with equal weights (as has been done historically), and once with heterogeneous-weights Though we are using real data, this is an illustrative example only and we select a weight matrix to that end.
We define a weight matrix to encode the hypothesis that island preference by species is higher for certain islands than others (see Figure 6(b)). This could result, for example, from several environmental and biological factors. We then proceed, as did Connor & Simberloff , to generate a sample of matrices under the null distribution, assuming equal weights (all 1’s) and then under . We generate samples from each weighting regime, after a burn-in of , and thinning of , using the weighted curveball algorithm.
The resulting distributions of C-scores are shown in Figure 8, along with the empirical upper-tailed p-values of the observed matrix under each model. We see significance at the level under uniform sampling, but not after accounting for the weights. The blocked structure of species-island preferences represented in favors more checkerboard structures between species in the top half and bottom half, such that the observed matrix no longer appears so exceptional. The lesson here is that by incorporating more information about the bird species into our null model, the evidence in favor of competitive exclusion changes, and thus researchers should carefully consider the appropriate weight matrix when performing such analyses.
Given the healthy debate surrounding the appropriate use of null matrix model analysis methods, we recommend that the practitioner weigh this approach against their particular situation and needs.
One advantage of this approach is that by constraining, equivalently conditioning on, the row and column sums of the matrix, we remove the need to explicitly model them.
However, without an estimation procedure for the weights, these values must be specified
Given the healthy debate surrounding the appropriate use of null matrix model analysis methods, we recommend that the practitioner weigh this approach against their particular situation and needs. One advantage of this approach is that by constraining, equivalently conditioning on, the row and column sums of the matrix, we remove the need to explicitly model them. However, without an estimation procedure for the weights, these values must be specifieda priori, and reasonable sensitivity analyses are warranted. Fixing row and column sums is only one way to perform null model analysis, and as Gotelli stated : “There is great value in exploring the results of several null models that incorporate different degrees of randomness”. These types of null models provide information about mechanism through the careful consideration of a summary statistic, but it often can be difficult to argue for the correspondence between mechanism and statistic directly. Therefore, if mechanism is of importance, a direct modeling approach may be more appropriate.
We have introduced weighted versions of the checkerboard swap and curveball algorithms, which extend binary matrix sampling to non-uniform probability distributions. Simulations demonstrate that the choice of weight matrix can impact the efficiency of sampling. In fact, when weights are equal to 0, this can slow the mixing process of the MCMC chain, or even prohibit sampling from the desired distribution. However, monotonic structural zeros can arise in time-evolved networks and are better behaved.
Structural zeros have also received attention from Rao et al  who noted that when structural zeros are on the matrix diagonal, the Markov space may be disconnected when limited to “alternating rectangles”(checkerboard swaps). Their solution was to also consider “alternating hexagons”, which they prove connects the space. Alternative swapping mechanisms like this can be incorporated into the non-uniform setting straightforwardly, by defining the swap probability as a ratio of the relevant weights.
Both the weighted swap and weighted curveball algorithms rely on local changes to perform an MCMC walk, which has implications for sampling efficiency (as seen in Figure 3). An unavoidable consequence is that large numbers of swaps are required to generate near independent samples. Of course, the computational cost of each swap is very small, and several chains can be run in parallel to produce larger samples. In general, the weighted checkerboard swaps is preferred when the matrix is very sparse or very dense, and weighted curveball is preferred when the matrix is closer to 50% dense (which will produce the most trade candidates). The choice of summary statistic also has implications for mixing. A single “global” statistic, like the C-score, mixes faster than a local statistic, like, for example, the marginal probability for an element of the matrix.
Another possible area of application is directed graphs, which are naturally represented as square binary adjacency matrices. Fixing marginal sums of these matrices is equivalent to fixing in-degree and out-degree for each vertex. Unweighted checkerboard swapping has been used to sample uniformly for vertex-labelled networks with self-loops, also known as the Configuration Model . However, we know of no implementation of graph sampling where edges are weighted, as in this paper.
In this work, we do not address how to select appropriate weights, as this should be informed by the scientific setting and specific hypothesis being tested. However, the following observation may help guide the construction of . Consider binary matrices and which differ only by a checkerboard swap in rows and , and columns and :
The ratio of weights communicates the relative probability of the differing elements of two matrices, conditioning on all other elements. Therefore eliciting or estimating the relative propensity of matrix elements could inform the weights specification.
Classic (unweighted) checkerboard swapping and curveball trading, while intuitive and simple to implement, lack the ability to incorporate heterogeneity of probabilities when sampling, rendering them useless in situations where such heterogeneity is known to exist. With our contribution, these algorithms are now equipped to incorporate easy to interpret weights, and to consider more plausible null models, which provides a valuable new tool for matrix analysis.
This work was partially supported by grant IOS-1856229 from the National Science Foundation.
- Besag and Clifford  Julian Besag and Peter Clifford. 1989. Generalized Monte Carlo significance tests. Biometrika 76, 4 (Dec. 1989), 633–642. https://doi.org/10.1093/biomet/76.4.633
Richard A. Brualdi.
Matrices of zeros and ones with fixed row and column sum vectors.Linear algebra and its applications 33 (1980), 159–231.
- Carstens  C. J. Carstens. 2015. Proof of uniform sampling of binary matrices with fixed row sums and column sums for the fast Curveball algorithm. Physical Review E 91, 4 (April 2015), 042812. https://doi.org/10.1103/PhysRevE.91.042812
- Chen et al.  Yuguo Chen, Persi Diaconis, Susan P. Holmes, and Jun S. Liu. 2005. Sequential Monte Carlo methods for statistical analysis of tables. J. Amer. Statist. Assoc. 100, 469 (2005), 109–120.
- Connor et al.  Edward F. Connor, Michael D. Collins, and Daniel Simberloff. 2013. The checkered history of checkerboard distributions. Ecology 94, 11 (Nov. 2013), 2403–2414. https://doi.org/10.1890/12-1471.1
- Connor and Simberloff  Edward F. Connor and Daniel Simberloff. 1979. The Assembly of Species Communities: Chance or Competition? Ecology 60, 6 (Dec. 1979), 1132. https://doi.org/10.2307/1936961
- Diamond  Jared M. Diamond. 1975. Assembly of Species Comunities. In Ecology and Evolution of Communities. Harvard University Press, 342–443.
- Diamond and Gilpin  Jared M. Diamond and Michael E. Gilpin. 1982. Examination of the “null” model of Connor and Simberloff for species co-occurrences on islands. Oecologia 52, 1 (1982), 64–74.
- Flegal and Jones  James M. Flegal and Galin L. Jones. 2010. Batch means and spectral variance estimators in Markov chain Monte Carlo. The Annals of Statistics 38, 2 (April 2010), 1034–1070. https://doi.org/10.1214/09-AOS735
- Fosdick et al.  Bailey K. Fosdick, Daniel B. Larremore, Joel Nishimura, and Johan Ugander. 2018. Configuring Random Graph Models with Fixed Degree Sequences. SIAM Rev. 60, 2 (Jan. 2018), 315–355. https://doi.org/10.1137/16M1087175 Publisher: Society for Industrial and Applied Mathematics.
- Gotelli  Nicholas J. Gotelli. 2000. Null Model Analysis of Species Co-Occurrence Patterns. Ecology 81, 9 (2000), 2606–2621. https://doi.org/10.1890/0012-9658(2000)081[2606:NMAOSC]2.0.CO;2
- Gotelli and Entsminger  Nicholas J. Gotelli and Gary L. Entsminger. 2001. Swap and fill algorithms in null model analysis: rethinking the knight’s tour. Oecologia 129, 2 (Oct. 2001), 281–291. https://doi.org/10.1007/s004420100717
- Gotelli and Graves  Nicholas J. Gotelli and Gary R. Graves. 1996. Null models in ecology. Smithsonian.
- Harrison and Miller  Matthew T. Harrison and Jeffrey W. Miller. 2013. Importance sampling for weighted binary random matrices with specified margins. arXiv:1301.3928 [math, stat] (Jan. 2013). http://arxiv.org/abs/1301.3928 arXiv: 1301.3928.
- Harvey et al.  P H Harvey, R K Colwell, J W Silvertown, and R M May. 1983. Null Models in Ecology. Annual Review of Ecology and Systematics 14, 1 (1983), 189–211. https://doi.org/10.1146/annurev.es.14.110183.001201
Peter D. Hoff.
A first course in Bayesian statistical methods. Springer, London ; New York. OCLC: ocn310401109.
- Jerrum and Sinclair  Mark Jerrum and Alistair Sinclair. 1996. The Markov chain Monte Carlo method: an approach to approximate counting and integration. In Approximation algorithms for NP-hard problems. PWS Publishing Co., Boston. Publisher: PWS Boston.
- Kannan et al.  Ravi Kannan, Prasad Tetali, and Santosh Vempala. 1999. Simple Markov‐chain algorithms for generating bipartite graphs and tournaments. Random Structures & Algorithms 13, 4 (1999), 16.
- Manly  Bryan F. J. Manly. 1995. A Note on the Analysis of Species Co-Occurrences. Ecology 76, 4 (1995), 1109–1115. https://doi.org/10.2307/1940919
- Miller and Harrison  Jeffrey W. Miller and Matthew T. Harrison. 2013. Exact sampling and counting for fixed-margin matrices. The Annals of Statistics 41, 3 (June 2013), 1569–1592. https://doi.org/10.1214/13-AOS1131
- Plummer et al.  Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. CODA: Convergence Diagnosis and Output Analysis for MCMC. R News 6, 1 (2006), 7–11. https://journal.r-project.org/archive/
- Rao et al.  A. Ramachandra Rao, Rabindranath Jana, and Suraj Bandyopadhyay. 1996. A Markov chain Monte Carlo method for generating random (0, 1)-matrices with given marginals. Sankhyā: The Indian Journal of Statistics, Series A (1996), 225–242.
- Roberts and Stone  Alan Roberts and Lewis Stone. 1990. Island-sharing by archipelago species. Oecologia 83, 4 (1990), 560–567.
- Sanderson et al.  James G. Sanderson, Michael P. Moulton, and Ralph G. Selfridge. 1998. Null matrices and the analysis of species co-occurrences. Oecologia 116, 1-2 (1998), 275–283.
- Stone and Roberts  Lewi Stone and Alan Roberts. 1990. The checkerboard score and species distributions. Oecologia 85, 1 (1990), 74–79.
- Strona et al.  Giovanni Strona, Domenico Nappo, Francesco Boccacci, Simone Fattorini, and Jesus San-Miguel-Ayanz. 2014. A fast and unbiased procedure to randomize ecological binary matrices with fixed row and column totals. Nature Communications 5, 1 (June 2014), 1–9. https://doi.org/10.1038/ncomms5114
- Wilson  J. Bastow Wilson. 1987. Methods for detecting non-randomness in species co-occurrences: a contribution. Oecologia 73, 4 (1987), 579–582.
- Zaman and Simberloff  Arif Zaman and Daniel Simberloff. 2002. Random binary matrices in biogeographical ecology—instituting a good neighbor policy. Environmental and Ecological Statistics 9, 4 (2002), 405–421.
Appendix A Proof of weighted checkerboard algorithm
Here we prove that the weighted checkerboard swap algorithm samples from defined in (1) when there are no structural zeros, by showing that the Markov chain implied by the algorithm is ergodic (irreducible and aperiodic) with stationary distribution [19, 18].
Irreducibility. Since we have , so irreducibility follows directly from the result of Brualdi , which showed that a series of swaps could change any matrix into any other matrix when row and column sums are preserved.
Aperiodicity. Since every checkerboard swap has probability less than 1, it is possible to remain in the same state after a step. This is enough to show aperiodicity.
Detailed Balance. Lastly, we show that is the stationary distribution by showing that it satisfies the detailed balance equation:
for any matrices , where and are transition probabilities. If and do not differ by a checkerboard swap, then both transition probabilities are 0, and detailed balance holds trivially. If and differ by a checkerboard swap, then let and be its rows and let and be its columns. Also, let be the sum of the row counts, and define the set , so that for . The transition probability is the probability of first selecting the two 1’s in the checkerboard, , times the probability of performing the swap . Starting from the left side of (10),
so detailed balance holds.
Appendix B Proof of weighted curveball algorithm
Here we prove that the weighted curveball algorithm samples from defined in (1) when there are no structural zeros.
Irreducibility. Any checkerboard swap can be viewed as a curveball trade, so the set of all weighted curveball trades contains the set of all weighted checkerboard swaps. Therefore any state transition in the checkerboard swap Markov chain is also a state transition in the curveball trade chain, hence irreducibility holds for curveball as well.
Aperiodicity. Again, since , then for any proposed trade, so the probability of remaining in the same state is non-zero, and the chain is aperiodic.
Detailed Balance. For detailed balance, consider matrices and which differ by a curveball trade in rows and . The curveball algorithm involves permuting the elements of the trade candidates in and to create a proposed trade, captured in . The probability of obtaining is
where is the cardinality operator. Note that is the same for all such partitions, that is, all partitions are equally likely. Let , and for a generic index and index set s, define . Again starting from the left side of (10):
Hence the detailed balance equation holds, where the stationary distribution is again given by (1).
Appendix C Irreducibility With Monotonic Structural Zeros
When structural zeros are present, that is, at least one of , the Markov chain implied by checkerboard swapping may become reducible. Here we prove that in the special case of monotonic structural zeros, irreducibility still holds. Specifically we show that between any two matrices and , there exists a checkerboard swap which reduces the Hamming distance between them, so that repeated application of such swaps eventually reduces the distance to 0.
Finding a Swap. Consider , and let be the Hamming distance between them. The structural zeros are monotonic, so assume the rows and columns have been arranged such that for . This places all the structural zeros in the bottom-right corner (See Figure 9).
Let denote the first element where and differ, where “first” means no rows above contain a difference, and in row , no column to the right of column contains a difference. Specifically, satisfies:
for and ( and agree above row )