# Pairwise Mahalanobis distance in R

I'm trying to calculate a Mahalanobis-type pairwise distance matrix in R. I have 33 individuals, each with 10 variables. The idea is to get a distance matrix D, where

$$D_{i,j}=(\mathbf{X}_i-\mathbf{X}_j)W^{-1}(\mathbf{X}_i-\mathbf{X}_j)^T$$

However I haven't been able build proper code for it.

-
Type ?mahalanobis in R and look the documentation. –  MYaseen208 Jul 19 '12 at 21:56
@MYaseen208: Wrong. Look below. –  user603 Aug 2 '12 at 8:32
–  rpierce Jul 27 at 21:56

## migrated from stackoverflow.comAug 2 '12 at 0:23

This question came from our site for professional and enthusiast programmers.

Here is the code to do it:

library("MASS")
library("ICSNP")

x0<-mvrnorm(33,1:10,diag(c(seq(1,1/2,l=10)),10))
x1<-pair.diff(x0) #C-implementation.
dM<-mahalanobis(x1,colMeans(x1),var(x1))


Following Roman Luštrik's suggestion, here are more details. The OP asked for pairwise Mahalanobis distance, which are multivariate U-statistics of distance. I have first seen them mentionned in Croux et al. 94 (below equation 6.4) but i'm sure others such as Oja have explored this concept.

-
It would be helpful if you expanded your answer with at least a short comment on your code. Why do you use it that way, why this approach is better compared to another approach... –  Roman Luštrik Aug 2 '12 at 9:13
It would be helpful if you expanded your comment with at least a short description of what that "another approach" could be. –  user603 Aug 2 '12 at 9:26
I leave that at your discretion. My point was that it would help others (including me) if you gave your answer some more context. –  Roman Luštrik Aug 2 '12 at 13:52
Here is a question from @user21060 who cannot leave comment at the moment: Should var(x1) be var(x0) instead as the $\sigma$ should be the variance of the original data, and I also question whether the empirical mean should be used here? –  chl May 31 at 11:17
thanks to user21060 for the questions. The answers are no and yes. both should be for x1. I'm correcting the oversight now. –  user603 May 31 at 11:19

The following worked for me in similar example where R is a dataframe of 54 individuals and 8 variables. Mahalanobis distance Ma between individuals X1 and X2 can be computed as ff:

# express difference (X1-X2) as atomic row vector
d <- as.matrix(X1-X2)[1,]

# find inverse of covariance matrix
Sx <- solve(cov(R))

# Mahalanobis calculation forced in two steps
Ma <- (d %\*% Sx) %\*% d

-
what are X1 and X2?...can you make this self contained? –  user603 Jan 26 at 21:48

You could try the gendistance function in the nbpMatching package

Here's a short example modified from the help page, with two variables instead of 10:

df <- data.frame(id=1:33, val1=rnorm(33), val2=rnorm(33))
df.dist <- gendistance(df, idcol=1)
df.dist\$dist


The distance matrix will have a 34th row/column-- this is for use in matching, and you can ignore it.

-
thanks for the extra info. From what i understand of the OP's description the answer should be a vector with 33*32/2 positive numbers in it... – –  user603 Jan 30 at 16:12

I have posted some theory giving the rationale behind the following code at Pairwise Mahalanobis distance

fastPwMahal = function(x1,invCovMat) {
SQRT = with(svd(invCovMat), u %*% diag(d^0.5) %*% t(v))
dist(x1 %*% SQRT)
}


At the link above I have also shown that the (currently top-voted) solution using the ICSNP package on this forum seems to be incorrect.

I joined stack exchange to answer this question, but initially did not have enough reputation to do so. I'm hoping a mod can arrange for these threads to be merged. Also, we need clarification or a disclaimer about the ICSNP solution on this page -- it's risky that the top-voted solution is essentially unchallenged and may be leading users to incorrect distance calculations.

-