In [3]:
%load_ext rpy2.ipython
In [4]:
%%R
library(coenocliner)
library(compositions)
library(dplyr)
library(tidyr)
library(ggplot2)
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
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
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 [134]:
%%R
n.taxa = 500
n.fracs = 500
In [142]:
%%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
tol <- runif(M, min=0.01, max=0.025) # species tolerances
h <- ceiling(rlnorm(M, meanlog = 12)) # max abundances
pars <- cbind(opt = opt, tol = tol, h = h) # put in a matrix
In [143]:
%%R -w 600 -h 350
t.1 = coenocline(locs, responseModel = "gaussian", params = pars, countModel = "poisson")
matplot(locs, t.1, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
xlab = "BD", ylab = "Abundance")
In [159]:
%%R
opt <- rnorm(M, mean=meang-0.01, sd=0.01) # species optima
h <- ceiling(rlnorm(M, meanlog = 12)) # max abundances
pars <- cbind(opt = opt, tol = tol, h = h) # put in a matrix
In [160]:
%%R -w 600 -h 350
t.2 = coenocline(locs, responseModel = "gaussian", params = pars, countModel = "poisson")
matplot(locs, t.2, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
xlab = "BD", ylab = "Abundance")
In [161]:
%%R -w 600 -h 350
to.rel = function(x){
y = sum(x)
if(y == 0){
return(x)
} else {
return(x / y)
}
}
t.1.rel = t(apply(t.1, 1, to.rel))
t.2.rel = t(apply(t.2, 1, to.rel))
matplot(locs, t.1.rel, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
xlab = "BD", ylab = "Relative abundance", main='comm1')
matplot(locs, t.2.rel, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
xlab = "BD", ylab = "Relative abundance", main='comm2')
In [162]:
%%R
t.1.rel %>% apply(1, var) %>% sum %>% print
t.2.rel %>% apply(1, var) %>% sum %>% print
In [163]:
%%R -w 600 -h 350
t.1.rel.var = t.1.rel %>% apply(1, var)
t.2.rel.var = t.2.rel %>% apply(1, var)
df.var = data.frame(t.1 = t.1.rel.var,
t.2 = t.2.rel.var)
matplot(locs, df.var, lty = "solid", type = "l", pch = 1:10, cex = 0.8,
xlab = "BD", ylab = "Variance")
In [ ]:
In [ ]: