Vizualing motif contributions to sites

We can take the output from the model fit and visualize the motif contributions matrix omega using a stacked bar chart.

east_west_dir <- factor(site_metadata$WorE, levels = c("W", "E"))
elevation_metadata <- site_metadata$Elevation
ecos_blocks(fit$omega,
                    blocker_metadata = east_west_dir,
                    order_metadata = elevation_metadata)

Because we are interested in species turnover across elevational and geographic locations, we block the two elevational gradients (western or eastern) and plot them side by side.

This blocking structure (blocker_metadata) can be used for any environmental gradient of interest, allowing replicates or experimental treatments (in our case, geographic location) to be easily visualized together. order_metadata determines the order of the stacked bars in each block. Other parameters can be customized to modify the axes, titles, etc.

See ecos_plot_pie below to plot the data in a spatially explicit context (i.e., on a map).

Identifying species contributions to motifs

A central component of this framework is assessing the species contributions to each motif. One may be tempted to merely examine the theta matrix component of the model fit, and identify the the species that differentially drive the clustering.

We use the function ExtractTopFeatures from the R package CountClust to compute and record the top contributing species to each species abundance motif.

features <- CountClust::ExtractTopFeatures(fit$theta, top_features = 5,
                                           method = "poisson", options = "max")
t(apply(features$indices, c(1,2), function(x) return(rownames(fit$theta)[x])))

##      [,1]                       [,2]                         
## [1,] "Phylloscopus_chloronotus" "Phylloscopus_xanthoschistos"
## [2,] "Phylloscopus_reguloides"  "Cyornis_rubeculoides"       
## [3,] "Seicercus_whistleri"      "Macronous_gularis"          
## [4,] "Phylloscopus_occipitalis" "Orthotomus_sutorius"        
## [5,] "Aethopyga_nipalensis"     "Copsychus_malabaricus"

Comparing observed fit to null expectation

One may also want to compare the clustering of the original data to that of randomly generated null matrices of species abundances. We provide a way to compare the observed model fit to null model fits of various types. Randomized matrices of are generated internally using the randomizeMatrix function in the package picante.

out <- ecos_nullmodel(species_counts, K=2, null.model = "richness",
                              iter_randomized= 10, option = "BF")

## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 20927.8, 10.5, 19.1, 4.6, 17.5, 5.2, 0.3, 0.9, 3.2, 0.5, 0.5, 0.2, done.
## log BF( 2 ) = 1353.48
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 1845.7, 6.7, 3.1, 16.9, 7.2, 3.5, 3.4, 0.6, done.
## log BF( 2 ) = 1440.07
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 18432.1, 11.1, 21.8, 1.1, 9.4, 0.4, 0.5, 2.1, 2, 0.2, 1, 0.2, done.
## log BF( 2 ) = 1455.82
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 1707.1, 46.3, 8.3, 5.7, 0.4, 0.3, 1.5, 2.8, 1, 0.1, done.
## log BF( 2 ) = 1324.43
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 1820, 14.1, 4.9, 3.6, 1.9, 3.3, 0.4, 1.4, 0.7, 0.9, 2.7, 0.4, done.
## log BF( 2 ) = 1473.37
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 19254.9, 24.9, 9, 3.5, 0.8, 9.2, 2, 3.4, 0.4, 0.8, done.
## log BF( 2 ) = 1270.9
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 1844.4, 28.1, 9.4, 6.8, 1.1, 1.6, done.
## log BF( 2 ) = 1243.54
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 1842.1, 18, 6.3, 11.9, 7.4, 1.3, 1.4, 8.2, 1.3, 0.5, done.
## log BF( 2 ) = 1360.54
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 20732.8, 41, 18.9, 5.5, 7.1, 0.5, 0.2, done.
## log BF( 2 ) = 1386.47
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 19997.9, 14.6, 47.5, 4.4, 0.7, 0.6, 0.6, 16.3, 1.8, 1.1, 0, done.
## log BF( 2 ) = 1279.4
## 
## Estimating on a 38 samples collection.
## Fit and Bayes Factor Estimation for K = 2
## log posterior increase: 14985.7, 110.4, 43.3, 2.5, 0.1, done.
## log BF( 2 ) = 1761.84

Vizualing species motif contributions

We can take the output from the model fit and again visualize the membership proportions omega. However, as the data may consist of thousands of spatially explicit locations, we do not use ecos_blocks as in the previous example, but instead use the function ecos_plot_pie, which incorporates latitudinal and longitudinal coordinates of the original sites. We plot the omega matrix in a spatial context as follows.

ecos_plot_pie(omega = pres_ab_fit$omega,
                      lat_lim = c(5,50),
                      long_lim = c(55,120),
                      path = "geostructure_plot.png",
                      color= c("red","blue"))

## reading background map shapefile from inst/extdata/ne_110m_coastline 
##             folder

This function will create a file geostructure_plot.png in the current folder (getwd() in R to figure out).

One can easily apply this to any region of interest. We provide an example of visualizing species presence-absence structure as generated by applying the binomial GoM model with \(K=6\) on the bird species presence absence matrix from Australia. Again these data have been processed from BirdLife International

data("australia_birds")
australia_model <- ecos_fit(australia_birds$pr_ab_data, 
                            K=6, tol = 0.1)

The results from this fit, which takes < 2 minutes, are saved as australia_model and can be visualized as follows

data("australia_model")
ecos_plot_pie(omega = australia_model$omega,
              coords = australia_birds$latlong, 
              long_lim = c(110,160),
              lat_lim = c(-50,-10),
              path = "geostructure_plot_aus.png",
              color= c("orange", "red", "yellow", 
                       "deepskyblue", "chartreuse", 
                       "blue"))

The plot should look as follows

Generating dispersion fields for local sites

For a given dataset of local species presences (can be abundances but dispersion fields are generated based on presence), we can generate dispersion fields for each site and again fit a GoM model. These methods require some types of GIS data source, either shapefiles of species distributions, a geodatabase of species distributions (as in the .gdb of global bird species ditrbutions provided by BirdLife International), or even point occurences of species presences in a region.

We generate dispersion fields for our 38 sites using the function dsp_create_from_survey.

## If GIS data are in the form of separate shapefiles (.shp)

local_him_disp <- dsp_create_from_survey(local_data = species_counts,
                              gis_data_type = "shp",
                              shp_dir = {PATH TO directory with species shapefiles})

## If GIS data are in the form of a geodatabase (.gdb)

local_him_disp <- dsp_create_from_survey(local_data = species_counts,
                              gis_data_type = "gdb",
                              gdb_dir = {PATH To Directory with Species Geodatabase})

dsp_create_from_survey takes a local data matrix and a specified GIS data type to produce dispersion fields in a specified regions. The function reads in shapefiles for each species (or a .gdb for all species) and rasterizes those distributions for a given resolution. There are many arguments to customize so be sure to examine the documentation for this function. The output ‘local_him_disp’ consists of 3 objects.

• raster : A list of rasters representing the dispersion fields for each site (length is the number of rows (sites) in local data) at the specified resolution.

• precise : A list of rasters representing the dispersion fields for each site (length is the number of rows (sites) in local data) at the resolution of the argument precise. These could be used for plotting high resolution dispersion fields.

• matrix : A list of the same rasters is raster but in matrix form. These are used for easily processing the dispersion fields into the structure required to fit a GoM model.

Visualizing dispersion fields

We can generate maps for all of the disperison fields using the function dsp_plot_map. We provide the dispersion fields for our 38 Himalyan community surveys, dispersion_field_ex, as a list of matrices.

data("dispersion_field_ex")
dsp_plots <- dsp_plot_map(dispersion_field_ex, 
                    raster_latlim = c(5, 50), raster_longlim = c(50, 120), 
                    scale = "percentage")
dsp_plots[[1]]

dsp_plots is now a list of plots, each of which visualizes the dispersion fields from our local sites.

Processing dispersion fields for GoM model fit

We can easily process the dispersion fields into one data matrix by using dsp_to_matrix

dsp_fields_matrix <- dsp_to_matrix(dispersion_field_ex)

The output from dsp_to_matrix is a data matrix with rows equivalent to each site’s dispersion field as a vector.

dim(dsp_fields_matrix)

## [1]     38 201600

Each column is a map cell in the region of interest, with the value being a count from each site of how many species from that site overlap the cell.

GoM model and visualization of dispersion fields

As these data are counts, we can now fit a multinomial GoM (using to examine geographic motifs. In this case, omega is a matrix of motif contributions to dispersion fields, and the theta matrix represents each map cells contrbution to each of the motifs. We can plot the thetas on a map, showing the contribution of each map cell to the motif. Some biologists may think of these motifs as analagous to the regional inputs into local communities.

The model fit is done using ecos_fit as above. Because this geographic data matrix may be particularly large, depending on the resolution chosen to produce the dispersion fields, the fit may take longer than the fit on the original data matrix.

fit <- ecos_fit(dsp_fields_matrix, K = 2, tol = 0.1, num_trials = 10)

We can plot the theta distributions of each geographical motif using the function dsp_motif. We provide an example theta matrix with the package called example_theta.

data("example_theta")
himalayan_geo_motifs <- dsp_motif(example_theta, 
                                  color_ramp = c("black", "darkseagreen3",
                                                 "orange","red"),
                                  dsp_fld_res = 8)
himalayan_geo_motifs$motif_maps[[1]]

We call this color palette the Ghostbusters theme !! :)

Phylogeny

We load the time-calibrated phylogenetic tree phylo_tree of the 304 Himalayan birds in species_counts. We then collapse the species into clades to based on certain time point cut off on the tree collapse_at. Below we cut the tree off at 10 million years ago.

data("phylo_tree")
phylo_counts <- ecos_prepare_by_phylo(species_counts, phylo_tree, collapse_at = 10)
dim(phylo_counts$outdat)

## [1] 38 63

Traits

We use morphological trait metadata from himalayan_birds data and then use the distance metric based on these traits to collapse the bird species as above. As in the case of the time component above, one can set prop_div to a collapse the tips of trait dendrogram to produce a specific number of functional groups (the proportion of the original diversity). Below we create \(91\) functional groups out of the original 304 species.

bill_traits <- as.matrix(dist(scale(species_metadata[,c(1:3)])))
counts_bill_traits <- ecos_prepare_by_trait(counts = species_counts, 
                                                    traits = bill_traits, 
                                                    prop_div=0.3)
dim(counts_bill_traits)

## [1] 38 91

Both phylo_counts and counts_bill_traits can now be passed to ecos_fit to generate either phylogenetic motifs or trait motifs.