Spatially partition your randomisations

In previous posts I described how the rand_structured algorithm works, and then how one can use spatially structured randomisations to compare observed patterns with stricter null models.   

This post is concerned with how one can spatially partition a randomisation.  Some of these details are in a previous blog post, but this post is a much more detailed description. 

In the standard implementation of a randomisation, in Biodiverse and I suspect more generally also, labels are randomly allocated to any group across the entire data set.  This works well for many studies and is very effective.  However, when one begins to scale analyses to larger extents, for example North America, Asia or globally, then the total pool of labels begins to span many different environments.  The randomisations could allocate a polar taxon to the tropics, and a desert taxon to a rainforest.  This does not make the randomisation invalid, but it is perhaps not as strict as it could be.  

The effect is perhaps best considered with a scenario using phylogenetic diversity where polar and tropical taxa are distinct clades on the tree.   For a randomisation that allocates labels anywhere across the data set, one will commonly end up with a random realisation containing many groups with a mixture of polar and tropical taxa.  The PD on the random data set will be higher than on the observed data set in almost all cases because both clades are sampled and thus more of the tree is represented.  Note that this does not make the result wrong - it means that, when compared with a random sample from the full pool of taxa, the observed PD is less than expected.  That is useful to know, but the next question to ask is "are the patterns in the tropics higher or lower than expected for the tropics?", and the same for the polar taxa.  

One could quite easily remove any groups outside a region of interest and then rerun the analysis.  However, then the ranges of the labels will be reduced and any endemism analyses will not be comparable with the larger analyses.  

This can be readily fixed in Biodiverse by specifying a spatial condition to define subsets. Or, if you only want to randomise a subset of your data while holding the rest as-observed, you can use a definition query.    The spatial condition approach was used in Mishler et al. (2020) to (very approximately) partition Canada, the US and Mexico to assess sensitivity of the results for the full data set.  The definition query was used in Allen et al. (2019) to only randomise values within Florida, while holding the adjacent regions constant.

How does it work?  Both approaches slice the data into subsets, apply the randomisation algorithm to each subset independently, and then reassemble the randomised subsets into a full randomised basedata that is then used for the comparisons.  The difference for the definition query is that the data is divided into two sets, and only the subset that passes the query is randomised.

The definition query is run first.  This means that, if you specify both a definition query and a spatial condition, then only the groups that pass the definition query are considered for randomisations, and those randomisations will be spatially partitioned.

This approach is general, and works for any of the randomisation algorithms in Biodiverse.  Note, however, it will not apply different randomisation algorithms in different subsets.  If there is a good reason to do so then we can look at it, but it will make the interface much more complex to implement.  

One other point to note is that, for the structured randomisations, the swapping algorithm is applied within the subsets.  Each subset is run to completion before they are all aggregated.

Some images will work better than text, so here are some clipped screenshots from Biodiverse.

The observed distribution of Acacia tenuissima (red cells), with outlines of Australian states and territories in blue.  Acacia data are from Mishler et al. 2014, and the cell size is 50 km.  

A tenuissima incidences in a randomisation constrained such that each incidence is allocated within the polygon that contains it. Note the absence of incidences in NSW, Victoria and Tasmania, and that there are only three in South Australia.  This is consistent with the observed data.

A tenuissima records randomised using a definition query so that only incidences within Western Australia are randomised.  All others are kept unchanged. Compare with the observed distribution above.   

So how does one use it?  It is done as part of the standard interface (this functionality has actually been in Biodiverse since version 1, but the interface was slightly reconfigured in version 2).  The user specifies spatial conditions using the same syntax as for the spatial analyses.   

Subsets to randomise independently are defined using a spatial condition (red arrow), while the definition query (blue arrow) is used to randomise a subset of the data.

The choice of condition is up to the user, and they can specify whatever condition they like (it is their analysis, after all...).  Generally speaking, though, it is better to use a condition that generates non-overlapping regions.  One that uses a shapefile would fit with many geographic cases.  Many conditions will work but might not be very sensible, for example overlapping circles.  Each group is allocated to only one subset, so in such overlapping cases it is the first one that contains the group "wins" the group.

One thing to watch for, and which is fixed for version 4, is that if the spatial condition does not capture all groups then the system will throw an error.  The workaround for version 3 and earlier is to use a definition query that matches (shadows) the spatial condition, although keep in mind that any groups outside the definition query will not be randomised.

To sum up, Biodiverse allows a high degree of complexity and fine grained control when using randomisations to assess the significance of observed patterns.  This post has described the spatial conditions and definition queries.  More features and approaches are described in other posts, and are grouped under the randomisation tag.  

Shawn Laffan

23-Nov-2020

For more details about Biodiverse, see http://shawnlaffan.github.io/biodiverse/

To see what else Biodiverse has been used for, see https://github.com/shawnlaffan/biodiverse/wiki/PublicationsList

You can also join the Biodiverse-users mailing list at http://groups.google.com/group/Biodiverse-users

Randomisations - modelling spatial structure

In the previous post, I described how the rand_structured algorithm works.  If you don't feel like reading that post right now, then a short description is that it uses a filling algorithm to randomly allocate the labels in a basedata onto the groups of a new basedata, with the end result being random assemblages in each group, and where each group has the same number of labels as the original basedata by default. A swapping algorithm is then used to allocate any labels that could not be allocated using the filling process, usually because there were no unfilled groups that did not already contain these labels.

The filling algorithm has some nice advantages over other methods that also match the richness patterns, such as the independent swaps algorithm.  The main one is that filling algorithms are easily generalised to model spatial processes in their allocation, for example random walks and diffusion processes.

The videos below show some of these spatially structured randomisations at work, using Acacia tetragonophylla incidences from Mishler et al. (2014).  Note that the shade of red is used only to indicate when in the sequence a cell was allocated to, with lighter shades being earlier and darker later.

The first video shows a random walk model, in which a seed location is chosen for the first allocation.  The next location is a randomly selected neighbour of the current location, and the process repeats.  If a location has no available neighbours (all are full or already contain this label) then it backtracks until it finds one it can allocate to (backtracking can be from the last allocated, the first, or randomly selected).   If there are no possible allocations then it will try a new seed location.  One can also set a probability to randomly reseed even if the window has not be fully allocated to.   

The second uses a diffusion process.  This also uses a seed location, but considers the neighbours of all of its current locations for the next allocation.  As with the random walk, this process is also constrained by the filled and previously allocated groups, and will reseed if necessary or at random.   There is no backtracking as it is not needed since there is no "path" being followed. 

The final model is a proximity allocation process.  A seed location is chosen, and a spatial condition is then used to select a window of neighbours.   The example here uses a circle of a three cell radius, but any supported condition could be used such as predefined regions represented using a shapefile.  Once all groups in the window have been allocated to, excluding fully allocated groups, a new seed location is chosen. 

As with the main rand_structured algorithm, any unallocated incidences are handled using a swapping algorithm (see details in the previous post).  This ensures richness targets are met, at the cost of some loss of spatial contiguity in the randomly generated distributions.  

The outliers for this proximity allocation randomisation are due to the swapping algorithm used to reach the per-group richness targets.

The above randomisations are all implemented under the rand_spatially_structured randomisation, for which parameters can be set to model each approach.  However, setting parameters can be complicated and error prone so there are variants called rand_diffusion and rand_random_walk that restrict the set of arguments needed and set some useful defaults, but otherwise are just wrappers around rand_spatially_structured (which itself is implemented in Biodiverse as a particular case of the rand_structured algorithm). There is not currently a wrapper for the proximity allocation, but one could be added if needed.  The next few screenshots show the parameter settings.   

The options to control the spatially structured randomisations are marked with the red and blue arrows.  The spatial condition (marked by the red arrow) determines which neighbours of the considered group will be allocated to next.  This example uses a square, but it could be any of the supported conditions.  The other parameters (marked with the blue arrow) control reseeding, backtracking, and allocation order.  

The rand_diffusion parameters are almost identical to those for rand_spatially_structured, as the former is just a special case of the latter and this simplifies setting it up.   

As with the rand_diffusion parameters, and for the same reason, the parameters for rand_random_walk are also almost identical to those for rand_spatially_structured.   

As an aside, if you want a good introduction to the general process of dynamic spatial modelling, then have a look at O'Sullivan and Perry (2013).  There is also a good introduction to the limits of the complete spatial randomness model in O'Sullivan and Unwin (2010).  

To sum up, there is more to randomisations than just randomly shuffling the data around.  Biodiversity patterns are normally spatially structured, so finding patterns that are significantly different from those that are completely random is often not much of a challenge.  Adding constraints such as matching species richness patterns is very useful, but one can go further and introduce even more spatial structure, making the tests "harder to pass" (an example is in Laffan and Crisp, 2003).  As noted by Gotelli (2001), trying multiple hammers is a useful thing, and to quote that paper: "The advantage of null models is that they provide flexibility and specificity that cannot often be obtained with conventional statistical analyses".  Using multiple approaches might be regarded as a fishing expedition, but it does allow one to generate a distribution of models and thus deeper understanding of the patterns one is analysing.  And fish are nutritious.

Shawn Laffan

19-Nov-2020

For more details about Biodiverse, see http://shawnlaffan.github.io/biodiverse/

To see what else Biodiverse has been used for, see https://github.com/shawnlaffan/biodiverse/wiki/PublicationsList

You can also join the Biodiverse-users mailing list at http://groups.google.com/group/Biodiverse-users

Randomisations: Curveball algorithm now in Biodiverse

Biodiverse supports a range of randomisations to assess significance of analysis results.  Most use cases in the published literature use the rand_structured algorithm, which is explained in this post, but several common algorithms are supported.  

One of the design principles of Biodiverse is to give the user choice.  To that end, the curveball algorithm is available from version 5.  

The publication describing Curveball is Strona et al. (2014).  The name is derived from a baseball card trading card pastime popular in North America.  

The curveball algorithm is applied to a data set of items (species, genera, words, or some other set of identifiers).  In the common biodiversity case this is a sites by species matrix, transformed to a list of lists, e.g. a list of site lists, where each site list comprises its species (or vice versa).  These lists can be considered as sets.  At each iteration, two lists (sets of items) are randomly selected.  Any items found in both sets are ignored.  The rest can be swapped between the two sets, with the number swapped limited by the smaller number of unique items in the two sets to ensure after swapping that each set retains the same number of items it started with.  As an example, consider the case where set 1 has ten items, set 2 has eight, and there are six common items found in both lists.  This means two items can be swapped between the two lists.

The general formula for the number of possible swaps at an iteration is (min (|A|,|B|) - |A ∩ B|), where A and B are the two sets being considered, and the pipes || denote the lengths of the sets (the numbers of items they contain).   If one prefers to think in terms of dissimilarity measures where a is the number of shared items, b the number unique to set 1 and c the number unique to set 2, then the formula is (min (b,c)).  Purely as an aside, this is also part of the denominator in Simpson's dissimilarity index.  

The curveball algorithm is related to the independent swaps algorithm.  The chief advantage of curveball over independent swaps is that, because it swaps as many items as it can at each iteration, it converges on a randomised result much faster.  Curveball also avoids the main pitfall of the independent swaps algorithm where a pair can be selected that cannot be swapped, thus "wasting" an iteration (swap attempt).  

Curveball does, however, have the same issue that independent swaps has in that the user needs to specify the number of iterations over which swaps will be attempted.  Too few and the resulting matrix will not be sufficiently random.  Too many and time will be "wasted".  This is addressed in Biodiverse by optionally tracking which of the original matrix entries have been swapped, and stopping when all have been done (the stop_on_all_swapped parameter).  This has some overhead in the tracking but generally this should be balanced by the time saved by running fewer iterations overall.  For those interested, the default number of swaps is the same as for the independent swaps algorithm, which is twice the number of non-zero matrix entries (twice the sum of the lengths of all lists).

Accessing the curveball algorithm in Biodiverse is the same as for any of the randomisations.  Open the Randomisation tab, select rand_curveball as the randomise function, select the number of randomisation iterations and any other algorithm specific parameters, then press Go (see image below).  The results are in the same format as always (e.g. see here, here and here).

Since it is just another algorithm, all the common options are available (another new change in version 5 is that more options are available across all algorithms in the GUI - see issue 946).  Users can define regions that are randomised separately before reassembly for analysis, including some that are not to be randomised.  One can also add some of the randomised results to the project to inspect them.

In terms of speed, curveball is faster than rand_structured.  This is largely due to there being less book-keeping required.  However, as with independent swaps, curveball can only be applied on a per-cell basis.  It does not extend to spatially structured randomisations like rand_structured does (one could ensure swap candidates come from within some local neighbourhood, but this is a different model to something like a diffusion process or a random walk.  Update 20241109: This has been implemented and will be available in V5).

All that is needed to run the curveball algorithm is to choose rand_curveball as the "Randomise function".  Other parameters are set as usual.

And that's pretty much it for the description.  If you want to read more randomisation related blog posts then check out the posts tagged with the randomisation label.  

----

Shawn Laffan

05-Nov-2024

For more details about Biodiverse, see http://shawnlaffan.github.io/biodiverse/  

For a list of some of the analyses Biodiverse has been used for, see https://github.com/shawnlaffan/biodiverse/wiki/PublicationsList 

You can also join the Biodiverse-users mailing list at https://groups.google.com/group/Biodiverse-users or start a discussion at https://github.com/shawnlaffan/biodiverse/discussions 

Biodiverse now includes the independent swaps randomisation algorithm

Randomisations are one of the key analyses in Biodiverse.  They give an assessment of whether an observed result is significantly different from a distribution of randomly generated versions.   Several recent posts have described how the spatially structured randomisations work, including the (default) rand_structured algorithm and its spatially structured variants. These can also be spatially constrained.  

These approaches have been available in Biodiverse for some time: rand_structured since the very beginning, spatial constraints since version 1, and the spatially structured variants since version 2.  

One of the general principles when developing Biodiverse is to provide options for users[1].  In that light, version 4 of Biodiverse will include the independent swaps algorithm.  

How does independent swaps work?  

A short summary is below, and more details and assessments are in Gotelli (2000) and Miklos & Podani (2004).  Code for the R implementation is available through the Picante package.  

1. Pick two different groups at random (call them group1 and group2)
2. Pick two different labels at random (call them label1 and label2)
3. If label1 is not in group1, or label2 is not in group2, or if label1 is already in group2, or label2 is already in group1  
1. Then they cannot be swapped, so switch label1 and label2.  
2. If the switched pairs cannot be swapped using the condition in #3 then go back to step #1
4. Swap the labels between groups
5. Increment the number of successful swaps by 1
6. Start back at step #1

This is a very simple algorithm, and simple has many advantages when it comes to implementation.  However, it is also a brute-force approach.  The fact that the pairs are selected at random means there is no consideration of "swappability".  Consider the case of a site by species matrix where most of the matrix entries are zero, such as where most taxa are highly range restricted.  In that case, there is a very high chance that the algorithm will select a pair that cannot be swapped and must retry with a different selection, e.g. because neither label is in either group.  This can lead to many wasted CPU cycles.  As an unflattering description one might call this a "throw enough mud at a wall and some of it will stick" approach.  Keep in mind, though, that for many cases this will work well enough, and the results are most definitely random - good enough is good enough.  

Another issue with the independent swaps approach is that one can never be quite certain how many iterations are needed to fully randomise a data set.  Miklos & Podani (2004) suggest twice the number of non-zero entries in the matrix, but I am not aware of whether that has been rigorously assessed.  Too few iterations and some of the original structure will be retained.  Too many and the analysis might take too long (even for very patient people)[2].

There is also no stopping criterion except the total number of iterations.  It is therefore possible for the algorithm to be trapped in an infinite loop when none of the group/label pairs can be swapped, as the iteration count is only incremented on a successful swap.   

In light of the above points, the implementation of independent swaps in Biodiverse has three parameters, two of which lead to early stopping.  

1. The default number of swaps is set to be twice the number of non-zero cells (if one thinks of the data as a site by species matrix).  This value follows Miklos & Podani.  Users can specify any positive integer here.  The algorithm will stop when this value is reached.  
2. The maximum number of iterations defaults to 100 times the number of swaps, but can also be specified as a positive integer.  If this number is exceeded then the algorithm stops, regardless of the number of swaps completed.  (Note that this value will likely be increased in the final Biodiverse 4 release, see below for why).  
3. The number of times each label/group pair in the original matrix is swapped out is tracked.  The algorithm ends once each label/group pair has been swapped at least once.  This is off by default.

The last case uses the assumption that, once a pair has been randomly swapped, there is no need to swap it any further as it will not be any more random.  If your random number generator is not good then more swaps can actually make things worse (but his is unlikely to be a problem as most systems use good random number generators now). 

A modified independent swaps algorithm

As noted in a previous post, a huge amount of effort has gone into optimising the randomisations in Biodiverse to make them run faster.  Armed with this knowledge, the independent swaps algorithm has also been optimised.  The primary aim here has been to reduce the search for swappable label/group pairs, at the cost of extra memory to store extra index lists.  

An overview of the algorithm is this:

1. Select a group using a multinomial probability distribution, with the probabilities being the richness scores of each group.  Call this group1.
2. Select label1 from the labels in group1 using a uniform random distribution (each label in group1 is equiprobable).
3. Select a group that does not already contain label1, again using a multinomial probability distribution.  Call this group2.
4. Select label2 from the labels in group2 using a uniform random distribution (each label in group2 is equiprobable).
1. If needed, repeat #4 until label2 is not one that is already in group1.
5. Swap label1 to group2, and label2 to group1.
6. Update the multinomial probability distributions.
7. Increment the swap count.
8. Go to step #1.

The multinomial distributions are used to replicate the chances of selecting a non-zero matrix entry when selecting entirely at random.  The use of the "absence" lists in step #3 also focuses the search to swappable pairs.  For example, a very wide ranging label (species) will occur in most groups so there is a high probability that the normal algorithm will select a group1 and group2 that both contain that label.

There is clearly a memory burden when one considers the "absence" tracking lists.  However, memory is cheap these days and the data sets must be relatively enormous for it to make too much difference.  Perl, in which Biodiverse is programmed, also has a copy-on-write mechanism to save memory when dealing with large strings (e.g. a copy of a 100MB string can be stored using a link to the original, and only needs to be really copied when it is changed in some way or is written out).  

The modified implementation has the same default stopping parameters as the unmodified implementation.

How long do they take?

So how fast are they?  The table below summarises a benchmark run with the Acacia data set from Mishler et al. (2014).  This has 506 species (labels) distributed across 3037 cells (groups).   

Both the original and modified Biodiverse implementations are used, as is that from the R Picante library and a variant implemented in pure R (and run using R version 4.03).  All were run on the same computer.  The Biodiverse runs used the GUI, the R runs were under RStudio.  Code for the R runs is available at the Biodiverse GitHub repo.  

Only one to five runs of each combination was used, and values are rounded off.  Ideally more would be run, but the results are consistent when re-run and the differences are of sufficient magnitude to use for a blog post.

There are 29,175 non-zero entries in the site by species matrix, so approximately 81% of the matrix is empty.  The target number of swaps for the independent swaps runs was set to twice the number of non-empty matrix cells, 58,350, which as noted above is the default in Biodiverse and follows Miklos and Podani (2004).

Values in the table are sorted by run time in seconds.  The "early stop" column indicates if swapping is stopped once all label/group pairs (matrix cells) have been swapped at least once.  "Num swaps" is the number of swaps actually run, "cells swapped" is the count of matrix cells swapped at least once, and is only tracked when early stopping is in use.  "Swap attempts" is how many iterations were used overall, including those where a swap was not possible.  "max attempts" is the maximum number of swap attempts before giving up, with the default being 100 times the target number of swaps.  The larger number is 2^31-1, which is the maximum value for a signed 32 bit integer.    

The main takeaway from the table is that the rand_structured algorithm is fastest of all those run, and is much faster than any of the unmodified independent swaps runs.  This difference will be more pronounced as the data set size increases.   That said, the rand_structured algorithm has been through many rounds of optimisation, so the difference might be reduced with further optimisations to the other algorithms. 

The times for the modified independent swaps are aggregates of iterations 2-5, with the value in brackets being the first iteration when many of the lists are set up and cached for later reuse.  

The modified independent swaps algorithm reaches its target with much less "wastage" of iterations, hence it is fastest of the independent swaps.  This is further aided by the early stopping condition.  

For the unmodified independent swaps algorithm, more than 24,600,000  swap attempts are needed to ensure all matrix cells have been swapped.  This is ~843 times the number of non-zero matrix cells.  The two fastest unmodified runs reach the maximum number of swap attempts when the number of actual swaps is still very low, at ~6% of the target number of swaps.

The unbounded independent swaps run in Biodiverse takes nearly seven minutes.  If used with 999 randomisation iterations in sequence then the run time will be nearly 4.8 days for the randomisations alone (not including running the various analyses for the randomised basedatas).  A faster computer could probably reduce it to 3 days sequential time...

It should also be noted that the Biodiverse implementations are all in Perl.  The Picante implementation is written in C++, which is much faster than languages like Perl, R and Python (they themselves are written in C or C++).  

Despite not being written in C, the Biodiverse modified independent swaps algorithm is two to five times faster than the Picante implementation for this data set.  The modified algorithms do more work per iteration, but run fewer iterations overall.  The Perl implementation could also use be reimplemented in C, but the flexibility of the rand_structured approach (see below) means any such efforts will likely be directed there. 

As a final note regarding speeds, the pure R implementation could be made faster.  Profiling shows that most of the processing time is spent in the sample.int() function.  Even so, it would still not be faster than the Picante implementation so the only advantage would be if early stopping were implemented.

Spatially structured independent swaps

As a final note, an issue with the independent swaps algorithm is that it is very difficult to add spatial constraints to model diffusion processes or random walks.  One can apply spatial constraints to the swapping, i.e. an incidence must be swapped with another within some distance, but that has not been implemented.  If there is sufficient demand, and perhaps funding or a code contribution, then it could be added.

Keep in mind, though, that you can apply spatial constraints so any swapping is done within sub-datasets, for example within regions.  This process is described in the spatially constrained randomisation post.  

Footnotes

[1]  One of the general principles when developing Biodiverse is to provide options for users.  This can sometimes lead to a huge array of options, for example the indices.  However, focusing on that example, many of the indices are provided to allow insights into how others are calculated, such as showing the relative contribution of each label to an endemism score. 

[2] Speaking of patient, Biodiverse is written for and by the impatient - impatience is one of the three virtues of programming, after all.

Shawn Laffan

07-Dec-2020

For more details about Biodiverse, see http://shawnlaffan.github.io/biodiverse/

To see what else Biodiverse has been used for, see https://github.com/shawnlaffan/biodiverse/wiki/PublicationsList

You can also join the Biodiverse-users mailing list at http://groups.google.com/group/Biodiverse-users