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
[email protected]
*
* 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/
