06-A-unconstrained-ordination.Rmd

---
layout: topic
title: "Unconstrained ordination"
author: Tad & Anna
output: html_document
---

**Assigned Reading:**

> Chapter 5 from: Borcard, D., Gillet, F. and Legendre, P. 2011. *Numerical Ecology with R.* Springer. [link](https://link.springer.com/book/10.1007/978-1-4419-7976-6)
>
> Ramette, A. 2007. Multivariate analyses in microbial ecology. *FEMS Microbiology Ecology* **62**: 142-160. [DOI: 10.1111/j.1574-6941.2007.00375.x](http:dx.doi.org/10.1111/j.1574-6941.2007.00375.x)


```{r include = FALSE}
# This code block sets up the r session when the page is rendered to html
# include = FALSE means that it will not be included in the html document

# Write every code block to the html document 
knitr::opts_chunk$set(echo = TRUE)

# Write the results of every code block to the html document 
knitr::opts_chunk$set(eval = TRUE)

# Define the directory where images generated by knit will be saved
knitr::opts_chunk$set(fig.path = "images/06-A/")

# Set the web address where R will look for files from this repository
# Do not change this address
repo_url <- "https://raw.githubusercontent.com/fukamilab/BIO202/master/"
```

### Key Points

#### Unconstrained vs. constrained ordination

Fig 3 in Ramette 2007 - left and right sides: Borcard et al. Ch 5 and 6, respectively

#### PCA

+ Scaling 1 -- use this if main interest is objects

"distance biplot" (distance among objects is Euclidean distance, angle among descriptors meaningless) - Fig 5.2

+ Scaling 2 -- use this if main interest is descriptors

"correlation biplot" (angles between descriptors reflect their correlations, distance among objects is not approxmate of Euclidean distance) - Fig 5.2

+ `vegan`'s default is scaling 2  (a bit strange?)

+ Combining PCA with clustering can be powerful for confirmation (Fig 5.3)

#### CA

+ Chi-square distance not influenced by double zeros, unlike Euclidian distance in PCA; data have to be frequency or freq-like

+ Scaling 1 vs. 2 – similar to PCA

+ post-hoc `envfit()` (Is constrained ordination better? Useful with NMDS?)

#### PCoA

+ Distance measured by any similarity measure

+ Q mode (objects, or sites) and R mode (variables, or species or env factors)

+ Rarely used (Ramette's Table 1) - used more now?

#### NMDS

+ Any distance matrix (commonly, Bray-Curtis dissimilarity)

+ May be sensitive to initial configuration (local optimum) - Could use PCoA ordination as a start. `metaMDS()` uses random starts

+ Needs fast computer – now not much of a problem

+ "A test winner, and a natural choice ..." (Jari Oksanen)

+ Bray-Curtis dissimilarity - Reasons why popular (Clarke and Warwick 2001):

1. 0 when identical
2. 1 when no species shared
3. Not affected by measurement unit
4. Not affected by inclusion or exclusion of species jointly absent
5. Not affected by inclusion of exclusion of a third sample
6. Differences in both total and relative abundances considered

+ Shepard diagram (Fig 5.11)

+ Stress (Clarke and Warwick 2001)

1. stress < 0.05: excellent representation
2. stress < 0.1: good representation
3. stress < 0.2: acceptable representation
4. stress > 0.3: unsatisfactory representation

#### Which method to use? (Criteria by Ramette 2007)

+ Linear (PCA) vs. unimodal (CA) - Use plot of species abundances along ordination axes as diagnosis

+ Gradient length short (PCA) vs. long (CA), as determined by DCA

+ Absolute abundance (PCA) vs. relative abundance (CA)

+ Many zeros (CA)

+ Linear (PCA) vs. non linear (NMDS, PCoA)

+ Shepard diagram to decide between NMDS and PCoA

+ Discrete groups of objects expected (cluster) vs. not expected (ordination)

#### Further reading: 

> "Vegan: an introduction to ordination" by Jari Oksanen 
> [link](https://cran.r-project.org/web/packages/vegan/vignettes/intro-vegan.pdf)


***

#### Unconstrained vs. constrained ordination

Fig 3 in Ramette 2007 - left and right: Borcard et al. Ch 5 and 6, respectively

#### PCA

+ Scaling 1 -- use this if main interest is objects

"distance biplot" (distance among objects is Euclidean distance, angle among descriptors meaningless) - Fig 5.2

+ Scaling 2 -- use this if main interest is descriptors

"correlation biplot" (angles between descriptors reflect their correlations, distance among objects is not approxmate of Euclidean distance) - Fig 5.2

+ `vegan`'s default is scaling 2  (a bit strange?)

+ Combining PCA with clustering can be powerful for confirmation (Fig 5.3)

#### CA

+ Chi-square distance not influenced by double zeros, unlike Euclidian distance in PCA; data have to be frequency or freq-like

+ Scaling 1 vs. 2 -- similar to PCA

+ post-hoc `envfit()` (Is constrained ordination better? Useful with NMDS?)

#### PCoA

+ Distance measured by any similarity measure

+ Q mode (objects, or sites) and R mode (variables, or species or env factors)

+ Rarely used (Ramette's Table 1) - used more now?

#### NMDS

+ Any distance matrix (commonly, Bray-Curtis dissimilarity)

+ May be sensitive to initial configuration (local optimum) - Could use PCoA ordination as a start. `metaMDS()` uses random starts

+ Needs fast computer -- now not much of a problem

+ "A test winner, and a natural choice ..." (Jari Oksanen)

+ Bray-Curtis dissimilarity - Reasons why popular (Clarke and Warwick 2001):

1. 0 when identical
2. 1 when no species shared
3. Not affected by measurement unit
4. Not affected by inclusion or exclusion of species jointly absent
5. Not affected by inclusion of exclusion of a third sample
6. Differences in both total and relative abundances considered

+ Shepard diagram (Fig 5.11) and stress

+ According to Clarke and Warwick 2001:

1. stress < 0.05: excellent representation
2. stress < 0.1: good representation
3. stress < 0.2: acceptable representation
4. stress > 0.3: unsatisfactory representation

#### Which method to use? (Criteria by Ramette 2007)

+ Linear (PCA) vs. unimodal (CA) - Use plot of species abundances along ordination axes as diagnosis

+ Gradient length short (PCA) vs. long (CA), as determined by DCA

+ Absolute abundance (PCA) vs. relative abundance (CA)

+ Many zeros (CA)

+ Linear (PCA) vs. non linear (NMDS, PCoA)

+ Shepard diagram to decide between NMDS and PCoA

+ Discrete groups of objects expected (cluster) vs. not expected (ordination)

#### Further reading: 

> "Vegan: an introduction to ordination" by Jari Oksanen 
> [link](https://cran.r-project.org/web/packages/vegan/vignettes/intro-vegan.pdf)



***


### Analysis Example

Unconstrained ordinations are used as exploratory methods to understand multidimensional data. Generally they use eigenvalues to represent new synthetic axes that explain the most variation in the data/cluster of samples (and are orthogonal to one another). Ordination can be thought of as a projection of multidimensional scatterplot onto new axes, ideally a few of which explain significant variation.

The most common unconstrained ordination methods are:

1. **PCA** (Principle component analysis) Often uses a correlation matrix of standardized response variables. Assumes linear relationships between variables and input data structure.  PCA is for shorter environmental gradients.

2. **PCoA** (Principle coordinate analysis) uses a distance matrix as input. Shows distance or dissimilarities between objects in ordination space.

3. **CA** (Correspondence anlysis) Maximizes the correspondence between sample and species variation scores. Better for longer gradients than PCA, assumes unimodal distribution of points along axes.

4. **NMDS** (Nonmetric multidimensional scaling) ranked order not original distances and based on distance measure.


```{r libraries, message=FALSE, warning=FALSE, include=FALSE}
library(dplyr)
library(tidyr)
```


```{r, message=FALSE, warning=FALSE, include=FALSE}
# Make a new phyloseq object
library(phyloseq)
otumat2 <- read.csv("~/Documents/Stanford/Eco_stats_202/data/OTUsnormalized.csv", sep= ",") #update with proper path
rownames(otumat2) <- paste(otumat2$X) #make OTU numbers the row names
otumat2 <- otumat2%>% 
  select(-X) #remove OTU column
OTU <- otu_table(otumat2, taxa_are_rows = TRUE)


system.order<-c("Fresh", "Oligohaline", "Mesohaline", "Polyhaline", "Euhaline")
stat.order<-as.character(c(657, 649, 3, 6, 9, 13, 15, 18, 21, 24, 27, 30, 34, 36))
envmat <- read.csv("~/Documents/Stanford/Eco_stats_202/data/envdat_phy.csv", header=T, row.names = 1) #update with proper path
envmat$Station.Number<-as.factor(envmat$Station.Number)
envmat$system<-factor(envmat$system, levels=system.order)
ENV <- sample_data(envmat)

TREE<-read_tree("~/Documents/Stanford/Eco_stats_202/data/otu.tre") #update with proper path

physeq2 <- phyloseq(OTU, ENV, TREE)
physeq2
```

#### Prepare the data

Center and scale environmental data (or variables, could be traits of individuals, etc)

Important notes for `scale` function:

 + center = TRUE means column means are subtracted
 + scale = TRUE means centered columns are divided by standard deviations
 
```{r, echo=TRUE, message=FALSE, warning=FALSE}
envdat.phy <- as.data.frame(sample_data(physeq2))
envdat.phy.scaled <- scale(envdat.phy[,-c(13:16)], # cols 13:16 are factors
                           center=TRUE, # center data by subtracting column means 
                           scale=TRUE) # divide centered data by standard deviation
```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
str(envdat.phy.scaled)
```

#### PCA

General steps:

 + Scale and center data (some programs can do for you in the ordination step)
 + Check eigenvalues (scree plots)
 + Look at ordination plots and appropriate axes

##### Vegan package
```{r, echo=TRUE, message=FALSE, warning=FALSE}
library(vegan)
library(gclus)
library(ape)

env.pca <- rda(envdat.phy[,-c(13:16)], scale=TRUE)    # rda function makes PCAs
                                                    # scale means divide by species unit variance
                                                    
env.pca
```

```{r, echo=FALSE, message=FALSE, warning=FALSE}
summary(env.pca, scaling=1)
```

Look at the eigenvalues. Numerical Ecology in R recommends keeping all axes that are above the mean. One can also use the "elbow" method, look for a strong decrease in eigenvalues followed by slower decreases like a bend in the elbow. If several PCs have similar variances they are ill defined. If you include one you should include all of them as arbitrarily choosing the first one may not be the correct interpretation.


```{r}
ev <- env.pca$CA$eig
ev>mean(ev)
```

Divide eigenvalues by total inertia to explain the proportion of variance explained. Inertia should be the sum of all the eigenvalues (total variance of data). The below plots show eigenvalues with the mean plotted as a horizonatl line(top) and the broken stick model (bottom) which caluclates the eigenvalues you would expect to get if you picked PCs randomly.

```{r, message=FALSE, warning=FALSE, include=FALSE}
n<- length(ev)
bsm <- data.frame(j=seq(1:n), p=0)
bsm$p[1] <- 1/n
for(i in 2:n){
  bsm$p[i] = bsm$p[i-1] + (1/(n+1-i))
}
bsm$p <- 100*bsm$p/n
```


```{r, echo=TRUE, message=FALSE, warning=FALSE}
par(mfrow=c(2,1))
barplot(ev, main="eigenvalues", col="bisque", las=2)
abline(h=mean(ev), col="red")
legend("topright", "Average eignenvalue", lwd=1, col=2)
barplot(t(cbind(100*ev/sum(ev), bsm$p[n:1])), beside=TRUE, main="% variance", col=c("bisque", 2), las=2)
```

You can look at projected points in two scalings. Scaling 1  is good for euclidean distances between objects (samples/sites). Scaling 2 is for understanding the correlation (angles between arrows) between species (variables).

```{r, echo=TRUE, message=FALSE, warning=FALSE}
par(mfrow=c(1,2))
biplot(env.pca, scaling=1, main="PCA scaling 1") # scaling 1 for distances between objects
biplot(env.pca, main="PCA scaling 2") # correlation biplot, for seeing correlation between response variables (species) see angles. 
```


##### Alternative PCA functions

There are lots of functions for creating PCAs in R including `princomp` (stats) and `dudi.pca` (ade4).

Notes:

+ Eigenvalue = measures the importance of variance of an axis. 
+ Inertia = "variation"; sum of variance of variables (covariance matrix) or sum of diagonal values (correlation matrix) 

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

# Just to prove that most should be doing the same things, look at single values for different methods (coefficients in front of variables to make new variables contain more information).

library(ade4)

svd(scale(envdat.phy[,-c(13:16)]))$v[,1] #single value decomposition matrix

prcomp(scale(envdat.phy[,-c(13:16)]))$rotation[,1] # principle components analysis 
#rotation is matrix of variable loadings (columns contain eigenvectors), uses svd on x to get singlular values

princomp(scale(envdat.phy[,-c(13:16)]))$loadings[,1] # uses eigen on correlation/covariance matrix and determined by cor

dudi.pca(envdat.phy[,-c(13:16)], nf=2, scannf=F)$c1[,1] # center is centering by mean, scale uniform row weights, scannf is to display screeplot, nf is the number of kept axes c1 is the column normed scores

# the signs of the columns of loadings and scores are arbitrary
```


```{r, echo=TRUE, message=FALSE, warning=FALSE}
library(ade4)
pr.env2 <- prcomp(envdat.phy.scaled)
prin.env2 <- princomp(envdat.phy.scaled)
dudi.env2 <-dudi.pca(envdat.phy.scaled, center=F, scale=F, scannf=FALSE, nf=4) # ade4 library
```

```{r}
pr.env2
```


```{r, echo=TRUE, message=FALSE, warning=FALSE}
prin.env2 # sv is the standard deviations of PCs
          # loadings is matrix of variable loadings, columns contain eigenvectors
          # Only for looking at scaling 2??
```

```{r, echo=TRUE}
dudi.env2
```

These should all be doing the same thing, so when we look at the eigenvalues they should be the same. And it looks like the first 4 axes are important.
```{r, echo=TRUE, message=FALSE, warning=FALSE}
meaneig <-mean(dudi.env2$eig) # .98

library(factoextra)
eig0 <- fviz_eig(pr.env2, geom="bar", width=.4, title="prcomp")+theme(panel.grid=element_blank())
eig1 <- fviz_eig(prin.env2, geom="bar", width=.4, title="princomp")+theme(panel.grid=element_blank())
eig2 <- fviz_eig(dudi.env2, geom="bar", width=.4, title="dudi.pca")+theme(panel.grid=element_blank())

library(gridExtra)
grid.arrange(eig0,eig1, eig2, nrow=1)
```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
p0<-fviz_pca_ind(pr.env2, geom="point",
             habillage= envdat.phy$system)

p0
```


```{r, echo=TRUE, message=FALSE, warning=FALSE}
fviz_pca_ind(pr.env2, geom="point",
             habillage= envdat.phy$system,
             addEllipses = TRUE, 
             ellipse.level=.69)
```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
prin.p <- fviz_pca_biplot(pr.env2, 
                label="var", 
                habillage=envdat.phy$system, 
                axes = 1:2)

prin.p
```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
fviz_pca_var(pr.env2, col.circle="black")+
  xlim(c(-1.2,1.2)) +
  ylim(c(-1.2,1.2))
```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
p1 <- fviz_pca_ind(pr.env2, geom="point",
             habillage= envdat.phy$system,
             axes= 2:3)

p2 <- fviz_pca_ind(pr.env2, geom="point",
             habillage= envdat.phy$system,
             axes= 3:4)

p3 <- fviz_pca_ind(pr.env2, geom="point",
             habillage= envdat.phy$system,
             axes= c(1,3))

p4 <- fviz_pca_ind(pr.env2, geom="point",
             habillage= envdat.phy$system,
             axes= c(1,4))

p5 <- fviz_pca_ind(pr.env2, geom="point",
             habillage= envdat.phy$system,
             axes= c(2,4))

grid.arrange(p0,p1,p2,p3,p4,p5)
```

#### CA

Should only be done on frequency data but just to show that you can use the wrong ordination if you want. The output doesn't make a whole lot of sense.
```{r}
dud.env.coa<- dudi.coa(envdat.phy[,-c(13:16)], 
                       nf=2,
                       scannf=FALSE)

```

```{r}
fviz_eig(dud.env.coa)
```

Joint plot shows samples/sites and variables/species. Species are near center of inertia of samples where they are most likely to occur

```{r}
fviz_ca_biplot(dud.env.coa, 
               geom=c("point", "text"),
               label="col")
```


#### PCoA

PCoA is common in microbial ecology. Here are some examples of PCoA on OTU counts using the phyloseq package (based on vegan).

*bray-curtis dissimilarity*

```{r, echo=FALSE}
library(phyloseq)
ord.pco3 <- ordinate(physeq2, method="PCoA", distance="bray")
p.b <- plot_ordination(physeq2, 
                ord.pco3, 
                color="Salinity", 
                shape="season")+
  geom_point(size=3)+
  scale_color_continuous(low="blue", high="goldenrod2")+
  labs(main="PCoA Bray")+
  theme(panel.grid=element_blank())

#ord.pco3
```

```{r}
plot_scree(ord.pco3)
```

```{r, echo=FALSE}
p.b
```

*jaccard distance*

```{r, message=FALSE, include=FALSE}
ord.pco1 <- ordinate(physeq2, method="PCoA", distance="jaccard")
p.j <-plot_ordination(physeq2, 
                ord.pco1, 
                color="Salinity", 
                shape="season")+
  geom_point(size=3)+
  scale_color_continuous(low="blue", high="goldenrod2")+
  labs(main="PCoA jaccard")+
  theme(panel.grid=element_blank())

plot_scree(ord.pco1)
```

```{r, echo=FALSE}
p.j
```


*wunifrac*

If you update to ape version 5.0, download the newest phyloseq version. Should be available on joey711 github if the newest bioconductor version still does not have the fix.

Unifrac is based on shared branch lengths between samples and calculated as (sum of unshared branch lengths)/(sum of all tree branch lengths (shared+unshared)). Weighted unifrac (wunifrac in phyloseq) incorporates the abundance of taxa as well.

```{r, eval=TRUE, message=FALSE, warning=FALSE, include=FALSE}

ord.pco2 <- ordinate(physeq2, method="PCoA", distance="wunifrac")
p.w <- plot_ordination(physeq2, 
                ord.pco2, 
                color="Salinity", 
                shape="season")+
  geom_point(size=3)+
  scale_color_continuous(low="blue", high="goldenrod2")+
  labs(main="PCoA weighted unifrac")+
  theme(panel.grid=element_blank())

plot_scree(ord.pco2)
```

```{r, eval=TRUE, include=FALSE}

p.w
```

A note on horseshoes and arch effects: if there is a lot of species turnover and samples at each end of a gradient have no overlapping species, then placing their true distance from each other is difficult. Since my gradient spans freshwater to saltwater, there is basically no overlap on the OTU level between the two ends of the gradient so a horshoe or arch is to be expected.

#### CA on taxa data

```{r}
ord.cca <- ordinate(physeq2, method="CCA")
ord.cca
```

```{r, echo=TRUE}
plot_scree(ord.cca)

```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
p.cca <-plot_ordination(physeq2, 
                ord.cca, 
                color="Salinity", 
                shape="season")+
  geom_point(size=3)+
  scale_color_continuous(low="blue", high="goldenrod2")+
  labs(main="CCA")+
  theme(panel.grid=element_blank())

p.cca
```

#### NMDS

While the original distance between samples is not preserved, NMDS is useful for understanding ranked distances in a 2D space.

And choice of distance matters a lot!

*jaccard*

```{r, echo=TRUE, message=FALSE, warning=FALSE}
ord.nmds1 <- ordinate(physeq2, method="NMDS", distance="jaccard")
ord.nmds1
```
```{r, echo=TRUE, message=FALSE, warning=FALSE}
plot_ordination(physeq2, 
                ord.nmds1, 
                color="Salinity", 
                shape="season")+
  geom_point(size=3)+
  scale_color_continuous(low="blue", high="goldenrod2")+
  scale_shape_manual(values=c(7,6,16,8))+
  labs(main="stress .10")+
  theme(panel.grid=element_blank())
```

*weighted unifrac*

```{r, eval=TRUE, message=FALSE, warning=FALSE, include=FALSE}

ord.nmds2 <- ordinate(physeq2, method="NMDS", distance="wunifrac")
ord.nmds2
```

```{r, eval=TRUE, message=FALSE, warning=FALSE, include=FALSE}
plot_ordination(physeq2, 
                ord.nmds2, 
                color="Salinity", 
                shape="season")+
  geom_point(size=3)+
  scale_color_continuous(low="blue", high="goldenrod2")+
  scale_shape_manual(values=c(7,6,16,8))+
  labs(main="stress= .07")+
  theme(panel.grid=element_blank())
```

#### Indirect Gradient Analysis

Environmental fit function `envfit` from vegan can be used to get explanatory variables. Looks at correlation between ordination and environmental variables. Arrow points in the "directin of the gradient" and length reflects strength of correlation.

```{r, echo=TRUE, message=FALSE, warning=FALSE}
ef.nmds <- envfit(ord.nmds2, envdat.phy.scaled, permutations = 999)
ef.nmds
```

```{r, echo=TRUE, message=FALSE, warning=FALSE}
plot(ord.nmds2, display = "sites", cex=1)
plot(ef.nmds, p.max = 0.1)
```


### Discussion Questions

+ Do you use unconstrained ordinations? 
+ What kind of data are unconstrained ordinations useful for?
+ Are they intuitive or do you need to thoroughly know linear algebra?