inst/doc/biodiversity-survey.R

In gradientForest: Random Forest functions for the Census of Marine Life synthesis project.

### R code from vignette source 'biodiversity-survey.rnw'

###################################################
### code chunk number 1: c1
###################################################
options(width=60) 
options(warn=-1) 


###################################################
### code chunk number 2: c2
###################################################
require(gradientForest)
# load the site by species matrix, object=Sp_mat
# includes transformed biomass of 110 species at 197 sites
load("GZ.sps.mat.Rdata") 
dim(Sp_mat)
# load the site physical data, object=Phys_site 
# includes 28 predictors at 197 sites
load("GZ.phys.site.Rdata") 
dim(Phys_site)


###################################################
### code chunk number 3: c3
###################################################
nSites <- dim(Sp_mat)[1]
nSpecs <- dim(Sp_mat)[2]
# set depth of conditional permutation
lev <- floor(log2(nSites*0.368/2))
lev 


###################################################
### code chunk number 4: c4 (eval = FALSE)
###################################################
## gf <- gradientForest(cbind(Phys_site,Sp_mat), predictor.vars=colnames(Phys_site), response.vars=colnames(Sp_mat), ntree=500, transform = NULL, compact=T, nbin = 201, maxLevel=lev, corr.threshold=0.5) 


###################################################
### code chunk number 5: c4b
###################################################
load("gf.Rdata") 


###################################################
### code chunk number 6: c4a
###################################################
gf
names(gf)


###################################################
### code chunk number 7: c5
###################################################
plot(gf,plot.type="O") 


###################################################
### code chunk number 8: c5
###################################################
# get a sorted list of the top 1:n most important predictors
most_important <- names(importance(gf))[1:25]
par(mgp=c(2, 0.75, 0))


###################################################
### code chunk number 9: c6
###################################################
plot(gf, plot.type="S", imp.vars=most_important, leg.posn="topright", cex.legend=0.4, cex.axis=0.6, cex.lab=0.7, line.ylab=0.9, par.args=list(mgp=c(1.5, 0.5, 0), mar=c(3.1,1.5,0.1,1)))


###################################################
### code chunk number 10: c7
###################################################
plot(gf, plot.type="C", imp.vars=most_important, show.overall=F, legend=T, leg.posn="topleft", leg.nspecies=5, cex.lab=0.7, cex.legend=0.4, cex.axis=0.6, line.ylab=0.9, par.args=list(mgp=c(1.5, 0.5, 0), mar=c(2.5,1.0,0.1,0.5), omi=c(0,0.3,0,0)))


###################################################
### code chunk number 11: c8
###################################################
plot(gf, plot.type="C", imp.vars=most_important, show.species=F, common.scale=T, cex.axis=0.6, cex.lab=0.7, line.ylab=0.9, par.args=list(mgp=c(1.5, 0.5, 0), mar=c(2.5,1.0,0.1,0.5), omi=c(0,0.3,0,0)))


###################################################
### code chunk number 12: c9
###################################################
plot(gf, plot.type="P", show.names=F, horizontal=F, cex.axis=1, cex.labels=0.7, line=2.5)


###################################################
### code chunk number 13: c10 (eval = FALSE)
###################################################
## plot(gf, plot.type="P", show.names=T, horizontal=F, cex.axis=1, cex.labels=0.7, line=2.5)
## plot(gf, plot.type="P", show.names=F, horizontal=T, cex.axis=1, cex.labels=0.6, line=2.5)
## plot(gf, plot.type="P", show.names=T, horizontal=T, cex.axis=1, cex.labels=0.6, line=2.5)


###################################################
### code chunk number 14: c11a
###################################################
load("GZ.phys.grid.Rdata") 


###################################################
### code chunk number 15: c11
###################################################
dim(Phys_grid) 
names(Phys_grid)


###################################################
### code chunk number 16: c12
###################################################
imp.vars <- names(importance(gf))
Trns_grid <- cbind(Phys_grid[,c("EAST","NORTH")], predict(gf,Phys_grid[,imp.vars])) 


###################################################
### code chunk number 17: c13
###################################################
Trns_site <- predict(gf) 


###################################################
### code chunk number 18: c14
###################################################
PCs <- prcomp(Trns_grid[,imp.vars]) 
# ensure BSTRESS points to upper right quadrant
sgn <- sign(PCs$rotation["BSTRESS",]) 
PCs$rotation <- sweep(PCs$rotation,2,sgn,"*")
PCs$x <- sweep(PCs$x,2,sgn,"*")
 
# set up a colour palette for the mapping
a1 <- PCs$x[,1]
a2 <- PCs$x[,2]
a3 <- PCs$x[,3]
r <- a1+a2
g <- -a2
b <- a3+a2-a1
r <- (r-min(r)) / (max(r)-min(r)) * 255
g <- (g-min(g)) / (max(g)-min(g)) * 255
b <- (b-min(b)) / (max(b)-min(b)) * 255


###################################################
### code chunk number 19: c15
###################################################
nvs <- dim(PCs$rotation)[1]  # number of variables
# choose predictor vectors to show
vec <- c("BSTRESS", "MUD", "SST_AV", "T_AV", "CHLA_AV", "SAND", "CRBNT", "GRAVEL") 
lv <- length(vec)
vind <- rownames(PCs$rotation) %in% vec

# choose a scaling factor to plot the vectors over the grid
scal <- 40 
xrng <- range(PCs$x[,1], PCs$rotation[,1]/scal)*1.1
yrng <- range(PCs$x[,2], PCs$rotation[,2]/scal)*1.1
plot((PCs$x[,1:2]), xlim=xrng, ylim=yrng, pch=".", cex=4, col=rgb(r,g,b, max = 255), asp=1)
# plot the other predictors with "+"
points(PCs$rotation[! vind,1:2]/scal, pch="+")  
# plot the chosen predictors as arrows
arrows(rep(0,lv), rep(0,lv), PCs$rotation[vec,1]/scal, PCs$rotation[vec,2]/scal, length = 0.0625) 
jit <- 0.0015
text(PCs$rotation[vec,1]/scal+jit*sign(PCs$rotation[vec,1]), PCs$rotation[vec,2]/scal+jit*sign(PCs$rotation[vec,2]), labels = vec)


###################################################
### code chunk number 20: c16a
###################################################
plot((PCs$x[,1:2]), xlim=xrng, ylim=yrng, pch=".", cex=4, col=rgb(r,g,b, max = 255), asp=1)
# plot the other predictors with "+"
points(PCs$rotation[! vind,1:2]/scal, pch="+")  
# plot the chosen predictors as arrows
arrows(rep(0,lv), rep(0,lv), PCs$rotation[vec,1]/scal, PCs$rotation[vec,2]/scal, length = 0.0625) 
jit <- 0.0015
text(PCs$rotation[vec,1]/scal+jit*sign(PCs$rotation[vec,1]), PCs$rotation[vec,2]/scal+jit*sign(PCs$rotation[vec,2]), labels = vec)
# first predict the PCs for the transformed site data
PCsites <- predict(PCs,Trns_site[,imp.vars])  
# plot all the sites as points on the biplot
points(PCsites[,1:2])
# calc & plot the weighted mean locations of species (from gf$Y) 
SpsWtd <- sweep(gf$Y,2,apply(gf$Y,2,min),"-")
SpsWtdPCs <- (t(SpsWtd) %*% (PCsites[,1:2]))/colSums(SpsWtd) 
points(SpsWtdPCs, col="red", pch="+") 
sp <- colnames(SpsWtd)[1]  
points(PCsites[,1:2], col="blue", cex=SpsWtd[,sp]/2)


###################################################
### code chunk number 21: c16 (eval = FALSE)
###################################################
## # first predict the PCs for the transformed site data
## PCsites <- predict(PCs,Trns_site[,imp.vars])  
## # plot all the sites as points on the biplot
## points(PCsites[,1:2])
## # calc & plot the weighted mean locations of species (from gf$Y) 
## SpsWtd <- sweep(gf$Y,2,apply(gf$Y,2,min),"-")
## SpsWtdPCs <- (t(SpsWtd) %*% (PCsites[,1:2]))/colSums(SpsWtd) 
## points(SpsWtdPCs, col="red", pch="+") 
## # interactively label some of the species, if desired
## #identify(SpsWtdPCs, labels = as.character(rownames(SpsWtdPCs)), col="blue")  


###################################################
### code chunk number 22: c17
###################################################
sp <- colnames(SpsWtd)[1]  
points(PCsites[,1:2], col="blue", cex=SpsWtd[,sp]/2)


###################################################
### code chunk number 23: c18
###################################################
plot(Trns_grid[,c("EAST","NORTH")],pch=".", cex=3, asp=1, col=rgb(r,g,b, max = 255))


###################################################
### code chunk number 24: c19
###################################################
require(cluster)
ncl <- 8 
clPCs <- clara(PCs$x,ncl,sampsize=1000) 
# get the medoid colour palette
medcolR <- r[clPCs$i.med]
medcolG <- g[clPCs$i.med]
medcolB <- b[clPCs$i.med]

# re-plot the biplot -- coloured by cluster medoids
plot((PCs$x[,1:2]), xlim=xrng, ylim=yrng,pch=".", cex=4, col=rgb(medcolR[clPCs$clustering], medcolG[clPCs$clustering], medcolB[clPCs$clustering],  max = 255), asp=1)
points(PCs$rotation[! vind,1:2]/scal, pch="+") 
arrows(rep(0,lv), rep(0,lv), PCs$rotation[vec,1]/scal, PCs$rotation[vec,2]/scal, length = 0.0625) 
text(PCs$rotation[vec,1]/scal+jit*sign(PCs$rotation[vec,1]), PCs$rotation[vec,2]/scal+jit*sign(PCs$rotation[vec,2]), labels = vec)
# plot the cluster medoids with ID number
text(clPCs$medoids[,1:2], labels = seq(1,ncl)) 
legend("bottomleft",as.character(seq(1,ncl)), pch=15, cex=1,col=rgb(medcolR,medcolG,medcolB, max = 255))


###################################################
### code chunk number 25: c20
###################################################
# re-plot the map --- coloured by cluster medoids
plot(Trns_grid[,c("EAST","NORTH")],pch=".", cex=3, asp=1, col=rgb(medcolR[clPCs$clustering], medcolG[clPCs$clustering], medcolB[clPCs$clustering], max = 255))
points(Trns_grid[clPCs$i.med,c("EAST","NORTH")], pch=as.character(seq(1,ncl)))
legend("bottomleft",as.character(seq(1,ncl)), pch=15, cex=1, col=rgb(medcolR,medcolG,medcolB, max = 255))


###################################################
### code chunk number 26: biodiversity-survey.rnw:395-396
###################################################
toLatex(sessionInfo())