R/plotBeta.R
In hmsc-r/HMSC: Hierarchical Model of Species Communities
Defines functions
plotBeta
Documented in
plotBeta
#' @title plotBeta
#'
#' @description Plots heatmaps of parameter estimates or posterior support values of species' environmental responses, i.e. how species in \code{Y} responds to covariates in \code{X}
#'
#' @param hM a fitted \code{Hmsc} model object
#' @param post posterior summary of Beta parameters obtained from \code{\link{getPostEstimate}}
#' @param param controls which parameter is plotted, current options include "Mean" for posterior mean estimate,
#' "Support" for the level of statistical support measured by posterior probability for a positive or negative
#' response, and "Sign" to indicate whether the response is positive, negative, or neither of these given the
#' chosen \code{supportLevel}
#' @param plotTree logical. Whether species' environmental responses is to be mapped onto the phylogeny used in
#' model fitting
#' @param SpeciesOrder controls the ordering of species, current options are "Original", "Tree", and "Vector".
#' If SpeciesOrder = "Vector", an ordering vector must be provided (see SpVector). If plotTree = TRUE, SpeciesOrder
#' is ignored
#' @param SpVector controls the ordering of species if SpeciesOrder = "Vector". If a subset of species are listed,
#' only those will be plotted. For alphabetic ordering, try
#' \code{match(1:hM$ns, as.numeric(as.factor(colnames(hM$Y))))}
#' @param covOrder controls the ordering of covariates, current options are "Original" and "Vector".
#' If covOrder = "Vector", an ordering vector must be provided (see covVector)
#' @param covVector controls the ordering of covariates if covOrder = "Vector". If a subset of covariates are listed,
#' only those will be plotted
#' @param spNamesNumbers logical of length 2, where first entry controls whether species names are added to axes,
#' and second entry controls whether species numbers are added
#' @param covNamesNumbers logical of length 2, where first entry controls whether covariate names are added to axes,
#' and second entry controls whether covariate numbers are added
#' @param supportLevel controls threshold posterior support for plotting
#' @param split if plotTree = TRUE, controls the division of the plotting window between the tree and the heatmap.
#' @param cex controls character expansion (font size). Three values, controlling covariate names, species names,
#' and color legend axis labels
#' @param colors controls the colors of the heatmap, default value \code{colorRampPalette(c("blue","white","red"))}
#' @param colorLevels number of color levels used in the heatmap
#' @param mar plotting margins
#' @param marTree plotting margins for phylogenetic tree
#' @param mgp can be used to set the location of the scale bar
#' @param main main title for the plot.
#' @param smallplot passed to \code{\link{image.plot}}
#' @param bigplot passed to \code{\link{image.plot}}
#' @param newplot set to false if the plot will be part of multi-panel plot initialized with par(mfrow)
#'
#'
#'
#' @examples
#' # Plot posterior support values of species' environmental responses
#' betaPost=getPostEstimate(TD$m, "Beta")
#' plotBeta(TD$m, post=betaPost, param="Support")
#'
#' # Plot parameter estimates of species' environmental responses together with the phylogenetic tree
#' betaPost=getPostEstimate(TD$m, "Beta")
#' plotBeta(TD$m, post=betaPost, param="Mean", plotTree=TRUE)
#'
#' @importFrom graphics par plot plot.new axis text title
#' @importFrom grDevices colorRampPalette
#' @importFrom ape keep.tip read.tree write.tree
#' @importFrom fields image.plot
#' @export
plotBeta = function(hM, post, param = "Support", plotTree = FALSE,
SpeciesOrder = "Original", SpVector = NULL,
covOrder="Original",
covVector=NULL, spNamesNumbers = c(TRUE, TRUE),
covNamesNumbers = c(TRUE, TRUE),
supportLevel = 0.9, split = 0.3, cex = c(0.7,0.7,0.8),
colors = colorRampPalette(c("blue","white","red")), colorLevels = NULL,
mar=NULL, marTree=c(6,0,2,0),mgp=c(3,2,0),
main = NULL,
smallplot=NULL, bigplot=NULL,newplot=TRUE)
{
## Check that text arguments are acceptable, and expand to full
## names if abbreviated.
param <- match.arg(param, c("Mean", "Support", "Sign"))
SpeciesOrder <- match.arg(SpeciesOrder, c("Original", "Vector", "Tree"))
covOrder <- match.arg(covOrder, c("Original", "Vector"))
if(is.null(colorLevels)){
if(param=="Sign"){
colorLevels=3} else {
colorLevels=200
}
}
if(is.null(mar)){
if(sum(spNamesNumbers)>0){
mar=c(6,6,2,0)
} else {
mar=c(6,0,2,0)
}
}
if(plotTree || SpeciesOrder == "Tree"){
untangle<-function(tree){
if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
obj<-attributes(tree)
tree <- ape::read.tree(text=ape::write.tree(tree))
ii <- !names(obj) %in% names(attributes(tree))
attributes(tree)<-c(attributes(tree),obj[ii])
tree
}
tree = keep.tip(hM$phyloTree,hM$spNames)
tree = untangle(tree)
}
spNames = character(hM$ns)
for (i in 1:hM$ns) {
sep = ""
if (spNamesNumbers[1]) {
spNames[i] = paste(spNames[i], hM$spNames[i], sep = sep)
sep = " "
}
if (spNamesNumbers[2]) {
spNames[i] = paste(spNames[i], sprintf("(S%d)", i),
sep = sep)
}
}
covNames = character(hM$nc)
for (i in 1:hM$nc) {
sep = ""
if (covNamesNumbers[1]) {
covNames[i] = paste(covNames[i], hM$covNames[i], sep = sep)
sep = " "
}
if (covNamesNumbers[2]) {
covNames[i] = paste(covNames[i], sprintf("(C%d)",
i), sep = sep)
}
}
if(plotTree){
order=tree$tip.label
} else {
order <-
switch(SpeciesOrder,
"Vector" = SpVector,
"Original" = rev(1:ncol(hM$Y)),
"Tree" = match(tree$tip.label,colnames(hM$Y)))
}
if(covOrder=="Vector"){covorder=covVector}
else if(covOrder=="Original"){covorder=1:hM$nc}
mbeta=post$mean
betaP=post$support
if(param=="Sign"){
toPlot = sign(mbeta)
toPlot = toPlot * ((betaP>supportLevel) + (betaP<(1-supportLevel))>0)
betaMat = matrix(toPlot, nrow=hM$nc, ncol=ncol(hM$Y))
}
else if(param=="Mean"){
toPlot = mbeta
toPlot = toPlot * ((betaP>supportLevel) + (betaP<(1-supportLevel))>0)
betaMat = matrix(toPlot, nrow=hM$nc, ncol=ncol(hM$Y))
}
else{ # param == "Support"
toPlot = 2*betaP-1
toPlot = toPlot * ((betaP>supportLevel) + (betaP<(1-supportLevel))>0)
betaMat = matrix(toPlot, nrow=hM$nc, ncol=ncol(hM$Y))
}
rownames(betaMat) = covNames
if(plotTree){colnames(betaMat) = gsub(" ", "_", hM$spNames)}
if(!plotTree){colnames(betaMat) = spNames}
X = t(betaMat[covorder,order])
old.par = par(no.readonly = TRUE)
colors = colors(colorLevels)
#With tree
if(plotTree){
if(newplot){
par(fig = c(0,split[1],0,1), mar=marTree)
} else {
plot.new()
cpar = par(no.readonly = TRUE)
cfig = cpar$fig
nfig = cfig
nfig[2] = cfig[1]+split[1]*(cfig[2]-cfig[1])
par(fig=nfig, mar=marTree, new=TRUE)
}
if(sum(!spNamesNumbers)==2){plot(tree, show.tip.label=FALSE)}
else{tree$tip.label[match(gsub(" ", "_", hM$spNames),tree$tip.label)]=spNames
plot(tree, show.tip.label=TRUE, adj=1, align.tip.label=TRUE, cex=cex[2])}
if(newplot){
par(fig = c(split[1],1,0,1), mar=mar, new=TRUE)
} else {
nfig = cfig
nfig[1] = cfig[1]+split[1]*(cfig[2]-cfig[1])
par(fig=nfig, new=TRUE)
par(old.par, mar=mar, new=TRUE)
}
START = .05
END = .7
ADJy=0
ADJx=1/(ncol(X)*4)
plot.new()
axis(1,seq((START+ADJx), (END-ADJx),
by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)), labels = FALSE)
text(x=seq((START+ADJx), (END-ADJx), by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)), par("usr")[3] - 0.05, srt = 90, adj = 1,cex=cex[1],
labels = covNames[covorder], xpd = TRUE)
print(par(no.readonly = TRUE)$fig)
}
#No tree
if(!plotTree){
START=0
END=.65
ADJy=1/(nrow(X)*4)
ADJx=1/(ncol(X)*4)
plot.new()
if(newplot){
par(fig = c(0,1,0,1), mar = mar)
} else {
cpar = par(no.readonly = TRUE)
cfig = cpar$fig
par(fig=cfig, mar=mar, new=TRUE)
}
axis(1,at = seq((START+ADJx), (END-ADJx),
by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)),
labels = FALSE)
text(x=seq((START+ADJx), (END-ADJx), by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)), par("usr")[3] - 0.05, srt = 90, adj = 1,cex=cex[1],
labels = covNames[covorder], xpd = TRUE)
names=gsub("_"," ",spNames[order])
text(y = seq(ADJy, 1-ADJy, length.out=nrow(X)),par("usr")[3] - 0.05, srt = 0, adj = 1,cex=cex[2],
labels = as.expression(lapply(names, function(names) bquote(italic(.(names))))), xpd = TRUE)
}
#Plot
if(all(is.na(X)) || sum(abs(X))==0){
warning("nothing to plot at this level of posterior support")
zlim = c(-1,1)
} else{
zlim = c(-max(abs(range(X))),max(abs(range(X))))
}
image.plot(x = seq(START+ADJx, END-ADJx, by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)),
y = seq(ADJy, 1-ADJy, length.out=nrow(X)),
z = t(X), add = TRUE, nlevel = colorLevels, box=TRUE,
legend.width = 2, legend.mar = NULL,
legend.cex=cex,
axis.args=if(param=="Sign")
{list(labels=c("+","0","-"),at=c(1,0,-1),cex.axis=cex[3],mgp=mgp,hadj=1)
} else {
list(cex.axis=cex[3],mgp=mgp,hadj=1)
},
graphics.reset = TRUE, horizontal = FALSE, bigplot = bigplot, smallplot = smallplot,
legend.only = FALSE, col = colors,zlim = zlim)
if (!is.null(main))
title(main = main)
if(newplot){
par(old.par)
} else {
par(mfrow=cpar$mfrow)
cmfg=cpar$mfg
if(cmfg[2]<cmfg[4]){
cmfg[2] = cmfg[2]+1
} else {
if(cmfg[1]<cmfg[3]){
cmfg[2]=1
cmfg[1] = cmfg[1]+1
}
else{
cmfg[1]=1
cmfg[2]=1
}
}
par(mfg=cmfg,mar=cpar$mar)
}
}
hmsc-r/HMSC documentation built on March 5, 2025, 10:52 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plotBeta
Documented in plotBeta
#' @title plotBeta
#'
#' @description Plots heatmaps of parameter estimates or posterior support values of species' environmental responses, i.e. how species in \code{Y} responds to covariates in \code{X}
#'
#' @param hM a fitted \code{Hmsc} model object
#' @param post posterior summary of Beta parameters obtained from \code{\link{getPostEstimate}}
#' @param param controls which parameter is plotted, current options include "Mean" for posterior mean estimate,
#' "Support" for the level of statistical support measured by posterior probability for a positive or negative
#' response, and "Sign" to indicate whether the response is positive, negative, or neither of these given the
#' chosen \code{supportLevel}
#' @param plotTree logical. Whether species' environmental responses is to be mapped onto the phylogeny used in
#' model fitting
#' @param SpeciesOrder controls the ordering of species, current options are "Original", "Tree", and "Vector".
#' If SpeciesOrder = "Vector", an ordering vector must be provided (see SpVector). If plotTree = TRUE, SpeciesOrder
#' is ignored
#' @param SpVector controls the ordering of species if SpeciesOrder = "Vector". If a subset of species are listed,
#' only those will be plotted. For alphabetic ordering, try
#' \code{match(1:hM$ns, as.numeric(as.factor(colnames(hM$Y))))}
#' @param covOrder controls the ordering of covariates, current options are "Original" and "Vector".
#' If covOrder = "Vector", an ordering vector must be provided (see covVector)
#' @param covVector controls the ordering of covariates if covOrder = "Vector". If a subset of covariates are listed,
#' only those will be plotted
#' @param spNamesNumbers logical of length 2, where first entry controls whether species names are added to axes,
#' and second entry controls whether species numbers are added
#' @param covNamesNumbers logical of length 2, where first entry controls whether covariate names are added to axes,
#' and second entry controls whether covariate numbers are added
#' @param supportLevel controls threshold posterior support for plotting
#' @param split if plotTree = TRUE, controls the division of the plotting window between the tree and the heatmap.
#' @param cex controls character expansion (font size). Three values, controlling covariate names, species names,
#' and color legend axis labels
#' @param colors controls the colors of the heatmap, default value \code{colorRampPalette(c("blue","white","red"))}
#' @param colorLevels number of color levels used in the heatmap
#' @param mar plotting margins
#' @param marTree plotting margins for phylogenetic tree
#' @param mgp can be used to set the location of the scale bar
#' @param main main title for the plot.
#' @param smallplot passed to \code{\link{image.plot}}
#' @param bigplot passed to \code{\link{image.plot}}
#' @param newplot set to false if the plot will be part of multi-panel plot initialized with par(mfrow)
#'
#'
#'
#' @examples
#' # Plot posterior support values of species' environmental responses
#' betaPost=getPostEstimate(TD$m, "Beta")
#' plotBeta(TD$m, post=betaPost, param="Support")
#'
#' # Plot parameter estimates of species' environmental responses together with the phylogenetic tree
#' betaPost=getPostEstimate(TD$m, "Beta")
#' plotBeta(TD$m, post=betaPost, param="Mean", plotTree=TRUE)
#'
#' @importFrom graphics par plot plot.new axis text title
#' @importFrom grDevices colorRampPalette
#' @importFrom ape keep.tip read.tree write.tree
#' @importFrom fields image.plot
#' @export
plotBeta = function(hM, post, param = "Support", plotTree = FALSE,
SpeciesOrder = "Original", SpVector = NULL,
covOrder="Original",
covVector=NULL, spNamesNumbers = c(TRUE, TRUE),
covNamesNumbers = c(TRUE, TRUE),
supportLevel = 0.9, split = 0.3, cex = c(0.7,0.7,0.8),
colors = colorRampPalette(c("blue","white","red")), colorLevels = NULL,
mar=NULL, marTree=c(6,0,2,0),mgp=c(3,2,0),
main = NULL,
smallplot=NULL, bigplot=NULL,newplot=TRUE)
{
## Check that text arguments are acceptable, and expand to full
## names if abbreviated.
param <- match.arg(param, c("Mean", "Support", "Sign"))
SpeciesOrder <- match.arg(SpeciesOrder, c("Original", "Vector", "Tree"))
covOrder <- match.arg(covOrder, c("Original", "Vector"))
if(is.null(colorLevels)){
if(param=="Sign"){
colorLevels=3} else {
colorLevels=200
}
}
if(is.null(mar)){
if(sum(spNamesNumbers)>0){
mar=c(6,6,2,0)
} else {
mar=c(6,0,2,0)
}
}
if(plotTree || SpeciesOrder == "Tree"){
untangle<-function(tree){
if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
obj<-attributes(tree)
tree <- ape::read.tree(text=ape::write.tree(tree))
ii <- !names(obj) %in% names(attributes(tree))
attributes(tree)<-c(attributes(tree),obj[ii])
tree
}
tree = keep.tip(hM$phyloTree,hM$spNames)
tree = untangle(tree)
}
spNames = character(hM$ns)
for (i in 1:hM$ns) {
sep = ""
if (spNamesNumbers[1]) {
spNames[i] = paste(spNames[i], hM$spNames[i], sep = sep)
sep = " "
}
if (spNamesNumbers[2]) {
spNames[i] = paste(spNames[i], sprintf("(S%d)", i),
sep = sep)
}
}
covNames = character(hM$nc)
for (i in 1:hM$nc) {
sep = ""
if (covNamesNumbers[1]) {
covNames[i] = paste(covNames[i], hM$covNames[i], sep = sep)
sep = " "
}
if (covNamesNumbers[2]) {
covNames[i] = paste(covNames[i], sprintf("(C%d)",
i), sep = sep)
}
}
if(plotTree){
order=tree$tip.label
} else {
order <-
switch(SpeciesOrder,
"Vector" = SpVector,
"Original" = rev(1:ncol(hM$Y)),
"Tree" = match(tree$tip.label,colnames(hM$Y)))
}
if(covOrder=="Vector"){covorder=covVector}
else if(covOrder=="Original"){covorder=1:hM$nc}
mbeta=post$mean
betaP=post$support
if(param=="Sign"){
toPlot = sign(mbeta)
toPlot = toPlot * ((betaP>supportLevel) + (betaP<(1-supportLevel))>0)
betaMat = matrix(toPlot, nrow=hM$nc, ncol=ncol(hM$Y))
}
else if(param=="Mean"){
toPlot = mbeta
toPlot = toPlot * ((betaP>supportLevel) + (betaP<(1-supportLevel))>0)
betaMat = matrix(toPlot, nrow=hM$nc, ncol=ncol(hM$Y))
}
else{ # param == "Support"
toPlot = 2*betaP-1
toPlot = toPlot * ((betaP>supportLevel) + (betaP<(1-supportLevel))>0)
betaMat = matrix(toPlot, nrow=hM$nc, ncol=ncol(hM$Y))
}
rownames(betaMat) = covNames
if(plotTree){colnames(betaMat) = gsub(" ", "_", hM$spNames)}
if(!plotTree){colnames(betaMat) = spNames}
X = t(betaMat[covorder,order])
old.par = par(no.readonly = TRUE)
colors = colors(colorLevels)
#With tree
if(plotTree){
if(newplot){
par(fig = c(0,split[1],0,1), mar=marTree)
} else {
plot.new()
cpar = par(no.readonly = TRUE)
cfig = cpar$fig
nfig = cfig
nfig[2] = cfig[1]+split[1]*(cfig[2]-cfig[1])
par(fig=nfig, mar=marTree, new=TRUE)
}
if(sum(!spNamesNumbers)==2){plot(tree, show.tip.label=FALSE)}
else{tree$tip.label[match(gsub(" ", "_", hM$spNames),tree$tip.label)]=spNames
plot(tree, show.tip.label=TRUE, adj=1, align.tip.label=TRUE, cex=cex[2])}
if(newplot){
par(fig = c(split[1],1,0,1), mar=mar, new=TRUE)
} else {
nfig = cfig
nfig[1] = cfig[1]+split[1]*(cfig[2]-cfig[1])
par(fig=nfig, new=TRUE)
par(old.par, mar=mar, new=TRUE)
}
START = .05
END = .7
ADJy=0
ADJx=1/(ncol(X)*4)
plot.new()
axis(1,seq((START+ADJx), (END-ADJx),
by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)), labels = FALSE)
text(x=seq((START+ADJx), (END-ADJx), by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)), par("usr")[3] - 0.05, srt = 90, adj = 1,cex=cex[1],
labels = covNames[covorder], xpd = TRUE)
print(par(no.readonly = TRUE)$fig)
}
#No tree
if(!plotTree){
START=0
END=.65
ADJy=1/(nrow(X)*4)
ADJx=1/(ncol(X)*4)
plot.new()
if(newplot){
par(fig = c(0,1,0,1), mar = mar)
} else {
cpar = par(no.readonly = TRUE)
cfig = cpar$fig
par(fig=cfig, mar=mar, new=TRUE)
}
axis(1,at = seq((START+ADJx), (END-ADJx),
by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)),
labels = FALSE)
text(x=seq((START+ADJx), (END-ADJx), by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)), par("usr")[3] - 0.05, srt = 90, adj = 1,cex=cex[1],
labels = covNames[covorder], xpd = TRUE)
names=gsub("_"," ",spNames[order])
text(y = seq(ADJy, 1-ADJy, length.out=nrow(X)),par("usr")[3] - 0.05, srt = 0, adj = 1,cex=cex[2],
labels = as.expression(lapply(names, function(names) bquote(italic(.(names))))), xpd = TRUE)
}
#Plot
if(all(is.na(X)) || sum(abs(X))==0){
warning("nothing to plot at this level of posterior support")
zlim = c(-1,1)
} else{
zlim = c(-max(abs(range(X))),max(abs(range(X))))
}
image.plot(x = seq(START+ADJx, END-ADJx, by = ((END-ADJx) - (START+ADJx))/(ncol(X) - 1)),
y = seq(ADJy, 1-ADJy, length.out=nrow(X)),
z = t(X), add = TRUE, nlevel = colorLevels, box=TRUE,
legend.width = 2, legend.mar = NULL,
legend.cex=cex,
axis.args=if(param=="Sign")
{list(labels=c("+","0","-"),at=c(1,0,-1),cex.axis=cex[3],mgp=mgp,hadj=1)
} else {
list(cex.axis=cex[3],mgp=mgp,hadj=1)
},
graphics.reset = TRUE, horizontal = FALSE, bigplot = bigplot, smallplot = smallplot,
legend.only = FALSE, col = colors,zlim = zlim)
if (!is.null(main))
title(main = main)
if(newplot){
par(old.par)
} else {
par(mfrow=cpar$mfrow)
cmfg=cpar$mfg
if(cmfg[2]<cmfg[4]){
cmfg[2] = cmfg[2]+1
} else {
if(cmfg[1]<cmfg[3]){
cmfg[2]=1
cmfg[1] = cmfg[1]+1
}
else{
cmfg[1]=1
cmfg[2]=1
}
}
par(mfg=cmfg,mar=cpar$mar)
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.