[R] Re: Where is a function for minimum spanning tree in R for win

[Yvonnick Noël]

I have written the following code for my own needs last year. You need to load the MVA (the DIST function is used)and SCATTERPLOT3D (for plotting high dimensional data) packages first.

For an example, try something like :

X<-matrix(runif(60),20,3)
N<-mst(X, plot.true=T)

A graph symmetric matrix is returned in N which contains 1s for maximum links and 0 elsewhere. Note that, if your data are p-dimensional, with p>3, the graph is plotted in the space of the 3 first principal components. In my practice, this very simple visualization approach often yields insightful pictures of the data structure.

Hope you'll find it useful. 

Yvonnick Noel, PhD.
U. of Lille 3
France

===
mst <- function(X,plot.true=F)
{
    size <- dim(X)
    n <- size[1]
    p <- size[2]

    D <- as.matrix(dist(X))             # distances matrix
    N <- matrix(0,n,n)                  # graph matrix
    tree <- NULL                        # sequence of nodes added to the graph

    large.value <- max(D) + 1           # set diagonal to large values (for min detection)
    diag(D) <- large.value

    index.i <- 1                        # construct graph
    for(i in 1:(n-1))
    {
        tree <- c(tree,index.i)
        m <- apply(as.matrix(D[,tree]),2,min)
        a <- sort.index(D[,tree])[1,]
        b <- sort.index(m)[1]
        index.j <- tree[b]
        index.i <- a[b]

        N[index.i,index.j] <- 1
        N[index.j,index.i] <- 1

        for(j in tree)
        {
            D[index.i,j] <- large.value
            D[j,index.i] <- large.value
        }
    }
    
    # Plots graph if required
    if(plot.true && (p>1))
    {
      if(p==2) 
      {
        plot(X[,1],X[,2],main="Minimum Spanning Tree",xlab="",ylab="",pch=19,col="blue")
        for(k in 2:n)
        {
            for(l in 1:(k-1))
            {
                if(N[k,l]==1) lines(rbind(X[k,],X[l,]),pch=19,col="blue")
            }
        }
      }

      else if(p==3)
      {
        scat.res <- scatterplot3d(X,highlight.3d=T)
        for(k in 2:n)
        {
            for(l in 1:(k-1))
            {
                if(N[k,l]==1) scat.res$points3d(rbind(X[k,],X[l,]),type="l",pch=19,col="blue")
            }
        }
      }

      else
      {
        C <- prcomp(X,retx=T)$x[,1:3]
        scat.res <- scatterplot3d(C,highlight.3d=T)
        for(k in 2:n)
        {
            for(l in 1:(k-1))
            {
                if(N[k,l]==1) scat.res$points3d(rbind(C[k,],C[l,]),type="l",pch=19,col="blue")
            }    
        }
    
      }
    }
    
    return(N)

}
sort.index <- function(X)
{
    # Function returning an index matrix for an increasing sort
    if(length(X)==1) return(1)                  # sorting a scalar?
    if(!is.matrix(X)) X <- as.matrix(X)         # force vector into matrix
    
    n <- nrow(X)

    # find the permutation
    #--> Rank transforming with ties breaking
    #rk <- apply(X,2,function(v) order(v,1:n))
    rk <- apply(X,2,rank)
    
    #--> Match()ing an integer sequence with the rank-transformed matrix
    apply(rk, 2, function(v) match(1:n,v))  
}

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._