#' Plot regions of the representative tree in 2D
#'
#' This function visualizes the regions of the representative tree
#' of the output of the \code{\link{mrs}} function.
#'
#' @param ans An \code{mrs} object.
#' @param type Different options on how to visualize the rectangular regions.
#' The options are \code{type = c("eff", "prob", "empty", "none")}.
#' Default is \code{type = "prob"}.
#' @param data.points Different options on how to plot the data points.
#' The options are \code{data.points = c("all", "differential", "none")}.
#' Default is \code{data.points = "all"}.
#' @param background Different options on the background.
#' The options are \code{background = c("smeared", "none") }.
#' @param group If \code{type = "eff"}, which group effect size is used.
#' Default is \code{group = 1}.
#' @param dim If the data are multivariate,
#' \code{dim} are the two dimensions plotted. Default is \code{dim = c(1,2)}.
#' @param levels Vector with the level of the regions to plot.
#' The default is to plot regions at all levels.
#' @param regions Binary vector indicating the regions to plot.
#' The default is to plot all regions.
#' @param legend Color legend for type. Default is \code{legend = FALSE}.
#' @param main Overall title for the legend.
#' @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
#' @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)
#' plot2D(ans, type = "prob", legend = TRUE)
#'
#' plot2D(ans, type="empty", data.points = "differential",
#' background = "none")
#'
#' plot2D(ans, type="none", data.points = "differential",
#' background = "smeared", levels = 4)
plot2D = function( ans,
type = "prob",
data.points = "all",
background = "none",
group = 1,
dim = c(1,2),
levels = sort(unique(ans$RepresentativeTree$Levels)),
regions = rep(1,length(ans$RepresentativeTree$Levels)),
legend = FALSE,
main = "default",
abs = TRUE )
{
if(class(ans)!="mrs")
{
print("ERROR: ans should be an mrs object.")
return(0)
}
if(ncol(ans$Data$X)<2)
{
print("ERROR: p>=2")
return(0);
}
.pardefault <- par(no.readonly = T)
size.tmp = ceiling(sqrt(length(levels)))
if( size.tmp * (size.tmp -1) == length(levels) )
{
mat = matrix( seq(1, size.tmp*(size.tmp-1)) , ncol=size.tmp, byrow=TRUE )
if(!legend)
layout( mat = mat )
else
{
mat = cbind(mat, size.tmp*(size.tmp - 1) + seq(1, size.tmp-1) )
layout( mat = mat, widths= c(rep(0.85, ncol(mat)-1)/(ncol(mat)-1), 0.15) )
}
}
else if( size.tmp * (size.tmp -1) > length(levels) )
{
mat = matrix( seq(1, size.tmp*(size.tmp-1)) , ncol=size.tmp, byrow=TRUE )
layout( mat = mat )
}
else if( size.tmp* size.tmp == length(levels) )
{
mat = matrix( seq(1, size.tmp^2), ncol=size.tmp, byrow=TRUE )
if(!legend)
layout( mat )
else
{
mat = cbind(mat, size.tmp^2 + seq(1, size.tmp) )
layout( mat = mat, widths= c(rep(0.85, ncol(mat)-1)/(ncol(mat)-1), 0.15) )
}
}
else
{
mat = matrix( seq(1, size.tmp*size.tmp) , ncol=size.tmp, byrow=TRUE )
layout(mat = mat)
# if(!legend)
# layout( mat = mat )
# else
# {
# mat = cbind(mat, size.tmp*(size.tmp - 1) + seq(1, size.tmp-1) )
# layout( mat = mat )
# }
}
par(mar = c(3.1, 3.1, 3.1, 1.1))
if(type == "prob")
{
names = round(ans$RepresentativeTree$AltProbs, digits=2)
col_range <- colorRampPalette(c("white","darkblue"))(100)
col = col_range[ ceiling( names*99 + 1) ]
if(main == "default")
main = "PMAPs"
}
else if(type == "eff")
{
if (abs == TRUE) {
names = abs(ans$RepresentativeTree$EffectSizes[,group])
col_range <- colorRampPalette(c("white","darkred"))(100)
col = col_range[ ceiling( names/max(names + 0.01)*99 + 1) ]
} else {
names = ans$RepresentativeTree$EffectSizes[, group]
col_range <- c(colorRampPalette(c("red","white"))(50),
colorRampPalette(c("white","dodgerblue"))(50))
col = col_range[ ceiling(names/max(abs(names)+0.01)*100/2+50) ]
}
if(main == "default")
main = paste("Eff. Size \n Group", group)
}
else
{
names = rep(1, length(ans$RepresentativeTree$AltProbs))
col = rep(NA, length(ans$RepresentativeTree$AltProbs))
if(main == "default")
main = ""
}
for( i in levels )
{
if(background == "smeared")
{
Lab.palette <- colorRampPalette(c("white", "black"), space = "rgb")
smoothScatter(ans$Data$X[,dim], nrpoints = 0, colramp = Lab.palette,
xlim = ans$Data$Omega[dim[1],], ylim=ans$Data$Omega[dim[2],],
main = paste("level",substitute(l, list(l = i))),
xlab = colnames(ans$Data$X)[dim[1]], ylab = colnames(ans$Data$X)[dim[2]])
}
else
{
plot(1, type="n", xlab = colnames(ans$Data$X)[dim[1]], ylab = colnames(ans$Data$X)[dim[2]],
xlim = ans$Data$Omega[dim[1],], ylim=ans$Data$Omega[dim[2],],
main = paste("level",substitute(l, list(l = i))) )
}
nodes = which(ans$RepresentativeTree$Levels==i & regions==1)
if( type =="eff")
nodes = nodes[sort.int( abs(ans$RepresentativeTree$EffectSizes[nodes,group]),
index.return= TRUE )$ix]
else if(type=="prob")
nodes = nodes[sort.int( ans$RepresentativeTree$AltProbs[nodes], index.return= TRUE )$ix]
for( j in nodes )
{
region = ans$RepresentativeTree$Regions[j,]
points_idx = ans$RepresentativeTree$DataPoints[[j]]
if(type=="empty")
rect(xleft=region[2*dim[1]-1] , ybottom=region[2*dim[2]-1] ,
xright=region[2*dim[1]] , ytop=region[2*dim[2]], col=NA, border = 1 )
else if(type!="none")
rect(xleft=region[2*dim[1]-1] , ybottom=region[2*dim[2]-1] ,
xright=region[2*dim[1]] , ytop=region[2*dim[2]], col=col[j], border = col[j] )
if(data.points=="all")
points(ans$Data$X[,dim], col=(ans$Data$G+1), pch = ans$Data$G )
else if(data.points=="differential")
{
points(ans$Data$X[points_idx,dim], col=(ans$Data$G[points_idx,]+1),
pch = ans$Data$G[points_idx,] )
}
}
if(legend && (i == levels[1]) && (data.points!="none") )
{
group_lab = rep(NA,max(ans$Data$G))
for(gg in 1:max(ans$Data$G))
group_lab[gg] = paste("sample", gg)
legend("bottomright", legend=group_lab, col=seq(2,max(ans$Data$G)+1),
pch = seq(1,max(ans$Data$G)), bg = "white" )
}
}
if(legend)
{
plot(NA,type="n",ann=FALSE,xlim=c(1,2),ylim=c(1,2),xaxt="n",yaxt="n",bty="n")
title(main, adj = 0)
if( (type == "prob") || (type == "eff") )
{
# mtext(main, side = 3, at = 1, cex = 1)
rect(xleft = 1,
ybottom = head(seq(1,2,1/100),-1),
xright = 1.2,
ytop = tail(seq(1,2,1/100),-1),
col=col_range, border=col_range )
rect(1, 1, 1.2, 2.0)
}
if(type == "prob")
mtext( format( round(seq(0,1,by=0.2), digits=1), nsmall=1) ,side=2,at=seq(1,2,by=.2),las=2,cex=1, line=0)
if(type == "eff") {
if (abs == TRUE) {
mtext( format(round(seq( 0, max(abs(names)) , length=5), digits=1), nsmall=1),
side=2,at=seq(1,2,length=5),las=2,cex=1, line=0)
} else {
mtext( format(round(seq(-max(abs(names)), max(abs(names)) , length=5), digits=1), nsmall=1),
side=2,at=seq(1,2,length=5),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.