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)

In [5]:

    

%%R
n.taxa = 500
n.fracs = 500


    

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')


    

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')


    

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')


    

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


    

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')


    

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')


    

In [89]:

    

%%R

t.01.rel %>% apply(1, var) %>% sum %>% print
t.025.rel %>% apply(1, var) %>% sum %>% print
t.05.rel %>% apply(1, var) %>% sum %>% print


    

    






[1] 0.09266403
[1] 0.004524045
[1] 0.003185228

In [91]:

    

%%R -w 600 -h 350
t.01.rel.var = t.01.rel %>% apply(1, var)
t.025.rel.var = t.025.rel %>% apply(1, var)
t.05.rel.var = t.05.rel %>% apply(1, var)

df.var = data.frame(t.01 = t.01.rel.var,
                    t.025 = t.025.rel.var,
                    t.05 = t.05.rel.var)

matplot(locs, df.var, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Variance")
matplot(locs, df.var[,c(2,3)], lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Variance")


    

In [68]:

    

%%R
t.01.rel.clr %>% apply(1, var) %>% sum %>% print
t.025.rel.clr %>% apply(1, var) %>% sum %>% print
t.05.rel.clr %>% apply(1, var) %>% sum %>% print


    

    






[1] 3490.541
[1] 973.4967
[1] 523.3454

In [93]:

    

%%R -w 600 -h 350
t.01.rel.clr.var = t.01.rel.clr %>% apply(1, var)
t.025.rel.clr.var = t.025.rel.clr %>% apply(1, var)
t.05.rel.clr.var = t.05.rel.clr %>% apply(1, var)

df.var = data.frame(t.01 = t.01.rel.clr.var,
                    t.025 = t.025.rel.clr.var,
                    t.05 = t.05.rel.clr.var)

matplot(locs, df.var, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Variance")
matplot(locs, df.var[,c(2,3)], lty = "solid", type = "l", pch = 1:10, cex = 0.8,
        xlab = "BD", ylab = "Variance")


    

In [ ]: