Matt Michalska-Smith
Zachary Miller
Stefano Allesina

This talk was originally prepared for a departmental seminar entitled Darwin's Weekly, at the University of Chicago Department of Ecology & Evolution on 2 February 2018. It was modified from this to be presented at the Annual Meeting of the Ecological Society of America on 6 August 2018. Text here will supplement slide material.

You can remove this pane (and get it back) by tapping v on your keyboard.

Some of the slides will have animations which require advancing the slide partway through these explanations. These will be indicated by "next slide" in these notes.

Slides can be advanced by tapping to display only the core slides, or Spc, PgDn, or n to also show some supplemental slides throughout the presentation.

Finally, this presentation will look best in full-screen, but should be more or less intelligible on smaller screens (no guarantees about phone screens). So, whenever you're ready...

next slide

Stefano Allesina

Zachary Miller

A motivating example
What is known
What we've done
What remains

This project stemmed from a statistical observation that arises when one is dealing with data strongly influenced by an underlying, one-dimensional gradient. Think, for example, of species distributions on the side of a mountain. For each species, there will be an elevation with the greatest abundance, and some probability of observing a member of the species outside of that elevation which decreases with distance.

Thus, if we were to sample species presence/absence data along this elvational gradient, we might collect a presence/absence matrix that looks something like the one here on the right. Here, each row is a site, for which we have recorded either the presence (red cell) or absence (white cell) of each species found anywhere on the mountain.

?

A motivating example
What is known
What we've done
What remains

Wartenberg, Ferson & Rohlf 1987

Ecologists have been interested in ordination for a long time and fairly early on, Principal Component Analysis (PCA) was identified as a powerful, agnostic tool in this effort. For instance, this figure, shows six (well, I'm not convinced by the bottom row, but at least four) examples of the arch pattern across a range of disciplines from Psychology (A) to Ecology (B) to Astronomy (C) to Archaeology (D). For a period of time, the arch observed in PCA of one dimensional data was viewed as a statistical artifact, distorting what was known to be an underlying linear gradient, but due in part to an increasing understanding of the mathematics uderlying the arch, efforts, starting with this publication by Wartenberg, Ferson, and Rohlf in 1987, have shifted away from "correcting" the ordination coordinates. Nevertheless, the presence of the arch still presents challenges for extracting the optimal ordering of the data.

Original Legend from Wartenberg Ferson, and Rohlf 1987:
FIG. 2.-A few examples of curvature of component scores from many fields of study. A, Psychology, Nonmetric, multidimensional scaling ordination of "judged similarities among 14 approximately spectral colors" (after Shepard 1978). The numbers correspond to wave- lengths. Subjects rated their perceived pairwise similarity between all possible pairs of 14 objects each with a different color. The derived arch shows a progression from long- wavelength colors (reds) to short-wavelength colors (blues). B, Ecology, Principal- components-analysis ordination of pitfalls based on species of lycosid spiders (after van der Aart 1973). The distribution of pitfalls along the observed arch reflects a sequence of biotopes with vegetation changing from bare sand and lichens at one end through grassy areas, meadows, parks with shrubs, and, finally, relatively dense woods. The symbols correspond to classes of pitfalls. C, Astronomy, Nonmetric, multidimensional scaling ordination of galaxy types based on Mahalonobis distance using morphological data such as size, color, brightness, etc. (after Nathanson 1971). The letters correspond to types of galaxies. The arch reflects a gradual transition in the physical properties of each galaxy. D, Archaeology, Nonmetric, multidimensional scaling ordination of graves in La Tene cemetery at Munsingen- Rain based on the presence or absence of artifacts found within (after Kendall 1971). The numbers represent Hodson's ordering of the graves. E, Anthropology, Procrustes ordination of humanoid populations based on skull characteristics (after Gower and Digby 1981). Groupings of points represent different populations. F, Sociology, Principal-components- analysis ordination of U.S. cities based on the frequency of different types of crimes (after Gower and Digby 1981).

For our purposes, we can take a very simple distance metric, where we first normalize the adjacency matrix and then multiply it by its transpose. This produces a matrix where high values indicate that the row and column species have many presences and absences in common. When we do this with our toy data from earlier, we get this new "projection" matrix. Interestingly, if species distributions are laid out as I suggested, then nearby sites should have similar species compositions, and if we were to re-order the species into the right configuration (e.g. by order of the elevational center of their distribution)

next slide

a pattern appears in both the occurance matrix and this new distance matrix. Note that there is some artifactually higher values in the lower left and upper right corners. This is due to an interaction between our simulating framework and the normalization we do on the occurance matrices prior to multiplying them. Basically, these species look more similar than they ought because they share many absences, even though their presences are nearly all unshared.

When looking at large datasets, we see this distance matrix approaching what is called a "Toeplitz" structure in mathematics. When the gradient loops back upon itself,

next slide

the matrix instead converges toward a structure known as "Circulant" and produces a corresponding circle in PC space.

The final piece of history was quite recently brought into the spotlight in Ecology in a paper by Novembre & Stephens. It turns out that the reason for this arch actually stems from the fact that PCA is an eigenvalue decomposition. An efficient means for compressing information is to sum a sequence of sinusoidal curves. This is the method used for compressing images into the .jpeg format, or audio files into .mp3's.

next slide

Likewise, the PCA tries to compress our occurance data and ends up using the exact same process, leading to what is termed the Discrete Cosine Transform for the bounded case (i.e. for the Toeplitz Matrix), and the Discrete Fourier Transform in the looping (i.e. Circulant) one. These are both sinusoidal decompositions and provide a nice intuition for why the PC plot produces the characteristic arch in these datasets.

A motivating example
What is known
What we've done
What remains

First, I'm going to lay out the simulation framework used to create all of the datasets seen so far in the talk. We can simulate a gradient quite simply as a line, where one can think of one end as the low end and the other as the high end of the gradient.

next slide

Along this gradient, we will sprinkle some sampling sites according to a prescribed distribution. For this talk, I am always using a uniform distribution, but this is easy to adjust.

next slide

Next, we sprinkle some species. These lines indicate what will be the centers of each species' distribution, which, in the case of normally distributed species, is the location of their most likely observance. Again, I am distributing these centers uniformly in this talk for convenience of presentation, but changing this and other assumptions in this simulation framework is a key area of further research.

next slide

Finally, each of these species will have some distribution around that center. For this figure and throughout the talk, I am just using a simple normal distribution. Within the framework, it is quite easy to change both the distribution used by each species, as well as to make it such that different species have different distribution kernels.

next slide

And, of course, these simulated datasets have more than 10 sites and 5 species, so there are actually very many more sites and species distributions at a time...

next slide

leading to occurance matrices like our example from earlier in the presentation. Importantly, because we know the locations of each site and the centers of each species' distribution, we know how to order this matrix to make the gradient apparent. When testing our algorithm, we will randomize this matrix prior to searching in order to test our ability to recover this known ordering.

From this occurance matrix, we can normalize and multiply by it's transpose to get the distance matrix, as mentioned before, and the PCA and eigenvectors of this matrix can likewise be plotted. Looking at this properly ordered case, it is apparent that there are multiple potential avenues of approach if we wanted to go backwards to extract this ordering.

First, one could look at the PCA plot iself, since this is order invariant, and try to trace along the arch to get the proper ordering. Though this looks easy in this case, this is a large simulated matrix with low variance leading to a very clean arch that is usually not the case in real data.

Likewise, one could try rearranging the eigenvector components to best fit the known theoretical transform. Again, we would run into some issues with noisy data, but it is an option.

Finally, one could work with the distance matrix, or, equivalently the occurance matrix, directly, rearranging rows or columns in order to put the most similar rows/columns next to one another. This is the path we've chosen as, in the early stages of prototyping and testing, it was shown to be both faster and more consistent at generating favorable results compared to fitting the eigenvectors.


# take a binary occurance matrix B
# for a decreasing neighborhood size d
for (d in ncol(B) / 2, ncol(B) / 4, ncol(B) / 8, ..., 2, 1) {
  # loop through each column to get the contribution to the total score
  for (each column_i in B) {
    for (each column_j in B) {
      # measure the euclidean distance between the two columns
      distance_i[j] <- sqrt(sum((column_i - column_j)^2))
    }
    # rank the distance to each other column
    rank_distance_i <- rank(distance_i)
    # add the ranks of the columns within d steps of the target column
    for (each column_j within d steps of column_i) {
      total_score <- total_score + sum(rank_distance_i[j])
    }
  }
}

We utilized a Metropolis-Coupled Markov Chain Monte Carlo (MCMCMC, or MC$^3$) stochastic search algorithm. The core functionality of this algorithm is the fitness function, which is what this pseudocode outlines. Basically, a good fit is one in which, for each column, the other columns that are most similar to the target column/row are also the nearest in position. To quantify this, we rank the distances to every other column from a target column and then sum those ranks for a local neighborhood around that target column, finally summing those values for every column/row in the matrix. Thus, the lower this value, the more low-rank distance columns/rows are nearby the target column/row and the better the fit. To improve the quality of the final fit, we narrow in the neighborhood of interest, starting with a broad neighborhood and getting down to just considering the immediate neighbors of each column/row.

So, how does this work? We can test it with our simulated data. This is how the results look when they are first generated, and ...

next slide

this is how the best solution found from a randomized starting configuration looks. If you can't see the difference, that is a good thing -- it means we were able to recover near the optimal ordering. This can be further demonstrated if we plot the two orderings against one another.

next slide

As you can see, there are a few swaps, but in general, this is a correlation of greater than 0.95. Also note the fitness scores shown in the axes labels -- we actually found a better fit (lower distance value) than the original. This could be due to a mismatch between our fitness measure and the optimal ordering, or it could be that, due to the stochastic nature of the data generation, this new order is actually more parsimonious than the original.

So, the algorithm works, at least with highly controlled, simulated data. What about in the real world? This map shows a subset of the experimental sites taking part in the Nutrient Network -- a replicated experiment which contains plots of land for which plant percent cover data has been collected for over a decade. If we build an occurance matrix analogous to thos shown earlier in the talk, i.e. species x plot, and perform the usual analysis of buiding a distance matrix and running a PCA on it, we see our familiar arch. Even more exciting, if we color the sites by longitude,

next slide

We see a clear pattern along the arch, which can be accentuated by a path drawn between the mean positions for each longitude:

next slide

The one extreme outlier in this case, the light blue in the upper left, is the site in Tyson, MO. Remarkably, this site is about to be closed down because it is suffering a rapid invasion of species more common in the Eastern US.

So how does our algorithm fare with this relatively nice, though clearly more noisy, empirical case? Here is the before picture, with the plots clustered by site, but unordered otherwise.

next slide

And here is the best found ordering of the plots. The change here is harder to interpret. the eigenvetors are smoother, but certainly not obviously related to the expected discrete cosine transformation. Likewise, the distance matrix looks a little closer to Toeplitz structure, but still not especially close.

next slide

If we again plot the new order agains the old, we tend to see clusters of points rearranged in the new order. These are whole sites being moved, rather than individual plots being picket out of sites. The next steps will be to see if this best-found order correlates with any of the known metainformation we have for the sites in question, for instance longitude.

A motivating example
What is known
What we've done
What remains

Identifying rank vs. identifying position
Gradients in more than one dimension
Types of variability

Species center distribution
Site location distribution
Species distribution parameters
Species distribution class

It is iportant to be clear about what this algorithm is actually doing; it is not identifying the position along the gradient of each species/site. Rather, it is identifying a rank-order for the species/sites along the gradient. These are different problems, but the former would be the more powerful result.

Second, are there similar patterns in the principal components or eigenvectors when a more complicated gradient underlies the data, such as a two or more dimensional gradient? This is something that has been under-explored, but which is quite easy to investigate using the framework explained in this talk.

A third area of expansion is more clearly defining the limits to our inference. Explicitly, how does introducing variation at different levels of the simulation affect our ability to consistently recover the proper orderings?

Questions?

Slides available on Michalska-Smith.com

Stefano Allesina

Zachary Miller

Outline

A motivating example

A motivating example

A motivating example

?

?

Outline

What is known

Wartenberg, Ferson & Rohlf 1987

What is known

What is known

What is known

What is known

Outline

What we've done

What we've done

What we've done

What we've done

What we've done

What we've done

Outline

What remains

Questions?

Gradients in higher dimensions

Increasing species distribution variance

Alternative Distance Metrics

Nestedness as Gradient