multi_comp.Rmd

---
title: "The data we have: multivariate compositional data"
author: "gg"
date: '`r format(Sys.time(), "%d %B, %Y")`'
bibliography: ~/Library/texmf/bibtex/bib/bibdesk_refs.bib
fig_caption: true
output:
    pdf_document:
        fig_caption: yes
---

To run this file:
Rscript -e "rmarkdown::render('multi_comp.Rmd')"

#Part 1: relative abundance data

Whenever we analyze microbiome data, we analyze the data as `relative abundance' data, with or without normalization. What are relative abundance data?

Relative abundance data are data where the absolute count of something has no meaning. Instead, the only thing we can assess are the relationships between datapoints. That is, their ratios. This is compositional data.

When analyzing compositional data, or CoDa, (as we shall see) it is customary to log transform in some way and log-transformed data requires a different way of thinking. In linear space, a ratio is $x/y$, while in log space a ratio is $log(x) - log(y)$. So whenever we are working on log-ratios, what we are really saying is that we are working with the differences between log-transformed datapoints, and that these differences are equivalent to ratios in linear space.

A second issue is that of correlation of log-transformed data [@Lovell:2015]. In linear space, perfectly correlated data will follow the familiar formula $y = m \times x _ b$, where y and x are variables, m is the slope of the line and b is the intercept. However, in log transformed data, m becomes the slope of the line, and b becomes a non-linear parameter. Consider the simple case of plotting y as a function of x from 1:10, with b set to 0, or b set to 2.

Here we can see that in the case of the intercept being 0, the data are linearly related in both the normal and logarithmic spaces. In normal space, the intercept is 0, and the slope of the line changes. However, the two parts are always in constant ratio. In log space the intercept of the line changes, but the slope stays at 1.

In the case of the intercept not being 0, in the normal space the lines appear to be perfectly correlated, but in log space we can see that the lines are now curved because of the effect of a non-0 intercept. This curvature indicates that the ratios between x and y are changing, and so are not associated in logarithmic space.

This problem with correlation is not corrected by using non-parametric correlation measures! It is also worth noting that the recently proposed DESeq count normalization [@McMurdie:2014a] is simply a linear transform of the relative abundance values: i.e., if the count being normlized to is $1x10^6$, then this is simply relative abundance $x 1x10^4$.


```{r, echo=FALSE, results='show', fig.width=8, fig.height=8, message=FALSE,fig.cap="Plots to show the effect of log-transformation on simple linear lines. The grey dashed lines in the log space plots are lines of constant ratio."}
y <- function(m,x,b){m * x + b}

par(mfrow=c(2,2))
x <- seq(1:10)
b <- 0
plot(x, y(1,x,b), col="red", pch=19, type="b", xlim=c(0,10), ylim=c(0,20), main="normal space")
points(x, y(2,x,b), col="blue", pch=19, type="b")
points(x, y(0.5,x,b), col="orange", pch=19, type="b")
text(1.5,15, labels="y=2x", col="blue")
text(1.5,12.5, labels="y=x", col="red")
text(1.5,10, labels="y=0.5x", col="orange")

plot(x, y(1,x,b), col="red", pch=19, type="b", xlim=c(1,10), ylim=c(0.5,20), log="xy", main="log space")
points(x, y(2,x,b), col="blue", pch=19, type="b")
points(x, y(0.5,x,b), col="orange", pch=19, type="b")
text(1.5,18, labels="y=2x", col="blue")
text(1.5,12, labels="y=x", col="red")
text(1.5,8, labels="y=0.5x", col="orange")
abline(0.1,1, lty=2, col=rgb(0,0,0,0.2))

b <- 2
plot(x, y(1,x,b), col="red", pch=19, type="b", xlim=c(0,10), ylim=c(0.5,20), main="normal space")
points(x, y(2,x,b), col="blue", pch=19, type="b")
points(x, y(0.5,x,b), col="orange", pch=19, type="b")
text(1.5,15, labels="y=2x + 2", col="blue")
text(1.5,12.5, labels="y=x + 2", col="red")
text(1.5,10, labels="y=0.5x + 2", col="orange")

plot(x, y(1,x,b), col="red", pch=19, type="b", xlim=c(1,10), ylim=c(2,20), log="xy", main="log space")
points(x, y(2,x,b), col="blue", pch=19, type="b")
points(x, y(0.5,x,b), col="orange", pch=19, type="b")
text(1.5,18, labels="y=2x + 2", col="blue")
text(1.5,12, labels="y=x + 2", col="red")
text(1.5,8, labels="y=0.5x + 2", col="orange")
abline(0.5,0.6, lty=2, col=rgb(0,0,0,0.2))

```

\newpage


#PART 2: multivariate compositional data

Univariate data have is one variable per sample. This is the typical kind of data generated by biologists. For example, you measure the height of people in a room, divide them into male and female and compare. In this case, we have one variable (height) measured for two groups (male and female). These data are unconstrained, within the bounds of human height, and are in general what traditional statistical tools were developed for.

Multivariate data have more than one variable per sample. Multivariate data can be independent or dependent. So for example, if the variables are truly independent (randomly chosen, not linked) then each may be treated as univariate. For example if we measure height, hair color and handedness of people, then each person is an observation that holds three variables. We can assume that these are (relatively) independent, and all statistical tests univariate and multivariate should be valid and importantly, multiple test corrections would be valid.

However, if the variables are dependent, then we have many unappreciated problems. Unfortunately, the typical high throughput sequencing dataset is multivariate and highly dependent [@gloor2016s;@Gloor:2016cjm]. This dependency is forced upon the data by the sequencing instrument itself.

We acknowledge this dependence when we call the data `relative abundance' data. What does this really mean?

Let us set up a thought experiment. Let us imagine we have a very simple dataset composed of 100 samples and three taxa. Note that everything about this example generalizes to datasets with more samples and  more taxa, it is just difficult to show this with an n-dimensional graphic. In this example, we will generate counts for 100 samples from a single experiment where the taxa abundance can range anywhere between 1 and an arbitrary number of counts. We will assume that we are actually counting the number of molecules belonging to each taxon, and that there is no practical upper limit on the number of molecules.

We will save the special case of 0 for our practical demonstration. Now we have 100 multivariate observations with three variables.

```{r, echo=FALSE, results='show', fig.width=8, fig.height=3, message=FALSE,fig.cap="Counts were generated randomly for three taxa labeled a, b and c, and were free to range from 1 to 1000 for a, up to 2000 for b, and up to 3000 for c. One hundred samples were generated. The correlation coefficient of each sample vs. each other is shown above the plots."}
library(compositions)
library(scatterplot3d)

a <- round(runif(100, min=1, max=1000))
b <- round(runif(100, min=1, max=2000))
c <- round(runif(100, min=1, max=3000))
abc <- rbind(a,b,c)
tot=1
abc.comp <- apply(abc, 2, function(x){x/sum(x)} * tot)
abc.acomp <- acomp(abc, tot=1)

par(mfrow=c(1,3))
plot(a,b, main=round(cor(a,b),3))
plot(a,c, main=round(cor(a,c),3))
plot(b,c, main=round(cor(b,c),3))
```

The plots show that the variables in each of the samples are randomly placed. Therefore, we can see that the data in each pair of samples are essentially uncorrelated, as we expect for randomly generated data. This is what we would expect for randomly generated data where each point is absolutely independent of each other: in other words this is the best case scenario that would be rarely seen in a biological context. We shall now constrain the data to a constant sum and see how this affects the shape of the data.

```{r, echo=FALSE, results='show', fig.width=8, fig.height=3, message=FALSE,fig.cap="The same data as above were constrained to sum to an arbitrary sum: in this case, the data were converted to proportions. Any arbitrary number (percentage, ppm) has the same effect. The correlation coefficient of each sample vs. each other is shown above the plots. We can see very clearly that the associations between each of the datapoints are now not independent."}
par(mfrow=c(1,3))
plot(abc.comp[1,], abc.comp[2,], main=round(cor(abc.comp[1,],abc.comp[2,]),3))
plot(abc.comp[1,], abc.comp[3,], main=round(cor(abc.comp[1,],abc.comp[3,]),3))
plot(abc.comp[2,], abc.comp[3,], main=round(cor(abc.comp[2,],abc.comp[3,]),3))
```

The simple act of converting each count to a proportion markedly skews the data. Here we can see that the data now appear much more correlated and constrained. While this is a simple example with three taxa, it is generalizable to any number of samples, and to any number of variables (taxa) in a sample.

Why is this so?

```{r, echo=FALSE, results='show', fig.width=8, fig.height=3, message=FALSE,fig.cap="Plotting these data in three dimensions shows the reason for the non-independence. The unconstrained data are on the left, the constrained data are in the middle, and the right side shows the data as if we were looking directly down onto the plane of data seen in the middle box. Data points are coloured from red to black to show their place on the b axis."}
par(mfrow=c(1,3))
scatterplot3d(t(abc), scale=1.5, highlight.3d=T)
scatterplot3d(t(abc.comp), xlim=c(0,tot), ylim=c(0,tot), zlim=c(0,tot), scale=1.5, highlight.3d=T)
plot(t(abc.acomp), pch=19, cex=0.5, col=rgb(0,0,0,0.5))

```

The essential problem is one of geometry. When the data are unconstrained, as shown in the box on the left, the data points are scattered uniformly within the box. This is a visualization of multivariate data that are unconstrained and uncorrelated. Essentially, in data of this type, knowing information about one datapoint gives no information about any other datapoint.

When the data are constrained to a constant sum as shown in the middle box, then the datapoints collapse to a flat plane within the box, with limits at the corners. This plane is called a simplex, and occurs because if we know information about all points but one, then we know the information about the last point as well. For example, if the data must sum to 1, then knowing that the previous 10 points have a sum of 0.9, we know that the last point must have a value of 0.1. Thus, we can see intuitively that the data are constrained. We can re-orient our view to look directly perpendicular to the simplex plane and observe how the three dimensions of the data map onto the two dimensional space on the right.

The data are clearly in a very different geometry than are the unconstrained data points. Moving from place to place in the box on the left is an additive process since the difference between points is a linear distance. In the simplex, moving from point to point is a multiplicative process, and we think of differential taxon abundance in terms of being multiples of the original abundance.

In order to analyze such data, we must place ourselves on the simplex, or modify the data to move it from a constrained to an unconstrained simplex. Thus, we need to recast our thinking into a compositional data analysis (CoDa) way of thinking.

#References and Further reading