This tutorial gives a basic introduction to the features built in to
library(sppairs). Before you go any further, please note that
sppairs is under active development, so the tools that are available – and the way that they work – may change slightly over time. However, I will update this page as those changes are made, so the examples given below should be reliable.
We begin by loading the required packages (using
devtools), along with some example data:
# prepare libraries library(devtools) install_github('mjwestgate/sppairs') library(sppairs) # then import a site by species matrix for example purposes library(cooccur) data(beetles) # remove rare species (<10% sites) dataset <-clean.dataset(beetles, cutoff.min=0.1)
Calculating pairwise associations
The first step we have to take is to choose a method for estimating pairwise associations between species (i.e. columns). Sppairs provides some options for this, e.g.:
or.symmetric(dataset[, 1:2]) # Symmetric odds ratio = 28.89 or.asymmetric(dataset[, 1:2]) # Asymmetric odds ratio = 11.45 mutual.information(dataset[, 1:2]) # MI = 0.235 cor(dataset[, 1], dataset[, 2]) # Pearson's correlation (from base R) = 0.62
Once we have chosen a method, we can apply it to all pairs of species in a dataset. For historical reasons, the default setting is to use asymmetric odds ratios as discussed by Lane et al (2014); but we don’t have to be limited to this.
# default is to calculate the symmetric odds ratio result <-spaa(dataset) # to calculate any asymmetric statistic, set asymmetric=TRUE result <-spaa(dataset, method='or.asymmetric', asymmetric=TRUE)
An unusual case is that given by Lane et al. (2014), where we have repeat visits to each site, and therefore wish to use mixed models to investigate pairwise associations:
# create grouping variable for example purposes ONLY groups <-as.factor(rep(c(1:10), length.out=dim(beetles))) # calculate odds ratio using mixed models result <-spaa(dataset, method='or.glmer', random.effect=groups, asymmetric=TRUE) # repeat the above, but give more information on the result result <-spaa(dataset, method='or.glmer', random.effect=groups, asymmetric=TRUE, complex=TRUE) # alternatively, return some raw coefficients from mixed models, rather than odds ratios result <-spaa(dataset, method='glmer.coef', random.effect=groups, asymmetric=TRUE)
method="or.glmer" will be much slower than other options, and will most likely give you a lot of error messages (passed from lme4); but these issues shouldn’t cause
spaa() to fail.
NOTE: As currently implemented,
spaa gives results that are marginally different from the example published by Lane et al. For example, comparing
spaa(method="or.glmer") to the +ve values in Table 3 of our paper in Ecology and Evolution using the same data shows a strong positive correlation, but not perfect concordance (lm gives: slope = 1.12 with SE= 0.06; intercept = 0.08 with SE= 0.32; adjusted r-squared = 0.87, n = 40). This is appears to be due to differences in calculation of logistic mixed models between GenStat (which we used for the paper) and R.
The code that we use to plot spaa results is for exploratory use only; if you want prettier diagrams, I’d recommend you use igraph (the more professional choice) or circleplot (my personal preference, but with fewer options).
# get some data result <-spaa(dataset, method='glmer.coef', random.effect=groups, asymmetric=TRUE)[, 1:3] # built-in plot code (very messy) plot.spaa(result) # plot using igraph graph.data <- graph_from_data_frame(result) plot(graph.data) # see ?igraph.plotting for more options # plot using circleplot result[, 3] <- log10(result[, 3]+1) # reduce effect of extreme values on the plot circleplot(result) # basic circleplot(result, plot.control=list(line.breaks=c(-2, -0.5, 0, 0.5, 2))) # prettier version