In [1]:
%load_ext rpy2.ipython
In [173]:
%%R
library(coenocliner)
library(compositions)
library(dplyr)
library(tidyr)
library(ggplot2)
In [3]:
%%R
set.seed(2)
M <- 3 # number of species
ming <- 3.5 # gradient minimum...
maxg <- 7 # ...and maximum
meang = mean(c(ming, maxg))
locs <- seq(ming, maxg, length = 100) # gradient locations
opt <- rnorm(M, mean=meang, sd=0.5) #runif(M, min = ming, max = maxg) # species optima
tol <- rep(0.25, M) # species tolerances
h <- ceiling(rlnorm(M, meanlog = 10)) # max abundances
pars <- cbind(opt = opt, tol = tol, h = h) # put in a matrix
In [4]:
%%R
mu1 <- coenocline(locs, responseModel = "gaussian", params = pars,
countModel = "poisson") #, countParams = list(alpha = 0.5))
head(mu1)
In [5]:
%%R -w 600 -h 350
matplot(locs, mu1, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [6]:
%%R
set.seed(2)
M <- 2 # number of species
ming <- 3.5 # gradient minimum...
maxg <- 7 # ...and maximum
meang = mean(c(ming, maxg)) + 0.3
locs <- seq(ming, maxg, length = 100) # gradient locations
opt <- rnorm(M, mean=meang, sd=0.5) #runif(M, min = ming, max = maxg) # species optima
tol <- rep(0.25, M) # species tolerances
h <- ceiling(rlnorm(M, meanlog = 9.5)) # max abundances
pars <- cbind(opt = opt, tol = tol, h = h) # put in a matrix
In [7]:
%%R
mu2 <- coenocline(locs, responseModel = "gaussian", params = pars,
countModel = "poisson") #, countParams = list(alpha = 0.5))
mu2 = cbind(mu1, mu2)
head(mu2)
In [8]:
%%R -w 600 -h 350
matplot(locs, mu2, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [9]:
%%R
to.rel = function(x){
y = sum(x)
if(y == 0){
return(x)
} else {
return(x / y)
}
}
mu1_rel = t(apply(mu1, 1, to.rel))
mu2_rel = t(apply(mu2, 1, to.rel))
# check
apply(mu1_rel, 1, sum)
In [10]:
%%R -w 600 -h 350
matplot(locs, mu1_rel, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [11]:
%%R -w 600 -h 350
matplot(locs, mu2_rel, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [12]:
%%R -w 600 -h 350
mu1_clr = clr(mu1_rel)
matplot(locs, mu1_clr, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [13]:
%%R -w 600 -h 350
mu2_clr = clr(mu2_rel)
matplot(locs, mu2_clr, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [14]:
%%R -w 600 -h 350
matplot(locs, mu2_clr[,c(3,5)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [15]:
%%R
apply(mu2, 2, sum)
In [16]:
%%R -w 600 -h 350
mu2_alr = alr(mu2_rel, ivar=5)
matplot(locs, mu2_alr, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [17]:
%%R -w 600 -h 350
mu1_rel_nZ = mu1_rel
mu1_rel_nZ[mu1_rel_nZ == 0] = min(mu1_rel[mu1_rel > 0]) / 2
mu1_clr = clr(mu1_rel_nZ)
matplot(locs, mu1_clr, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [18]:
%%R -w 600 -h 350
mu2_rel_nZ = mu2_rel
mu2_rel_nZ[mu2_rel_nZ == 0] = min(mu2_rel[mu2_rel > 0]) / 2
mu2_clr = clr(mu2_rel_nZ)
matplot(locs, mu2_clr, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [19]:
%%R -w 600 -h 350
mu2_alr = alr(mu2_rel_nZ, ivar=5)
matplot(locs, mu2_alr, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [20]:
%%R -w 600 -h 350
ss_fit = function(y, x){
#ss = smooth.spline(x,y)
#return(ss$y)
ss = pspline::smooth.Pspline(x,y,method=3)
#plot(ss)
return(ss$ysmth)
}
mu2_clrSsInv = clrInv(apply(mu2_clr, 2, ss_fit, x=locs))
matplot(locs, mu2_clrSsInv, lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [21]:
%%R
center.mass = function(y,x){
wm = weighted.mean(x, y)
return(wm)
}
apply(mu2, 2, center.mass, x=locs) %>% print
apply(mu2_rel, 2, center.mass, x=locs) %>% print
apply(mu2_clr, 2, center.mass, x=locs) %>% print
apply(mu2_clrSsInv, 2, center.mass, x=locs) %>% print
In [ ]:
In [22]:
%%R -w 600 -h 350
ss3 = smooth.spline(locs, mu2_clr[,3])
ss5 = smooth.spline(locs, mu2_clr[,5])
matplot(locs, mu2_clr[,c(3,5)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
lines(ss3, lty=2, col='black')
lines(ss5, lty=2, col='red')
In [23]:
%%R -w 600 -h 350
plot_clr_spline = function(col, mu, mu_rel, mu_clr, locs){
ss = smooth.spline(locs, mu_clr[,col])
plot(locs, mu[,col])
plot(locs, mu_rel[,col])
plot(locs, mu_clr[,col])
plot(locs, ss$y)
plot(locs, clrInv(ss$y))
}
plot_clr_spline(1, mu2, mu2_rel, mu2_clr, locs)
In [24]:
%%R -w 600 -h 350
plot_clr_spline(2, mu2, mu2_rel, mu2_clr, locs)
In [25]:
%%R -w 600 -h 350
plot_clr_spline(3, mu2, mu2_rel, mu2_clr, locs)
In [26]:
%%R -w 600 -h 350
plot_clr_spline(4, mu2, mu2_rel, mu2_clr, locs)
In [27]:
%%R -w 600 -h 350
plot_clr_spline(5, mu2, mu2_rel, mu2_clr, locs)
In [29]:
%%R
center.mass = function(locs, mu, col, ss=FALSE){
if(ss==TRUE){
y = smooth.spline(locs, mu[,col])
y = clrInv(y$y)
} else {
y = mu[,col]
}
wm = weighted.mean(locs, y)
print(wm)
}
col = 1
center.mass(locs, mu2, col)
center.mass(locs, mu2_rel, col)
center.mass(locs, mu2_clr, col, ss=TRUE)
In [30]:
%%R
col = 2
center.mass(locs, mu2, col)
center.mass(locs, mu2_rel, col)
center.mass(locs, mu2_clr, col, ss=TRUE)
In [31]:
%%R
col = 3
center.mass(locs, mu2, col)
center.mass(locs, mu2_rel, col)
center.mass(locs, mu2_clr, col, ss=TRUE)
In [32]:
%%R
col = 4
center.mass(locs, mu2, col)
center.mass(locs, mu2_rel, col)
center.mass(locs, mu2_clr, col, ss=TRUE)
In [33]:
%%R
col = 5
center.mass(locs, mu2, col)
center.mass(locs, mu2_rel, col)
center.mass(locs, mu2_clr, col, ss=TRUE)
In [34]:
%%R
col = 6
center.mass(locs, mu2, col)
center.mass(locs, mu2_rel, col)
center.mass(locs, mu2_clr, col, ss=TRUE)
In [35]:
%%R
# mean
x = acomp(mu2_rel)
mean(x)
In [36]:
%%R
mvar(x) %>% print
msd(x) %>% print
In [37]:
%%R
variation(x)
In [38]:
%%R -h 300
matplot(locs, mu2[,c(1,4)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
matplot(locs, mu2_rel[,c(1,4)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
matplot(locs, mu2_clr[,c(1,4)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [39]:
%%R -h 300
matplot(locs, mu2[,c(1,3)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
matplot(locs, mu2_rel[,c(1,3)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
matplot(locs, mu2_clr[,c(1,3)], lty = "solid", type = "p", pch = 1:10, cex = 0.8,
xlab = "pH", ylab = "Abundance")
In [156]:
%%R
library(vegan)
In [210]:
%%R -w 800
locs.ap = expand.grid(locs, locs)
locs.ap$diff = locs.ap %>% apply(1, diff) %>% abs %>% as.vector
locs.ap = locs.ap %>%
spread(Var1, diff)
rownames(locs.ap) = locs.ap$Var2
locs.ap$Var2 = NULL
locs.dist = locs.ap %>% as.matrix
locs.dist[upper.tri(locs.dist, diag=TRUE)] = 0
locs.dist = locs.dist %>% as.dist
plot(hclust(locs.dist))
In [226]:
%%R
mu1.dist = vegdist(mu1)
mu1.dist[is.nan(mu1.dist)] = 0
stopifnot((mu2.dist %>% length) == (locs.dist %>% length))
mu1.mcorr = mantel.correlog(mu1.dist, locs.dist, n.class=24)
plot(mu1.mcorr)
In [227]:
%%R
mu2.dist = vegdist(mu2)
mu2.dist[is.nan(mu2.dist)] = 0
stopifnot((mu2.dist %>% length) == (locs.dist %>% length))
mu2.mcorr = mantel.correlog(mu2.dist, locs.dist, n.class=24)
plot(mu2.mcorr)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [40]:
%%R
library(robCompositions)
In [41]:
%%R
data(expenditures)
x = xOrig = expenditures
head(x)
In [42]:
%%R
to.rel = function(x){
x[is.na(x)] = 0
y = sum(x)
if(is.na(y) | y == 0){
x = rep(0, length(x))
return(x)
} else {
return(x / y)
}
}
x_rel = apply(x, 1, function(x) to.rel(x))
x_rel = t(x_rel)
head(x_rel)
In [43]:
%%R
## calculate the relative Aitchsion distance between xOrig and xImp:
aDist(x_rel[1,], x_rel[2,])
In [107]:
%%R
.aDist = function(x, df){
i = x[1]
ii = x[2]
aDist(df[i,], df[ii,])
}
aDist.all = function(df, meta=NULL){
# pairwise aDist on all rows
i.ii = combn(1:nrow(df), m=2)
vals = apply(i.ii, 2, .aDist, df=df)
vals = cbind(i.ii %>% t, vals)
vals = as.data.frame(vals)
colnames(vals) = c('V1', 'V2', 'aDist')
# adding metadata if provided
if(! is.null(meta)){
if(class(meta) == 'matrix'){
meta = as.data.frame(meta)
}
if(class(meta) != 'data.frame'){
meta = data.frame(1:length(meta), meta)
colnames(meta)[1] = c('V1')
}
name1 = colnames(vals)[1]
name2 = colnames(meta)[1]
vals = inner_join(vals, meta, c(V1 = name2))
vals = inner_join(vals, meta, c(V2 = name2))
}
return(vals)
}
test = mu2_rel_nZ[1:5,]
test.locs = data.frame(1:5, locs[1:5])
colnames(test.locs) = c('V1', 'pH')
aDist.all(test, test.locs)
In [108]:
%%R
mu2_rel_nZ_aDist = aDist.all(mu2_rel_nZ)
ggplot(mu2_rel_nZ_aDist, aes(V1, V2, fill=aDist)) +
geom_tile()
In [ ]:
In [142]:
%%R
library(gstat)
library(sp)
data(meuse)
coordinates(meuse) = ~x+y
lzn.vgm = variogram(log(zinc)~1, meuse)
lzn.fit = fit.variogram(lzn.vgm, model = vgm(1, "Sph", 900, 1))
plot(lzn.vgm, lzn.fit)
In [143]:
%%R
?vgram.sph
In [110]:
%%R
df.locs = data.frame(1:length(locs), locs)
colnames(df.locs) = c('V1', 'pH')
mu2_rel_nZ_aDist = aDist.all(mu2_rel_nZ, df.locs) %>%
mutate(delta.pH = abs(pH.x - pH.y))
ggplot(mu2_rel_nZ_aDist, aes(delta.pH, aDist)) +
geom_point()
In [146]:
%%R
data(SimulatedAmounts)
meanCol(sa.lognormals)
variation(acomp(sa.lognormals))