#' Plot nodes of the representative tree
#'
#' This function visualizes the representative tree of the output of the \code{\link{mrs}} function.
#' For each node of the representative tree, the posterior probability of difference (PMAP) or the effect size is plotted.
#' Each node in the tree is associated to a region of the sample space.
#' All non-terminal nodes have two children nodes obtained by partitiing the parent region with a dyadic cut along a given direction.
#' The numbers under the vertices represent the cutting direction.
#'
#' @param ans A \code{mrs} object.
#' @param type What is represented at each node.
#' The options are \code{type = c("eff", "prob")}.
#' @param group If \code{type = "eff"}, which group effect size is used.
#' @param legend Color legend for type. Default is \code{legend = FALSE}.
#' @param main Main title. Default is \code{main = ""}.
#' @param node.size Size of the nodes. Default is \code{node.size = 5}.
#' @param abs If \code{TRUE}, plot the absolute value of the effect size.
#' Only used when \code{type = "eff"}.
#' @references Soriano J. and Ma L. (2016).
#' Probabilistic multi-resolution scanning for two-sample differences.
#' \emph{Journal of the Royal Statistical Society: Series B (Statistical Methodology)}.
#' \url{http://onlinelibrary.wiley.com/doi/10.1111/rssb.12180/abstract}
#' @references Ma L. and Soriano J. (2016).
#' Analysis of distributional variation through multi-scale Beta-Binomial modeling.
#' \emph{arXiv}.
#' \url{http://arxiv.org/abs/1604.01443}
#' @export
#' @note The package \pkg{igraph} is required.
#' @examples
#' set.seed(1)
#' p = 2
#' n1 = 200
#' n2 = 200
#' mu1 = matrix( c(9,9,0,4,-2,-10,3,6,6,-10), nrow = 5, byrow=TRUE)
#' mu2 = mu1; mu2[2,] = mu1[2,] + 1
#'
#' Z1 = sample(5, n1, replace=TRUE)
#' Z2 = sample(5, n2, replace=TRUE)
#' X1 = mu1[Z1,] + matrix(rnorm(n1*p), ncol=p)
#' X2 = mu2[Z2,] + matrix(rnorm(n2*p), ncol=p)
#' X = rbind(X1, X2)
#' colnames(X) = c(1,2)
#' G = c(rep(1, n1), rep(2,n2))
#'
#' ans = mrs(X, G, K=8)
#' plotTree(ans, type = "prob", legend = TRUE)
plotTree <- function(ans, type="prob", group = 1, legend = FALSE, main = "", node.size=5, abs = TRUE)
{
if(class(ans)!="mrs")
{
print("ERROR: ans should be an mrs object.")
return(0)
}
.pardefault <- par(no.readonly = T)
if(legend)
{
layout(matrix(c(1,2), nrow=1),widths=c(.85,.15))
par(mar=c(1.1, 1.1, 3.1, 0.1))
}
else
{
layout(1)
}
if(type == "prob")
{
box.size = (ans$RepresentativeTree$AltProbs)*node.size + .2*node.size
name = round(ans$RepresentativeTree$AltProbs, digits=2)
col_range <- colorRampPalette(c("white","darkblue"))(100)
col = col_range[ ceiling( name/max(name + 0.01)*99 + 1) ]
}
else if(type == "eff")
{
box.size = abs(ans$RepresentativeTree$EffectSizes[,group])/max(abs(ans$RepresentativeTree$EffectSizes[,group]))*node.size + .2*node.size
if (abs == TRUE) {
name = abs(ans$RepresentativeTree$EffectSizes[,group])
col_range <- colorRampPalette(c("white","darkred"))(100)
col = col_range[ ceiling( name/max(name + 0.01)*99 + 1) ]
} else {
name = ans$RepresentativeTree$EffectSizes[, group]
col_range <- c(colorRampPalette(c("red","white"))(50),
colorRampPalette(c("white","dodgerblue"))(50))
col = col_range[ ceiling(name/max(abs(name)+0.01)*100/2+50) ]
}
}
else
{
box.size = 5
name = NA
col = "white"
}
M = list()
it = 0
for(i in 1:length(ans$RepresentativeTree$Levels))
{
for(j in 1:length(ans$RepresentativeTree$Levels))
{
if(ans$RepresentativeTree$Levels[i] == (ans$RepresentativeTree$Levels[j] -1)
&& ( 2*ans$RepresentativeTree$Ids[i] == ans$RepresentativeTree$Ids[j]
|| (2*ans$RepresentativeTree$Ids[i]+1) == ans$RepresentativeTree$Ids[j] ))
{
it = it + 1
M[[it]] = c(i,j)
}
}
}
M = unlist(M)
vertex.label = ans$RepresentativeTree$Directions
vertex.label[vertex.label==0] = NA
vertex.label[ans$RepresentativeTree$Levels==max(ans$RepresentativeTree$Levels)]= NA
G = igraph::graph(M, directed=FALSE)
co <- igraph::layout.reingold.tilford(G, params=list(root=which(ans$RepresentativeTree$Levels==0)))
igraph::plot.igraph(G, layout=co,
vertex.size = box.size, edge.label=NA,
vertex.color = col, vertex.label.color = "black",
vertex.label = vertex.label, asp=0, main=main, vertex.label.family = "Helvetica",
vertex.label.dist = 1.5, vertex.label.degree = pi/2)
if(legend)
{
plot(NA,type="n",ann=FALSE,xlim=c(1,2),ylim=c(1,2),xaxt="n",yaxt="n",bty="n")
rect(xleft = 1.25,
ybottom = head(seq(1,2,1/100),-1),
xright = 1.75,
ytop = tail(seq(1,2,1/100),-1),
col=col_range, border = col_range )
rect(1.25, 1, 1.75, 2.0)
if(type == "prob")
mtext(formatC(seq(0,1,.1), format = "f", digits = 1),side=2,at=seq(1,2.,by=.1),las=2,cex=1, line=0)
if(type == "eff") {
if (abs == TRUE) {
mtext(formatC(seq( 0, max(name) , length=11), format = "f", digits = 1),
side=2,at=seq(1,2.,by=.1),las=2,cex=1, line=0)
} else {
mtext(formatC(seq( -max(abs(name)), max(abs(name)) , length=11), format = "f", digits = 1),
side=2,at=seq(1,2.,by=.1),las=2,cex=1, line=0)
}
}
}
par(.pardefault)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.