Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: testing equality of covariance matrices


From   May Boggess <mboggess@stata.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: testing equality of covariance matrices
Date   07 Oct 2004 10:12:31 -0500

On Wednesday, Dan asked:

> Does anyone know of an effective method for testing the equivalence of
> covariance matrices in Stata?
> 

For multivariate normal data, I would use Box's M-test to test for
equality of covariance matrices. Here is an example of how to perform
this test in Stata. I will use the turtle dataset from Mardia et. al.
"Multivariate Analysis", reproducing their example on page 141.

First let's input the data (it is only 48 lines):

clear
input str1 sex length width height
f 98 81 38
f 103 84 38
f 103 86 42
f 105 86 40
f 109 88 44
f 123 92 50
f 123 95 46
f 133 99 51
f 133 102 51
f 133 102 51
f 134 100 48
f 136 102 49
f 137 98 51
f 138 99 51
f 141 105 53
f 147 108 57
f 149 107 55
f 153 107 56
f 155 115 63
f 155 117 60
f 158 115 62
f 159 118 63
f 162 124 61
f 177 132 67
m 93 74 37
m 94 78 35
m 96 80 35
m 101 84 39
m 102 85 38
m 103 81 37
m 104 83 39
m 106 83 39
m 107 82 38
m 112 89 40
m 113 88 40
m 114 86 40
m 116 90 43
m 117 90 41
m 117 91 41
m 119 93 41
m 120 89 40
m 120 93 44
m 121 95 42
m 125 93 45
m 127 96 45
m 128 95 45
m 131 95 46
m 135 106 47
end

Now, I will summarize the data, just to verify I have close to the same
numbers as they do:

*----- look at the data -----------------------------------
 sum  length width height if sex=="m"
 sum  length width height if sex=="f"
 mat acc Sm = length width height if sex=="m", dev nocons
 local nm=r(N)
 mat Sm= (1/(`nm')) * Sm
 mat list Sm
 mat acc Sf = length width height if sex=="f", dev nocons
 local nf=r(N)
 mat Sf= (1/(`nf')) * Sf
 mat list Sf
 local n=`nm'+`nf'
 mat S = (`nm'*Sm+`nf'*Sf)*(1/`n')
 mat list S
 di det(Sm)
 di det(Sf)
 di det(S)

Since this test requires samples from multivariate normal distributions,
I should test this first. I am going to use two user-written commands
for this (these can be installed by using -findit- and then following
the instructions). 

*-------- tests for multivariate normality------------------


 omninorm length width height if sex=="m"
 preserve
 keep if sex=="m"
 multnorm length width height
 restore
 
 omninorm length width height if sex=="f"
 preserve
 keep if sex=="f"
 multnorm length width height
 restore
 
The turtle data passes the tests just fine, so now we calculate
Box's M using the formulas on page 140 in Mardia:

*-------- Box's M-test, page 140-141 Mardia ---------------------
 mat acc Sm = length width height if sex=="m", dev nocons
 local nm=r(N)
 mat Sm= (1/(`nm'-1)) * Sm
 mat acc Sf = length width height if sex=="f", dev nocons
 local nf=r(N)
 mat Sf= (1/(`nf'-1)) * Sf
 local n=`nm'+`nf'
 mat S = (`nm'*Sm+`nf'*Sf)*(1/`n')

 di det(Sm) 
 di det(Sf)
 di det(S)

 local k=2
 local p=3
 local gamma =  1-(2*`p'^2+3*`p'-1)/(6*(`p'+1)*(`k'-1))*
(1/(`nm'-1)+1/(`nf'-1)-1/(`n'-`k'))
 di "gamma =  `gamma' "

 mat SminvS=inv(Sm)*S
 local detm=det(SminvS)
 mat SfinvS=inv(Sf)*S
 local detf=det(SfinvS)
 local M = `gamma'*((`nm'-1)*ln(`detm')+(`nf'-1)*ln(`detf'))
 di "    M =  `M' "
 local df=(1/2)*`p'*(`p'+1)*(`k'-1)
 di "   df = `df' "
 di "p-value= " 1- chi2(`df', `M')

Note that I recalculated the covariance matrices because we use the
unbiased estimators to obtain M. Our results agree with those obtained
by Mardia, that the covariance matrices are significantly different for
male and female turtles.

--May
mmb@stata.com


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index