# Analysing patterns in species composition

```{r,warning=FALSE,message=FALSE,echo=FALSE}
#source("http://bioconductor.org/biocLite.R")
#biocLite("limma")
# library(devtools)
# install_github(qecol)
library(qecol)
library(vegan)
library(mgcv)

```

In many ecological settings multivariate analysis is concerned with looking at gradients of species composition or finding methods to classify sites based on their species composition. 

If you think about the methods that you have seen so far you may "reinvent" one of the earliest methods of analysing species composition. 

You could think about finding the correlation (or rank correlation) between the occurrence of one species and another.

If you do this for all possible pairs then you will obtain a matrix like this



```{r}


```


```{r}
data(BCI)
cmat<-cor(BCI,method="spearman")

```

This is data "overload" for a data set such as the BCI plots. Although you could possibly tease out some relationships from this for a few species that you know well, it is difficult to spot any pattern. A data reduction technique that can be applied to any correlation matrix is principal components analysis.  This could be used if the number of species were small in relation to the number of sites. The technique is often useful when dealing with multi-colinearity of environmental variables. However principal components analysis will not work in a meaningful sense with so many  variables.

Let's try PCA using only the most common species.

```{r}
selected<-names(sort(apply(BCI,2,sum),dec=T)[1:5])
selected
BCIcommon<-BCI[,names(BCI)%in%selected]

```


```{r,message=FALSE,warning=FALSE}

cor(BCIcommon)
Xpairs(BCIcommon)
```

```{r}
pca<-princomp(BCIcommon,cor=F,scores=T)
plot(pca)
summary(pca)

```

Let's take this output a step at a time. The scree plot shows how much of the variability is captured by each principal component. If we want to reduce the number of variables effectively we would like to see a steep scree plot. This would show a high degree of muktiple colinearity, which would in this sort of case be useful (the abundance of one species is a good predictor of the abundance of others).

The summary shows how much of the "explainable" variability is captured by each axes.  When we have as many axes as variables the cumulative variance explained will be 1. This does NOT mean that all the variability in the data is captured by five axes, simply that all the explainable variability has been.

If you have variables that are measured on different scales you should either transform them onto a common scale before running PCA or use the version of PCA that takes the correlation matrix. In this case I used cor=FALSE. PCA provides a set of loadings for each variable. If cor=FALSE the loadings correspond to a simple linear equation of the form

$score1=loading1*var1+loading2*var2 .. etc$

```{r}
loadings(pca)
```

```{r}
normalise<-function(x){
  (max(x)-x)/(max(x)-min(x))
}
normalise(scores(pca)[,1])
```

The data provided in the vegan package does not include any information
on measured local environmental variability (micro topography, soil, moisture etc). However the position of the plots are provided, as they were placed on a rectangular grid.  If we design a plot to show one of the PCA axes we might pick up a spatial pattern (based only on the five commonest species)


```{r}
x<- rep(seq(625754, 626654, by=100),each=5)
y<- rep(seq(1011569,  1011969, by=100),len=50) 
coords<-data.frame(x,y)
plot(coords,pch=21,bg=2,cex=3*normalise(scores(pca)[,1]))


```



```{r,eval=FALSE}
str(pca)
library(rgl)
spheres3d(pca$scores[,1:3],radius=5,col="red")
rgl.bbox()
```


## Differences between sites


Another way of condensing the sites by species matrix into a simpler form is based on an investigation of total compositional difference. This can be conducted on completely*a priori* grounds. A simple example is when species have been classified into taxonomic or functional groups. The proportional abundance in each group can then form a response variable. An example may be pioneer species and species of mature forest, or soft bodied invertebrates compared to hard bodied invertebrates.

However there is often a need to look at compositional differences that are derived from the data themselves. The starting point for many analyses involves calculating differences (or similarities) between sites. As we have seen, PCA does this, but it assumes linearity and is based on the same assumptions as simple linear regresion, which may not be appropriate in many cases. It is not a suitable method for data sets representing long ecological gradients with high species turnover.

The issue of how best to measure differences between sites is a potentially difficult one. It all depends upon what you mean be a difference.

Let's first take the simplest case in which only presence or absence
of species are considered.

Imagine just two sites. Site one has species A, B, C and D present. But site 2 has species A,B, E, F and G.

```{r,results='hide'}
library(limma)


Venn1<-function(vec1,vec2,names=c("Site 1","Site 2")){
  d1<-data.frame(id=sitea,abun1=1)
  d2<-data.frame(id=siteb,abun2=1)
	d<-merge(d1,d2,all=T)
	d[is.na(d)]<-0
	g <- cbind(d$abun1, d$abun2)
	vennDiagram(vennCounts(g),names=names)
	t(g)
}


Venn2<-function(g){
	vennDiagram(vennCounts(t(g)))
}

```

```{r}
#Example
sitea<-c("a","b","c","d")
siteb<-c("a","b","e","f","g")
Venn1(sitea,siteb)

```




One way of measuring the difference (or similarity) between the two sites would therfore be to look at the ratio of the number of shared species (2) against the total number of records made of species presences (4+5). This is known as Sorenson's index. 

$Sorenson=\frac{2c}{A+B}$

Where c= number of shared species, A = number of presence records for site A and B =number of presence records for site B. 

So the difference is 2*2/4+5

This can be extended to count data by simply weighting the values by the number of individuals. The index is then commonly known as the Bray Curtis index.

There are some important issues to be aware of when thinking about differences. What happens if site 1 has species A and B, while site 2 has species A,B,C,D, E and F?

```{r}
sitea<-c("a","b")
siteb<-c("a","b","c","d","e","f")
Venn1(sitea,siteb)
```

In one respect the sites could be considered to be identical, as the species in site 1 are just a nested sub set of those in site 2. In another sense they are quite different because species richness is clearly not the same. There are a range of other difference indices available in addition to Sorenson's. These weight such situations in different ways, and can result in quite different outcomes. You
should be aware of this issue if you go on to use these methods in a dissertation. 

You may need to look in more depth at other indices, depending on the questions you wish to adress. For our purposes we will use Sorenson's (Bray-Curtis) dissimilarity as an acceptable measure of differences between sites.

Assuming that we have decided on a way of measuring the difference between
two sites another question then arises. Why would we want to do this at all? If we have fifty sites there are 1225 paired comparisons between them. How can we possibly interpret the results?

A common reason for calculating differences between sites is to look for patterns of change in species composition. This can be achieved in many different ways. 

The BCI data appeared to come from a rather homogeneous area of tropical
forest, with high species diversity. So we might expect that any differences
between sites may be just the result of sampling issues. The plots
simply are not large enough to include all the species. Can we look
at this? 



## Mantel tests

One way to look for a pattern in species composition is to relate differences between sites in species composition to other differences between sites. The simplest measure of a difference between two sites is simply the distance between them. The so called first law of geography (that has been attributed to Tobler) is *Everything is related to everything else, but near things are more related than distant things.* This is generally true, but can we **test** this?

Well first let's extract the Bray Curtis difference between sites in the BCI data.

```{r}
compdist<-vegdist(BCI)
geodist<-dist(coords,"euclidean")
plot(compdist~geodist)
```

There seems to be a relationship, although there is a lot of noise.
However we can't simply use correlation or regression analysis to
test the significance of the difference because a difference matrix
has a lot of duplicated information. Instead we use a permutation
based procedure called a Mantel test.

```{r}
mantel(compdist,geodist)
```


There is a significant, and quite high, correlation between differences in composition and geographical distance.


## Correlation with environmental variables

The previous example did not use any environmental variables. It only allowed us to conclude that the sites showed some spatial dependence. Mapping the position of sites with similar composition may help to show the pattern. We will come on to that later. How could we use Mantel tests to look at the relationship with environmental variability? One issue is that if we want to use Euclidean distances the environmental variables should be on the same scale. The function bioenv cleverly finds the best subset of environmental variables, so that the Euclidean distances of scaled environmental variables have the maximum (rank) correlation with community dissimilarities. 

The example below uses the RIKZ data discussed by Zuur. You should look carefully at the data and Zuur's example analysis.

```{r}

RIKZsp<-read.csv("/home/aqm/data/Examples/RIKZsp.csv")
RIKZenv<-read.csv("/home/aqm/data/Examples/RIKZenv.csv")
bioenv(RIKZsp ~., RIKZenv)
```

This technique may be particularly useful for screening variables for further work.


## Ordination

One of the goals when analysing a sites by species matrix is to attempt to place sites in some sort of logical order based on their composition. The expectation is that the order of sites will reflect one or more environmental gradients. Therefore the techniques are used in direct and indirect gradient analysis. Indirect gradient analysis attempts to ordinate the sites first, and then explain the results in terms
of environmental gradients. Traditionally the techniques used for this include correspondence analysis and detrended correspondence analysis. These are forms of eigenanalyses that are conducted on the sites by species matrix. One of the issues with these technqies is the so called arch effect. This is an artefect that is removed by detrending. As detrending is a rather ``ad hoc'' procedure some influential authors and editors (particularly Jari Oksanen) recommend multi dimensional scaling instead. So we will look at this first.


## Non-metric multi dimensional scaling NMDS

Difference matrices are hard to understand because of the large number of dimensions. Multidimensional scaling tries to find one or more axes for the data that produce a configuration of the points in such a way that the distances between them are as close as possible to the actual distances in the distance matrix. For non-metric MDS, only the rank order of entries in the data matrix (not the actual dissimilarities) is assumed to contain the key information. So in this case, the distances of the final configuration should as far as possible be in the same rank order as the original data. 


## Single gradient NMDS

A straightforward way of using NMDS is to ask the algorithm to line up the species along one single axis. Some information is lost by doing this, and the technique may introduce some artefacts. However it is then easy to use the variable in other models (such as regression, or GAM's) as a simple measure of change in species composition.

```{r,warning=FALSE,message=FALSE}
library(mgcv)
nmds<-metaMDS(RIKZsp,k=1)
grad<-scores(nmds,"sites")
mod<-gam(grad~s(height,k=4)+s(exposure,k=3),data=RIKZenv)
anova(mod)
par(mfcol=c(1,2))
plot(mod)
```


The stress of an NMDS is a measure of how close to the true rank order the scaled values are. The lower the stress the better. There is no clear rule governing how much stress is tolerable, but Kruskall advised that ordinations with stress values of over 0.2 are likely to be unreliable. This ordination just falls just above this level.

It can sometimes be useful to look at a Shephard diagram for the ordination. This shows graphically how the stress arises by plotting the distance between the points in ordination space against observed dissimilarity.

```{r}
stressplot(nmds)
```



## Multiple axis NMDS

The more axes that are used in an NMDS the lower the overall stress. One of the reasons to use several axes is to look at the configuration of the points with respect to environmental variables and attempt to explain the resulting patterns. Understanding the results may require some further thought and work, as it is harder to see patterns in two or more dimensions.

```{r}
nmds<-metaMDS(RIKZsp,k=2)
stressplot(nmds)
```


One way to look at the pattern of dissimilarites in two dimensional space that is produced by NMDS is to plot contours representing the smoothed (using gam) mean values of environmental variables. For example we can plot the NMDS results with contours for height above sea level (blue) and exposure value (red).

```{r}
ordisurf(nmds,RIKZenv$height,col="blue")
ordisurf(nmds,RIKZenv$exposure,col="red",add=T)
```


A more traditional way of looking at this is to assume linearity and show the strength and direction of the linear relationship as vectors on the ordination plot. This may hide the shape of the relationship, but emphasises the size of the linear correlation.

```{r}
fit <- envfit(nmds~height+exposure, data=RIKZenv, perm = 999)
plot(nmds)
plot(fit)
```


## Clustering

Another way of working with the distance matrix is to attempt to group the data. Looking at the BCI data again.

```{r}
vdist<-vegdist(BCI)
clus<-hclust(vdist)
plot(clus)
```


There seems to be two main divisions. This may be useful in tracing the spatial pattern that we observed. We can cut the tree to find a grouping variable.

```{r}
grp<-cutree(clus,2)
plot(coords,pch=21,bg=grp,cex=2)
```


The results may be related to a slight ridge that appears along the middle of the plot. Or it may indicate disturbance of some kind.

![alt text](figs/googleBCI.png)




## Clam test

The clam function may help in finding which species are specialists within the groups of sites.

```{r}
clam<-clamtest(BCI, grp==1, alpha=0.005)
##clam
plot(clam)
```


In this case the division is not very convincing. However notice that the group contains comparatively high abundance of *Socratea exorrhiza*, the so called ``walking palm''. This may be indicative of recent disturbance. Thus the technique may have located a past forest opening due perhaps to logging. *Poulsenia armata* is another pioneer species found in this group.


## Simper

Another way of looking at which species are responsible for the grouping is to run a Simper analysis. This provides measures of the contribution of individual species to the overall Bray-Curtis dissimilarity between two sites. If the species are placed in order the analysis helps to find species that discriminate between the groups.

```{r}
simp<-simper(BCI, grp)
## summary(simp,ordered=T)
```


```{r}
head(summary(simp,ordered=T)[[1]],15)
```


## Clustering by species

Species can also be clustered by transposing the matrix. This should cluster species that co-occur in the same sites together. In the case of the large number of species found in the BCI data interpretation is a major task! It is only worth looking at the commonest species.

```{r}
BCI2<-BCI[,apply(BCI,2,sum)>100]
spdist<-vegdist(t(BCI2))
clus<-hclust(spdist)
plot(clus)
```



## Anosim

Although we may be able to find ways to group sites, we should use a statistical test of some kind in order to establish whether the differences between groups are significant.

If we can divide the data by groups (factors), such as sediment type, we can run an equivalent of a multiple analysis of variance (MANOVA) using the distance matrix as input. This technique is known as Anosim. It is particularly widely used by marine ecologists as a result of early implementation in the software Primer.

Anosim tests whether compositional dissimilarities between groups is greater than those within the groups. Plotting the results in R produces a notched boxplot that can be interpreted in the usual way (if the notches do not overlap then the median dissmilarity is significantly different).

We could use Anosim to test whether there really is a significant difference between the two groups of sites we obtained by cutting the cluster analysis of the BCI data.

```{r}
aa<-anosim(BCI,as.factor(grp))
aa
plot(aa)
```



## Adonis

Using Adonis is rather like fitting a multiple linear models to all the species in the matrix. It is also based on permutation tests and partitions the sums of squares of a multivariate data set in an analogous to MANOVA (multivariate analysis of variance). The inputs are linear predictors with the response being a matrix of more than two columns. 

Adonis tends to be very (too) powerful. In other words it will show significant relationships where other methods fail. This is because it only really needs one of the multiple species in the data set to have a strong relationship for the method to produce significance. Whether you can intepret the results easily is another matter.

```{r}
mod<-adonis(RIKZsp~.,data=RIKZenv)
mod
```



Asking for coefficents will produce a long list of the coefficents for each species (not shown).

```{r}
#coefficents(mod)
```


```{r}
library(rgdal)
dd<-data.frame(coords,grp)
coordinates(dd)<-~x+y

```


## Canonical Correspondence Analysis

Canonical correspondence analysis is a rather complex technique that combines ordination with regression in one algorithm. The technique ordinates sites and species along axes that are constrained to be linear combinations of environmental variables. For this reason the proximity of sites on a bi-plot may not reflect similarity in species composition as well as NMDS. 

A screeplot shows the eigenvalues for each axis. The total ``inertia''
o  the ordination is the sum of all the eigenvalues, not just those shown on the screeplot. In order to 

```{r}
cano<-cca(RIKZsp~.,data=RIKZenv)
screeplot(cano)
```


The technique can be used instead of Adonis in order to look at which species-environment relationships are most important. One way to do this is to treat the model rather like a multiple regression.

In R the AIC is calculated for the model. Thus backwards stepwise methods can be used to drop variables from the model. Burnham and Anderson type approaches could also be used for model selection, although there is no package in R that automates this.

```{r}
cano<-step(cano)
cano
```


After deciding on the model we can get a measure of the variance explained by looking at the proportion of the inertia that is attributed to the constrained axis. In this case it is 17.5\%.

Notice that because the axes are constrained to be **linear** combinations of environmental variables, they will tend to show clear orthogonality. In the case of the RIKZ data we have already seen that the relationship with height above sea level is not truly linear, thus the CCA has introduced some distortion.

```{r}
ordisurf(cano,RIKZenv$height,col="blue")
ordisurf(cano,RIKZenv$exposure,col="red",add=T)
```

The ordiplot function can be used to produce a dynamic biplot which allows labelling of species or sites. This avoids clutter.


```{r,eval=FALSE}
fig <- ordiplot(cano)
identify(fig, "spec") 
```


![""](figs/cca.png)



