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

Yvonnick Noël ynoel at nordnet.fr
Sun Dec 10 08:35:18 CET 2000

```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
p <- size

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

```