[R] fixing a problem in the subtree code

Wiener, Matthew matthew_wiener at merck.com
Mon Feb 11 16:28:10 CET 2002


Hi, all.  Last week I posted some code for extracting subtrees of trees in
hclust format.  Petra Steiner quickly found an example for which the code
breaks, and sent it to me.

The problem seems to be that I had not considered the possibility of trees
with unlabeled nodes.  In the new version of f.make.subtree (below), I steal
some code from plot.hclust to assign labels if there are none.  That fixes
the problem.  At the end I drop the labels again if the original tree didn't
have any.

Thanks to Petra for the help.  If anyone has other examples, or other
suggestions for how to improve the code, please let me know.

Regards,

Matthew Wiener
Applied Computer Science & Mathematics Dept.
Merck Research Labs
Rahway, NJ 07090
732-594-5303 

----------------------------------------------------------------------------
----

"f.make.subtree" <-
  function (tree, groups, n = 1, add.element = NULL) 
{
  which.nodes <- which(groups == n)
  names(which.nodes) <- NULL
  ## steal some code from plot.hclust
  use.labels <-
    if (is.null(tree$labels)) 
      paste(1:(nrow(tree$merge) + 1))
    else as.character(tree$labels)

  subtree.labels <- use.labels[which.nodes]
  
  subtree.order <- tree$order[use.labels[tree$order] %in% 
                              subtree.labels]
  which.merge.elements <- which(-tree$merge %in% which.nodes)
  which.merge.rows <- sort(unique(which.merge.elements%%nrow(tree$merge)))
  which.merge.rows[which.merge.rows == 0] <- nrow(tree$merge)
  old.length <- 0
  new.length <- length(which.merge.rows)
  if (length(which.merge.elements) == 1) {
    res <- list(merge = NULL, heights = tree$height[which.merge.rows], 
                order = 1,
                labels = {if(is.null(tree$labels)) NULL else
subtree.labels},
                method = tree$method, 
                call = tree$call, dist.method = tree$dist.method)
  }
  else {
    while (new.length - old.length > 0) {
      old.length <- new.length
      new.rows <- numeric(0)
      more.new.rows <- numeric(0)
      row.elements <- as.vector(tree$merge[which.merge.rows, 
                                           ])
      pos.elements <- row.elements[row.elements > 0]
      if (length(pos.elements) > 0) 
        new.rows <- pos.elements[apply(tree$merge[pos.elements, 
                                                  , drop = F], 1,
function(x) {
                                                    any(x > 0) & all(-x[x <
0] %in% which.nodes)
                                                  })]
      more.new.rows <- which(apply(tree$merge, 1, function(x) {
        all(x %in% which.merge.rows)
      }))
      which.merge.rows <- sort(unique(c(which.merge.rows, 
                                        new.rows, more.new.rows)))
      new.length <- length(which.merge.rows)
    }
    merge.list <- tree$merge[which.merge.rows, , drop = F]
    merge.height <- tree$height[which.merge.rows]
    pos.elements <- sort(merge.list[merge.list > 0])
    for (i in seq(along = pos.elements)) merge.list[merge.list == 
                    pos.elements[i]] <- which(pos.elements[i] ==
which.merge.rows)
    neg.elements <- merge.list[merge.list < 0]
    minus.neg.elements <- -neg.elements
    num.elements <- length(neg.elements)
    change.table <- cbind(sort(minus.neg.elements), 1:num.elements)
    for (i in 1:num.elements) {
      merge.list[merge.list == -change.table[i, 1]] <- -change.table[i, 
                                                                     2]
    }
    old.order <- subtree.order
    subtree.order <- match(use.labels[old.order], subtree.labels)
    res <- list(merge = merge.list, height = merge.height, 
                order = subtree.order,
                labels = {if(is.null(tree$labels)) NULL else
subtree.labels},
                method = tree$method, 
                call = list(match.call(), tree$call), dist.method =
tree$dist.method)
    class(res) <- "hclust"
    for (i in seq(along = add.element)) {
      res[[add.element[i]]] <- tree[[add.element[i]]][which.nodes]
    }
  }
  res
}

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list