st: Calculating Mahalanobis distance

[lahiff@berkeley.edu]

I came across the postings from last spring.

estat mahal is useful to get Mahalanobis distances to group means after
discriminant analysis.

the stata code below may be helpful for those who need Mahalanobis
distances for a single group.

it's in the context of q-q-plot for 8 variables.

I've left in the output; you may want to trim it.

Maureen Lahiff
Division of Biostatistics
School of Public Health
University of California, Berkeley
lahiff@berkeley.edu


. summ  x1 - x8

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
          x1 |       150    99.78417    3.306238   91.35384   108.0844 x2
|       150     79.9418    3.386843   70.37208   87.93929 x3 |  
    150    100.0108     2.59557   92.85123   105.7917 x4 |      
150    49.95506    2.434154    43.8062   55.11219 x5 |       150
   154.6239    7.046493   136.5407   170.4589
-------------+--------------------------------------------------------
          x6 |       150    179.9501    9.416581   154.9588   202.6351 x7
|       150    75.11949    2.144913    69.8177   81.19572 x8 |  
    150    167.7471    5.586031   153.3667   182.8943

. gen x1cen = x1 - 99.78417

. gen x2cen = x2 - 79.9418

. gen x3cen = x3 - 100.0108

. gen x4cen = x4 - 49.95506

. gen x5cen = x5 - 154.6239

. gen x6cen = x6 - 179.9501

. gen x7cen = x7 - 75.11949

. gen x8cen = x8 - 167.7471

. summ x1cen - x8cen

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       x1cen |       150   -2.08e-06    3.306238  -8.430334   8.300241
x2cen |       150    4.54e-06    3.386843  -9.569722   7.997493
x3cen |       150   -3.39e-06     2.59557  -7.159573   5.780857
x4cen |       150   -3.83e-06    2.434154  -6.148858    5.15713
x5cen |       150    .0000253    7.046493  -18.08316   15.83501
-------------+--------------------------------------------------------
       x6cen |       150   -.0000144    9.416581  -24.99127     22.685
x7cen |       150   -3.51e-06    2.144913  -5.301795   6.076234
x8cen |       150    4.72e-06    5.586031  -14.38042   15.14717


. matrix accum SSCP = x1 - x8, nocons dev
(obs=150)

. matrix S = SSCP / ( r(N) - 1 )

. matrix list S

symmetric S[8,8]
            x1          x2          x3          x4          x5          x6
x1   10.931208
x2   8.0582002   11.470707
x3   7.2391807   6.6085483   6.7369822
x4  -1.4297405  -.45207946   -.2219191   5.9251058
x5  -4.5134066  -1.5641428  -1.7979711   15.222238   49.653068
x6   -5.162709  -1.7072405  -2.0620991   20.013281   62.783416   88.671992
x7    -.484696  -.23621137  -.14138931   .36964288   1.3110053   1.1950953
x8  -1.3335759  -.88431477  -.53214212   .85896683   3.1507676   2.6020501

            x7          x8
x7   4.6006533
x8   11.774218   31.203743


. matrix Sinv = invsym( S )


. mkmat x1cen - x8cen, matrix(Xcen)

. matrix Xcentran = Xcen'

. matrix mahaldist = Xcen * Sinv * Xcentran


. matrix distsvec = vecdiag(mahaldist)



. mat mahalanobis = distsvec'


. svmat mahalanobis, name(distances)




. matrix dir
    mahalanobis[150,1]
     distsvec[1,150]
    mahaldist[150,150]
     Xcentran[8,150]
         Xcen[150,8]
         Sinv[8,8]
            S[8,8]
         SSCP[8,8]



. qchi distances1, df(8) title(chi square plot of Mahalanobis distances)



*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/