Non-Uniform Sampling of Fixed Margin Uniform Matrices

by   Alex Fout, et al.
Colorado State University

Data sets in the form of binary matrices are ubiquitous across scientific domains, and researchers are often interested in identifying and quantifying noteworthy structure. One approach is to compare the observed data to that which might be obtained under a null model. Here we consider sampling from the space of binary matrices which satisfy a set of marginal row and column sums. Whereas existing sampling methods have focused on uniform sampling from this space, we introduce modified versions of two elementwise swapping algorithms which sample according to a non-uniform probability distribution defined by a weight matrix, which gives the relative probability of a one for each entry. We demonstrate that values of zero in the weight matrix, i.e. structural zeros, are generally problematic for swapping algorithms, except when they have special monotonic structure. We explore the properties of our algorithms through simulation studies, and illustrate the potential impact of employing a non-uniform null model using a classic bird habitation dataset.



There are no comments yet.


page 12

page 14

page 19


fastball: A fast algorithm to sample binary matrices with fixed marginals

Many applications require randomly sampling binary graphs with fixed deg...

By chance is not enough: Preserving relative density through non uniform sampling

Dealing with visualizations containing large data set is a challenging i...

Reinforcement Learning with Non-uniform State Representations for Adaptive Search

Efficient spatial exploration is a key aspect of search and rescue. In t...

Weighted matrix completion from non-random, non-uniform sampling patterns

We study the matrix completion problem when the observation pattern is d...

A Fast Scheme for the Uniform Sampling of Binary Matrices with Fixed Margins

Uniform sampling of binary matrix with fixed margins is an important and...

Phase transition in random contingency tables with non-uniform margins

For parameters n,δ,B, and C, let X=(X_kℓ) be the random uniform continge...

ClustGeo: an R package for hierarchical clustering with spatial constraints

In this paper, we propose a Ward-like hierarchical clustering algorithm ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In 1975, Diamond [8] 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 [7] 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 [25]. This accompanied theoretical work in binary matrices with fixed marginal sums [3]

and Markov Chain Monte Carlo (MCMC) sampling 

[2], 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 

[6], 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 [7] 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 [26] adapt the fill method of Connor & Simberloff to check whether the proposed fill will lead to a dead end.

Sanderson et al [25] 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 [5]

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 [21] 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 [15] 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 [7]. 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 [3] 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 [19] 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 [13] 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 [29], 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 [27] 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 [4] 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 [9] 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 [7] 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 [15], 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 [8] and Connor & Simberloff [7] 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 [15], 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

1 0
0 1

   0 1 1 0

Figure 1: A submatrix of from rows and , and columns and , which are not necessarily consecutive. Swapping the positions of 1’s and 0’s preserves the row and column sums of .

We propose a swapping algorithm similar to that of Brualdi [3], 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.

1 Let be the observed matrix;
2 for  in 1:N do
3        Set ;
4        Sample two elements and uniformly at random from the non-zero elements of ;
5        if  or  then
6               next  (no checkerboard);
7        end if
8       Sample ;
9        if  then
10               set and  (perform swap);
11        end if
13 end for
Algorithm 1 Sampling Via Weighted Checkerboard Swaps

If all weights are positive, then we have the following result:

Theorem 1.

Given binary matrix and weight matrix with , then Algorithm 1 generates a Markov chain with stationary distribution given by .

Proof of Theorem 1 is given in Appendix A. If weights are allowed to equal zero, the Markov chains may mix more slowly and even become reducible as we demonstrate in the next section.

2.3 Structural Zeros

As mentioned, structural zeros have received little attention since originally introduced by Connor & Simberloff [7]. 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 [3] 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 [23] 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:

Definition 2.1.

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

represent 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:

Theorem 2.

Given binary matrix and weight matrix where any structural zeros are monotonic, then Algorithm 1 generates a Markov chain with stationary distribution given by .

The only modification of the proof for Theorem 1 needed to prove Theorem 2 is to show irreducibility of the chain given structural zeros (see Appendix C).

2.4 Weighted Curveball Algorithm

Strona et al. [27] 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.

1 Let be the observed matrix;
2 for  in 1:N do
3        Set ;
4        Sample two rows and (such that ) uniformly at random from the non-zero elements of ;
5        Let ;
6        Let ;
7        if  or  then
8               next  (no trade);
9        end if
10        Let be the ordered list resulting from concatenating and ;
11        Permute the elements of uniformly at random ;
12        Let be the first elements of , and let be the last elements of ;
13        Let } (0 to 1 in row ) ;
14        Let } (0 to 1 in row ) ;
15        Compute ;
16        Sample ;
17        if  then
18               Set and for ;
19               Set and for ;
21        end if
23 end for
Algorithm 2 Sampling Via Weighted Curveball

Carstens [4] proves that the unweighted curveball algorithm can be used to sample uniformly, and we contribute the following theorem for the weighted case:

Theorem 3.

Given binary matrix and weight matrix where any structural zeros are monotonic, then Algorithm 2 generates a Markov chain with stationary distribution given by .

In Appendix B, we provide a proof of this theorem when there are no structural zeros, and in Appendix C this is extended to the case of monotonic structural zeros. R code implementing Algorithms 1 and 2 is provided on Github at (

3 Simulations

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.

1 0 0 0 1 0 0 0 1    0 1 0 1 0 0 0 0 1    1 0 0 0 0 1 0 1 0

0 0 1 0 1 0 1 0 0    0 1 0 0 0 1 1 0 0    0 0 1 1 0 0 0 1 0

Figure 2: All matrices with row and column sums equal to 1.

Table 1 shows agreement between the true and empirical probability distributions. The KL-Divergence between the empirical and true distributions is .

State Probability Empirical Probability
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)
Table 1: Theoretical vs. empirical probabilities for matrices in Figure 2 with binary matrix and weight matrix as defined in (6). Empirical probabilities also show Tukey-Hanning [10]

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 [17]. 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 [22].

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.

Figure 3: Average effective sample size for the simulations described in Section 3.2, for varying weighting schemes and percentage of structural zeros, where structural zeros are either randomly placed (blue), or arranged monotonically (red). Each line represents the average ESS from 10 separate chains using the same starting matrix.
Figure 4: Detailed view of effective sample size for the simulations described in Section 3.2, using a weight matrix of all 1’s and with 50% structural zeros. Average ESS is shown in bold, with each iteration shown partially transparent.

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.

(a) Initial matrix,
(b) Base weight matrix,
Figure 5: Initial binary matrix and base weight matrix for the simulation described in Section 3.3. Both use a shading scale where white denotes 0 and black denotes 1. The minimum value in is 0.5, achieved in the top-right and bottom-left most elements. Note the unfavorable state of , since it contains 1’s corresponding to lower weighted regions of .
Figure 6: Sampling distribution of the diagonal divergence statistic, under different weighting matrices, defined by raising the base weight matrix to the power .

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

is 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 [8], 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 [24]). Accordingly, Stone & Roberts [26] 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 [7], 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.

(a) Species-Island occurrence,
(b) Weight matrix,
Figure 7: Observed New Hebrides bird species combination  [8], and non-uniform weight matrix for the data analysis described in Section 4. For , white denotes 0 and grey denotes 1. For , light grey denotes 1 and grey denotes 15. With higher values indicating species (rows) preference for certain island (columns) habitats.
Figure 8: Sampling distribution of the C-score statistic, under uniform (red) and non-uniform (blue) weights. The empirical p-value under non-uniform weights is no longer significant, demonstrating how incorporating species habitat preference could impact whether the observed combination is relatively unlikely.

5 Discussion

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

a 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 [12]: “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 [23] 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 [11]. 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.

6 Acknowledgements

This work was partially supported by grant IOS-1856229 from the National Science Foundation.


  • [1]
  • Besag and Clifford [1989] Julian Besag and Peter Clifford. 1989. Generalized Monte Carlo significance tests. Biometrika 76, 4 (Dec. 1989), 633–642.
  • Brualdi [1980] Richard A. Brualdi. 1980.

    Matrices of zeros and ones with fixed row and column sum vectors.

    Linear algebra and its applications 33 (1980), 159–231.
  • Carstens [2015] 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.
  • Chen et al. [2005] 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. [2013] Edward F. Connor, Michael D. Collins, and Daniel Simberloff. 2013. The checkered history of checkerboard distributions. Ecology 94, 11 (Nov. 2013), 2403–2414.
  • Connor and Simberloff [1979] Edward F. Connor and Daniel Simberloff. 1979. The Assembly of Species Communities: Chance or Competition? Ecology 60, 6 (Dec. 1979), 1132.
  • Diamond [1975] Jared M. Diamond. 1975. Assembly of Species Comunities. In Ecology and Evolution of Communities. Harvard University Press, 342–443.
  • Diamond and Gilpin [1982] 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 [2010] 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.
  • Fosdick et al. [2018] 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. Publisher: Society for Industrial and Applied Mathematics.
  • Gotelli [2000] Nicholas J. Gotelli. 2000. Null Model Analysis of Species Co-Occurrence Patterns. Ecology 81, 9 (2000), 2606–2621.[2606:NMAOSC]2.0.CO;2
  • Gotelli and Entsminger [2001] 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.
  • Gotelli and Graves [1996] Nicholas J. Gotelli and Gary R. Graves. 1996. Null models in ecology. Smithsonian.
  • Harrison and Miller [2013] 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). arXiv: 1301.3928.
  • Harvey et al. [1983] 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.
  • Hoff [2009] Peter D. Hoff. 2009.

    A first course in Bayesian statistical methods

    Springer, London ; New York. OCLC: ocn310401109.
  • Jerrum and Sinclair [1996] 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. [1999] 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 [1995] Bryan F. J. Manly. 1995. A Note on the Analysis of Species Co-Occurrences. Ecology 76, 4 (1995), 1109–1115.
  • Miller and Harrison [2013] 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.
  • Plummer et al. [2006] 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.
  • Rao et al. [1996] 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 [1990] Alan Roberts and Lewis Stone. 1990. Island-sharing by archipelago species. Oecologia 83, 4 (1990), 560–567.
  • Sanderson et al. [1998] 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 [1990] Lewi Stone and Alan Roberts. 1990. The checkerboard score and species distributions. Oecologia 85, 1 (1990), 74–79.
  • Strona et al. [2014] 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.
  • Wilson [1987] J. Bastow Wilson. 1987. Methods for detecting non-randomness in species co-occurrences: a contribution. Oecologia 73, 4 (1987), 579–582.
  • Zaman and Simberloff [2002] 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 [3], 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

(A) at (-5-0.5, 5); (A2) at (); (B) at (0.5,5); (B2) at (); (size) at (5,-5); () rectangle (); () rectangle (); at () ; at () ; () – (); () – (); () – (); () – (); at () j; at () j; () – (); () – (); () – (); () – (); at () j’; at () j’; () – (); () – (); () – (); () – (); at () i; at () i; () – (); () – (); () – (); () – (); at () i’; at () i’; () – (); () – (); () – (); () – (); at () c; at () c; [fill=black] () – () – () – () – () – () – () – () – () – () – () – () – () – () – cycle; [fill=black] () – () – () – () – () – () – () – () – () – () – () – () – () – () – cycle; [pattern=north west lines] () – () – () – () – () – () – cycle; [pattern=north west lines] () – () – () – () – () – () – cycle; [fill=orange, line width=0] () – () – () – () – cycle; [fill=orange, line width=0] () – () – () – () – cycle; [color=red, line width=1.5] () – () – () – () – cycle; [color=blue, line width=1.5] () – () – () – () – cycle; [color=red, line width=1.5] () – () – () – () – cycle; [color=blue, line width=1.5] () – () – () – () – cycle; at () 0; at () 1; at () 1; at () 1; at () 0; at () 0;

Figure 9: Illustration of irreducibility proof in Appendix C. and have structural zeros (in black), and agree for elements in the hashed region. Their partials sums must agree in the red regions, and their partial sums must agree in the blue regions as well. Four cases are possible for the orange element.

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:

  1. for and ( and agree above row )

  2. for (