Goal

  • Find compositional method to measure component distribution shapes in Euclidean space

Init


In [3]:
%load_ext rpy2.ipython

In [4]:
%%R
library(coenocliner)
library(compositions)
library(dplyr)
library(tidyr)
library(ggplot2)


/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: This is coenocliner 0.2-0

  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: Loading required package: tensorA

  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: 
Attaching package: ‘tensorA’


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: The following object is masked from ‘package:base’:

    norm


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: Loading required package: robustbase

  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: Loading required package: energy

  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: Loading required package: bayesm

  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: Welcome to compositions, a package for compositional data analysis.
Find an intro with "? compositions"


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: 
Attaching package: ‘compositions’


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: The following objects are masked from ‘package:stats’:

    cor, cov, dist, var


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: The following objects are masked from ‘package:base’:

    %*%, scale, scale.default


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: 
Attaching package: ‘dplyr’


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: The following objects are masked from ‘package:stats’:

    filter, lag


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: 
Attaching package: ‘tidyr’


  res = super(Function, self).__call__(*new_args, **new_kwargs)
/opt/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.py:106: UserWarning: The following object is masked from ‘package:coenocliner’:

    expand


  res = super(Function, self).__call__(*new_args, **new_kwargs)

Simulate comunities


In [5]:
%%R
n.taxa = 500
n.fracs = 500

tolerance = 0.01


In [80]:
%%R 
set.seed(2)
M <- n.taxa                                 # number of species
ming <- 1.66                                # gradient minimum...
maxg <- 1.78                                # ...and maximum
meang = mean(c(ming, maxg))
locs <- seq(ming, maxg, length = n.fracs)   # gradient locations
opt <- rnorm(M, mean=meang-0.01, sd=0.01)        # species optima
tol  <- rep(0.01, M)                         # species tolerances
h    <- ceiling(rlnorm(M, meanlog = 12))    # max abundances
pars <- cbind(opt = opt, tol = tol, h = h)  # put in a matrix

In [81]:
%%R -w 600 -h 350
t.01 = coenocline(locs, responseModel = "gaussian", params = pars, countModel = "poisson")

matplot(locs, t.01, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Abundance", main='tol=0.01')


tollerance = 0.025


In [82]:
%%R
tol  <- rep(0.025, M)                         # species tolerances
pars <- cbind(opt = opt, tol = tol, h = h)  # put in a matrix

In [83]:
%%R -w 600 -h 350
t.025 = coenocline(locs, responseModel = "gaussian", params = pars, countModel = "poisson")

matplot(locs, t.025, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Abundance", main='tol=0.025')


tollerance = 0.05


In [84]:
%%R
tol  <- rep(0.05, M)                         # species tolerances
pars <- cbind(opt = opt, tol = tol, h = h)  # put in a matrix

In [85]:
%%R -w 600 -h 350
t.05 = coenocline(locs, responseModel = "gaussian", params = pars, countModel = "poisson")

matplot(locs, t.05, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Abundance", main='tol=0.05')


Fill in missing values


In [86]:
%%R
fill.val = 0.001
t.01[t.01 <= 0] = fill.val
t.025[t.025 <= 0] = fill.val
t.05[t.05 <= 0] = fill.val

Convert to relative abundances


In [87]:
%%R -w 600 -h 350

to.rel = function(x){
    y = sum(x)
    if(y == 0){
        return(x)
    } else {
        return(x / y)
    }
}

t.01.rel = t(apply(t.01, 1, to.rel))
t.025.rel = t(apply(t.025, 1, to.rel))
t.05.rel = t(apply(t.05, 1, to.rel))

matplot(locs, t.01.rel, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Relative abundance", main='tol=0.01')
matplot(locs, t.025.rel, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Relative abundance", main='tol=0.025')
matplot(locs, t.05.rel, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Relative abundance", main='tol=0.05')


CLR transformation


In [88]:
%%R -w 600 -h 350
t.01.rel.clr = clr(t.01.rel)
t.025.rel.clr = clr(t.025.rel)
t.05.rel.clr = clr(t.05.rel)

matplot(locs, t.01.rel.clr, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "CLR-trans abundance", main='tol=0.01')
matplot(locs, t.025.rel.clr, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "CLR-trans abundance", main='tol=0.025')
matplot(locs, t.05.rel.clr, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "CLR-trans abundance", main='tol=0.05')