From walton.green at yale.edu  Sun Oct 23 02:12:52 2005
From: walton.green at yale.edu (Walton A. Green)
Date: Sat, 22 Oct 2005 20:12:52 -0400 (EDT)
Subject: [Rd] Coloring leaves, twigs and labels in plot.dendrogram()
Message-ID: <Pine.LNX.4.44.0510222005200.25034-100000@ares.its.yale.edu>


Core developers,

I couldn't find any simple way to send a vector of colors to apply to each
terminal in plot.dendrogram() or plot.hclust()---I asked R-help about it a
few weeks ago and didn't get any response---so I hacked that functionality
into the plot.dendrgram code (see below for hacked function plus 
examples)....
    
Is there any chance this functionality could be added to the base code? 
Unquestionably it can be done more elegantly than below, but surely it 
would be a useful feature to have....

Or is there an easier way of doing this that I have missed?

Walton Green

*	*	*

my.dend.plot <- function (x, type = c("rectangle", "triangle"),
    center = FALSE, edge.root = is.leaf(x) || !is.null(attr(x,
    "edgetext")), nodePar = NULL, edgePar = list(),
    leaflab = c("perpendicular", "textlike", "none"), dLeaf = NULL,
    xlab = "", ylab = "", xaxt = "n", yaxt = "s", horiz = FALSE,
    frame.plot = FALSE, twig.col = NULL, lab.col = NULL, leaf.col = NULL, 
...){

############################## 80 Characters Wide 
##############################

plotNode <- function (x1, x2, subtree, type, center, leaflab, dLeaf, 
nodePar, 
    edgePar, horiz = FALSE, leafcounter = leafcounter, ...) 
{
	#cat('.') uncomment this for debugging...
    inner <- !is.leaf(subtree) && x1 != x2
    yTop <- attr(subtree, "height")
    bx <- plotNodeLimit(x1, x2, subtree, center)
    xTop <- bx$x
    hasP <- !is.null(nPar <- attr(subtree, "nodePar"))
    if (!hasP) 
        nPar <- nodePar
    if (getOption("verbose")) {
        cat(if (inner) 
            "inner node"
        else "leaf", ":")
        if (!is.null(nPar)) {
            cat(" with node pars\n")
            str(nPar)
        }
        cat(if (inner) 
            paste(" height", formatC(yTop), "; "), "(x1,x2)= (", 
            formatC(x1, wid = 4), ",", formatC(x2, wid = 4), 
            ")", "--> xTop=", formatC(xTop, wid = 8), "\n", sep = "")
    }
    Xtract <- function(nam, L, default, indx) rep(if (nam %in% 
        names(L)) 
        L[[nam]]
    else default, length.out = indx)[indx]
    asTxt <- function(x) if (is.character(x) || is.expression(x) || 
        is.null(x)) 
        x
    else as.character(x)
    i <- if (inner || hasP) 
        1
    else 2
    if (!is.null(nPar)) {
        pch <- Xtract("pch", nPar, default = 1:2, i)
        cex <- Xtract("cex", nPar, default = c(1, 1), i)
        col <- Xtract("col", nPar, default = par("col"), i)
        bg <- Xtract("bg", nPar, default = par("bg"), i)
		if(is.null(leaf.col) || inner){
		    leaf.col <- rep(col, x2)
		}
		if(!(is.null(leaf.col)) && (is.null(lab.col)) && 
(is.null(twig.col))){
		    leafcounter <- leafcounter + 1
		}
        points(if (horiz) 
            cbind(yTop, xTop)
        else cbind(xTop, yTop), pch = pch, bg = bg, col = 
leaf.col[leafcounter - 1], 
            cex = cex)
    }
    if (leaflab == "textlike") 
	p.col <- Xtract("p.col", nPar, default = "white", i)
    lab.cex <- Xtract("lab.cex", nPar, default = c(1, 1), i)
    lab.font <- Xtract("lab.font", nPar, default = par("font"), i)
    if (is.leaf(subtree)) {
        if (leaflab == "perpendicular") {
            if (horiz) {
                X <- yTop + dLeaf * lab.cex
                Y <- xTop
                srt <- 0
                adj <- c(0, 0.5)
            }
            else {
                Y <- yTop - dLeaf * lab.cex
                X <- xTop
                srt <- 90
                adj <- 1
            }
            nodeText <- asTxt(attr(subtree, "label"))
			if(!is.null(lab.col)){
			    lab.col <- lab.col[leafcounter - 1]
				if(is.null(twig.col) && 
is.null(leaf.col)){
				    leafcounter <- leafcounter + 1
				}
            }
			text(X, Y, nodeText, xpd = TRUE, srt = srt, adj = 
adj, 
                cex = lab.cex, col = lab.col, font = lab.font)
        }
    }
    else if (inner) {
        segmentsHV <- function(x0, y0, x1, y1, ...) {
            if (horiz) 
                segments(y0, x0, y1, x1, ...)
            else segments(x0, y0, x1, y1, ...)
        }
        for (k in 1:length(subtree)) {
            child <- subtree[[k]]
            yBot <- attr(child, "height")
            if (getOption("verbose")) 
                cat("ch.", k, "@ h=", yBot, "; ")
            if (is.null(yBot)) 
                yBot <- 0
            xBot <- if (center) 
                mean(bx$limit[k:(k + 1)])
            else bx$limit[k] + .midDend(child)
            hasE <- !is.null(ePar <- attr(child, "edgePar"))
            if (!hasE) 
                ePar <- edgePar
            i <- if (!is.leaf(child) || hasE) 
                1
            else 2
            col <- Xtract("col", ePar, default = par("col"), 
                i)
            lty <- Xtract("lty", ePar, default = par("lty"), 
                i)
            lwd <- Xtract("lwd", ePar, default = par("lwd"), 
                i)
            if (type == "triangle") {
                segmentsHV(xTop, yTop, xBot, yBot)
            }
            else {
			    if(is.null(twig.col)){
                    segmentsHV(xTop, yTop, xBot, yTop)
                    segmentsHV(xBot, yTop, xBot, yBot)
				}else{
				    if(is.leaf(child)){
                        # cat(leafcounter) uncomment this for debugging...
					    segmentsHV(xTop, yTop, xBot, 
yTop)
						segmentsHV(xBot, yTop, 
xBot, yBot, col = twig.col[leafcounter])
						leafcounter <- leafcounter 
+ 1
					}else{
					    segmentsHV(xTop, yTop, xBot, 
yTop)
						segmentsHV(xBot, yTop, 
xBot, yBot)
					}
				}
            }
            vln <- NULL
            if (is.leaf(child) && leaflab == "textlike") {
                nodeText <- asTxt(attr(child, "label"))
                if (getOption("verbose")) 
                  cat("-- with \"label\"", format(nodeText))
                hln <- 0.6 * strwidth(nodeText, cex = lab.cex)/2
                vln <- 1.5 * strheight(nodeText, cex = lab.cex)/2
                rect(xBot - hln, yBot, xBot + hln, yBot + 2 * 
                  vln, col = p.col)
                text(xBot, yBot + vln, nodeText, xpd = TRUE, 
                  cex = lab.cex, col = lab.col, font = lab.font)
            }
            if (!is.null(attr(child, "edgetext"))) {
                edgeText <- asTxt(attr(child, "edgetext"))
                if (getOption("verbose")) 
                  cat("-- with \"edgetext\"", format(edgeText))
                if (!is.null(vln)) {
                  mx <- if (type == "triangle") 
                    (xTop + xBot + ((xTop - xBot)/(yTop - yBot)) * 
                      vln)/2
                  else xBot
                  my <- (yTop + yBot + 2 * vln)/2
                }
                else {
                  mx <- if (type == "triangle") 
                    (xTop + xBot)/2
                  else xBot
                  my <- (yTop + yBot)/2
                }
                p.col <- Xtract("p.col", ePar, default = "white", 
                  i)
                p.border <- Xtract("p.border", ePar, default = par("fg"), 
                  i)
                p.lwd <- Xtract("p.lwd", ePar, default = lwd, 
                  i)
                p.lty <- Xtract("p.lty", ePar, default = lty, 
                  i)
                t.col <- Xtract("t.col", ePar, default = col, 
                  i)
                t.cex <- Xtract("t.cex", ePar, default = 1, i)
                t.font <- Xtract("t.font", ePar, default = par("font"), 
                  i)
                vlm <- strheight(c(edgeText, "h"), cex = t.cex)/2
                hlm <- strwidth(c(edgeText, "m"), cex = t.cex)/2
                hl3 <- c(hlm[1], hlm[1] + hlm[2], hlm[1])
                polygon(mx + c(-hl3, hl3), my + sum(vlm) * c(-1:1, 
                  1:-1), col = p.col, border = p.border, lty = p.lty, 
                  lwd = p.lwd)
                text(mx, my, edgeText, cex = t.cex, col = t.col, 
                  font = t.font)
			}
	            leafcounter <- 1 + plotNode(bx$limit[k], bx$limit[k + 
1], subtree = child, 
                type, center, leaflab, dLeaf, nodePar, edgePar, 
                horiz, leafcounter = leafcounter, ...)
        } # End of loop through k
	}
return(leafcounter - 1)
}

############################## 80 Characters Wide 
##############################

plotNodeLimit <- function (x1, x2, subtree, center) 
{
    inner <- !is.leaf(subtree) && x1 != x2
    if (inner) {
        K <- length(subtree)
        mTop <- .memberDend(subtree)
        limit <- integer(K)
        xx1 <- x1
        for (k in 1:K) {
            m <- .memberDend(subtree[[k]])
            xx1 <- xx1 + (if (center) 
                (x2 - x1) * m/mTop
            else m)
            limit[k] <- xx1
        }
        limit <- c(x1, limit)
    }
    else {
        limit <- c(x1, x2)
    }
    mid <- attr(subtree, "midpoint")
    center <- center || (inner && !is.numeric(mid))
    x <- if (center) 
        mean(c(x1, x2))
    else x1 + (if (inner) 
        mid
    else 0)
    list(x = x, limit = limit)
}

############################## 80 Characters Wide 
##############################

.memberDend <- function (x) 
{
    r <- attr(x, "x.member")
    if (is.null(r)) {
        r <- attr(x, "members")
        if (is.null(r)) 
            r <- 1:1
    }
    r
}

.midDend <- function (x) 
if (is.null(mp <- attr(x, "midpoint"))) 0 else mp

############################## 80 Characters Wide 
##############################
########## Body follows

    type <- match.arg(type)
    leaflab <- match.arg(leaflab)
    hgt <- attr(x, "height")
    if (edge.root && is.logical(edge.root)) 
        edge.root <- 0.0625 * if (is.leaf(x)) 
            1
        else hgt
    mem.x <- .memberDend(x)
    yTop <- hgt + edge.root
    if (center) {
        x1 <- 0.5
        x2 <- mem.x + 0.5
    }
    else {
        x1 <- 1
        x2 <- mem.x
    }
    xlim <- c(x1 - 1/2, x2 + 1/2)
    ylim <- c(0, yTop)
    if (horiz) {
        xl <- xlim
        xlim <- rev(ylim)
        ylim <- xl
        tmp <- xaxt
        xaxt <- yaxt
        yaxt <- tmp
    }
    plot(0, xlim = xlim, ylim = ylim, type = "n", xlab = xlab, 
        ylab = ylab, xaxt = xaxt, yaxt = yaxt, frame.plot = frame.plot, 
        ...)
    if (is.null(dLeaf)) 
        dLeaf <- 0.75 * (if (horiz) 
            strwidth("w")
        else strheight("x"))
    if (edge.root) {
        x0 <- plotNodeLimit(x1, x2, x, center)$x
        if (horiz){
            segments(hgt, x0, yTop, x0)
        }else{
            segments(x0, hgt, x0, yTop)
        }
        if (!is.null(et <- attr(x, "edgetext"))) {
            my <- mean(hgt, yTop)
            if (horiz) 
                text(my, x0, et)
            else text(x0, my, et)
        }
    }
    plotNode(x1, x2, x, type = type, center = center, leaflab = leaflab, 
        dLeaf = dLeaf, nodePar = nodePar, edgePar = edgePar, 
        horiz = horiz, leafcounter = 1, ...)
} # End of function

############################## 80 Characters Wide 
##############################

## Examples
hc <- hclust(dist(USArrests), "ave")
NewEnglandStates <- rep(0, 50)
NewEnglandStates[c(7, 19, 21, 29, 39, 45)] <- 1
my.dend.plot(as.dendrogram(hc, hang = 0.02),
             twig.col = (NewEnglandStates + 1)[hc$order],
			 lab.col = (NewEnglandStates + 1)[hc$order])
my.dend.plot(as.dendrogram(hc), twig.col = (NewEnglandStates + 
1)[hc$order],
             leaf.col = (NewEnglandStates + 1)[hc$order],
			 nodePar = list(pch = c(NA,1)))



[Rd] Coloring leaves, twigs and labels in plot.dendrogram()
Walton A. Green