 |
|
| The GLIMMIX Procedure |
Example 38.9 Testing Equality of Covariance and Correlation Matrices
Fisher’s iris data are widely used in multivariate statistics. They comprise measurements in millimeters of four flower attributes, the length and width of sepals and petals for 50 specimens from each of three species, Iris setosa, I. versicolor, and I. virginica (Fisher 1936).
When modeling multiple attributes from the same specimen, correlations among measurements from the same flower must be taken into account. Unstructured covariance matrices are common in this multivariate setting. Species comparisons can focus on comparisons of mean response, but comparisons of the variation and covariation are also of interest. In this example, the equivalence of covariance and correlation matrices among the species are examined.
The following DATA step creates the iris data in multivariate format—that is, each observation contains multiple response variables. The subsequent DATA step creates a data set in univariate form, where each observation corresponds to a single response variable. This is the form needed by the GLIMMIX procedure.
proc format;
value specname
1='Setosa '
2='Versicolor'
3='Virginica ';
value variable
1='Sepal Length'
2='Sepal Width '
3='Petal Length'
4='Petal Width';
run;
data iris;
input Y1 Y2 Y3 Y4 Species @@;
format Species specname.;
label Y1 = 'Sepal Length (mm)'
Y2 = 'Sepal Width (mm)'
Y3 = 'Petal Length (mm)'
Y4 = 'Petal Width (mm)';
datalines;
50 33 14 02 1 64 28 56 22 3 65 28 46 15 2 67 31 56 24 3
63 28 51 15 3 46 34 14 03 1 69 31 51 23 3 62 22 45 15 2
59 32 48 18 2 46 36 10 02 1 61 30 46 14 2 60 27 51 16 2
65 30 52 20 3 56 25 39 11 2 65 30 55 18 3 58 27 51 19 3
68 32 59 23 3 51 33 17 05 1 57 28 45 13 2 62 34 54 23 3
77 38 67 22 3 63 33 47 16 2 67 33 57 25 3 76 30 66 21 3
49 25 45 17 3 55 35 13 02 1 67 30 52 23 3 70 32 47 14 2
64 32 45 15 2 61 28 40 13 2 48 31 16 02 1 59 30 51 18 3
55 24 38 11 2 63 25 50 19 3 64 32 53 23 3 52 34 14 02 1
49 36 14 01 1 54 30 45 15 2 79 38 64 20 3 44 32 13 02 1
67 33 57 21 3 50 35 16 06 1 58 26 40 12 2 44 30 13 02 1
77 28 67 20 3 63 27 49 18 3 47 32 16 02 1 55 26 44 12 2
50 23 33 10 2 72 32 60 18 3 48 30 14 03 1 51 38 16 02 1
61 30 49 18 3 48 34 19 02 1 50 30 16 02 1 50 32 12 02 1
61 26 56 14 3 64 28 56 21 3 43 30 11 01 1 58 40 12 02 1
51 38 19 04 1 67 31 44 14 2 62 28 48 18 3 49 30 14 02 1
51 35 14 02 1 56 30 45 15 2 58 27 41 10 2 50 34 16 04 1
46 32 14 02 1 60 29 45 15 2 57 26 35 10 2 57 44 15 04 1
50 36 14 02 1 77 30 61 23 3 63 34 56 24 3 58 27 51 19 3
57 29 42 13 2 72 30 58 16 3 54 34 15 04 1 52 41 15 01 1
71 30 59 21 3 64 31 55 18 3 60 30 48 18 3 63 29 56 18 3
49 24 33 10 2 56 27 42 13 2 57 30 42 12 2 55 42 14 02 1
49 31 15 02 1 77 26 69 23 3 60 22 50 15 3 54 39 17 04 1
66 29 46 13 2 52 27 39 14 2 60 34 45 16 2 50 34 15 02 1
44 29 14 02 1 50 20 35 10 2 55 24 37 10 2 58 27 39 12 2
47 32 13 02 1 46 31 15 02 1 69 32 57 23 3 62 29 43 13 2
74 28 61 19 3 59 30 42 15 2 51 34 15 02 1 50 35 13 03 1
56 28 49 20 3 60 22 40 10 2 73 29 63 18 3 67 25 58 18 3
49 31 15 01 1 67 31 47 15 2 63 23 44 13 2 54 37 15 02 1
56 30 41 13 2 63 25 49 15 2 61 28 47 12 2 64 29 43 13 2
51 25 30 11 2 57 28 41 13 2 65 30 58 22 3 69 31 54 21 3
54 39 13 04 1 51 35 14 03 1 72 36 61 25 3 65 32 51 20 3
61 29 47 14 2 56 29 36 13 2 69 31 49 15 2 64 27 53 19 3
68 30 55 21 3 55 25 40 13 2 48 34 16 02 1 48 30 14 01 1
45 23 13 03 1 57 25 50 20 3 57 38 17 03 1 51 38 15 03 1
55 23 40 13 2 66 30 44 14 2 68 28 48 14 2 54 34 17 02 1
51 37 15 04 1 52 35 15 02 1 58 28 51 24 3 67 30 50 17 2
63 33 60 25 3 53 37 15 02 1
;
data iris_univ;
set iris;
retain id 0;
array y{4};
format var variable.;
id+1;
do var=1 to 4;
response = y{var};
output;
end;
drop y:;
run;
The following GLIMMIX statements fit a model with separate unstructured covariance matrices for each species:
ods select FitStatistics CovParms CovTests;
proc glimmix data=iris_univ;
class species var id;
model response = species*var;
random _residual_ / type=un group=species subject=id;
covtest homogeneity;
run;
The mean function is modeled as a cell-means model that allows for different means for each species and outcome variable. The covariances are modeled directly (R-side) rather than through random effects. The ID variable identifies the individual plant, so that responses from different plants are independent. The GROUP=SPECIES option varies the parameters of the unstructured covariance matrix by species. Hence, this model has 
covariance parameters: 
unique parameters for a 
covariance matrix for each of three species.
The COVTEST statement requests a test of homogeneity—that is, it tests whether varying the covariance parameters by the group effect provides a significantly better fit compared to a model in which different groups share the same parameter.
| Fit Statistics | |
|---|---|
| -2 Res Log Likelihood | 2812.89 |
| AIC (smaller is better) | 2872.89 |
| AICC (smaller is better) | 2876.23 |
| BIC (smaller is better) | 2963.21 |
| CAIC (smaller is better) | 2993.21 |
| HQIC (smaller is better) | 2909.58 |
| Generalized Chi-Square | 588.00 |
| Gener. Chi-Square / DF | 1.00 |
The "Fit Statistics" table shows the 
restricted (residual) log likelihood in the full model and other fit statistics (Output 38.9.1). The "-2 Res Log Likelihood" sets the benchmark against which a model with homogeneity constraint is compared. Output 38.9.2 displays the 30 covariance parameters in this model.
There appear to be substantial differences among the covariance parameters from different groups. For example, the residual variability of the petal length of the three species is 12.4249, 26.6433, and 40.4343, respectively. The homogeneity hypothesis restricts these variances to be equal and similarly for the other covariance parameters. The results from the COVTEST statement are shown in Output 38.9.3.
| Covariance Parameter Estimates | ||||
|---|---|---|---|---|
| Cov Parm | Subject | Group | Estimate | Standard Error |
| UN(1,1) | id | Species Setosa | 12.4249 | 2.5102 |
| UN(2,1) | id | Species Setosa | 9.9216 | 2.3775 |
| UN(2,2) | id | Species Setosa | 14.3690 | 2.9030 |
| UN(3,1) | id | Species Setosa | 1.6355 | 0.9052 |
| UN(3,2) | id | Species Setosa | 1.1698 | 0.9552 |
| UN(3,3) | id | Species Setosa | 3.0159 | 0.6093 |
| UN(4,1) | id | Species Setosa | 1.0331 | 0.5508 |
| UN(4,2) | id | Species Setosa | 0.9298 | 0.5859 |
| UN(4,3) | id | Species Setosa | 0.6069 | 0.2755 |
| UN(4,4) | id | Species Setosa | 1.1106 | 0.2244 |
| UN(1,1) | id | Species Versicolor | 26.6433 | 5.3828 |
| UN(2,1) | id | Species Versicolor | 8.5184 | 2.6144 |
| UN(2,2) | id | Species Versicolor | 9.8469 | 1.9894 |
| UN(3,1) | id | Species Versicolor | 18.2898 | 4.3398 |
| UN(3,2) | id | Species Versicolor | 8.2653 | 2.4149 |
| UN(3,3) | id | Species Versicolor | 22.0816 | 4.4612 |
| UN(4,1) | id | Species Versicolor | 5.5780 | 1.6617 |
| UN(4,2) | id | Species Versicolor | 4.1204 | 1.0641 |
| UN(4,3) | id | Species Versicolor | 7.3102 | 1.6891 |
| UN(4,4) | id | Species Versicolor | 3.9106 | 0.7901 |
| UN(1,1) | id | Species Virginica | 40.4343 | 8.1690 |
| UN(2,1) | id | Species Virginica | 9.3763 | 3.2213 |
| UN(2,2) | id | Species Virginica | 10.4004 | 2.1012 |
| UN(3,1) | id | Species Virginica | 30.3290 | 6.6262 |
| UN(3,2) | id | Species Virginica | 7.1380 | 2.7395 |
| UN(3,3) | id | Species Virginica | 30.4588 | 6.1536 |
| UN(4,1) | id | Species Virginica | 4.9094 | 2.5916 |
| UN(4,2) | id | Species Virginica | 4.7629 | 1.4367 |
| UN(4,3) | id | Species Virginica | 4.8824 | 2.2750 |
| UN(4,4) | id | Species Virginica | 7.5433 | 1.5240 |
Denote as 
the covariance matrix for species 
with elements 
. In processing the COVTEST hypothesis 
, the GLIMMIX procedure fits a model that satistfies the constraints
 |
 |
|||
 |
 |
|||
 |
 |
|||
 |
 |
|||
 |
 |
where is the covariance between the 
th and 
th variable for the 
th species. The 
2 restricted log likelihood of this restricted model is 2959.55 (Output 38.9.3). The change of 146.66 compared to the full model is highly significant. There is sufficient evidence to reject the notion of equal covariance matrices among the three iris species.
Equality of covariance matrices implies equality of correlation matrices, but the reverse is not true. Fewer constraints are needed to equate correlations because the diagonal entries of the covariance matrices are free to vary. In order to test the equality of the correlation matrices among the three species, you can parameterize the unstructured covariance matrix in terms of the correlations and use a COVTEST statement with general contrasts, as shown in the following statements:
ods select FitStatistics CovParms CovTests;
proc glimmix data=iris_univ;
class species var id;
model response = species*var;
random _residual_ / type=unr group=species subject=id;
covtest 'Equal Covariance Matrices' homogeneity;
covtest 'Equal Correlation Matrices' general
0 0 0 0 1 0 0 0 0 0
0 0 0 0 -1 0 0 0 0 0,
0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 -1 0 0 0 0 0,
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 -1 0 0 0 0,
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 -1 0 0 0 0,
0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 -1 0 0 0,
0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 -1 0 0 0,
0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 -1 0 0,
0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 -1 0 0,
0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 -1 0,
0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 -1 0,
0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 -1,
0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 -1 / estimates;
run;
The TYPE=UNR structure is a reparameterization of TYPE=UN. The models provide the same fit, as seen by comparison of the "Fit Statistics" tables in Output 38.9.1 and Output 38.9.4. The covariance parameters are ordered differently, however. In each group, the four variances precede the six correlations (Output 38.9.4). The first COVTEST statement tests the homogeneity hypothesis in terms of the UNR parameterization, and the result is identical to the test in Output 38.9.3. The second COVTEST statement restricts the correlations to be equal across groups. If 
is the correlation between the th and th variable for the th species, the 12 restrictions are
 |
 |
|||
 |
 |
|||
 |
 |
|||
 |
 |
|||
 |
 |
|||
 |
 |
The ESTIMATES option in the COVTEST statement requests that the GLIMMIX procedure display the covariance parameter estimates in the restricted model (Output 38.9.4).
| Fit Statistics | |
|---|---|
| -2 Res Log Likelihood | 2812.89 |
| AIC (smaller is better) | 2872.89 |
| AICC (smaller is better) | 2876.23 |
| BIC (smaller is better) | 2963.21 |
| CAIC (smaller is better) | 2993.21 |
| HQIC (smaller is better) | 2909.58 |
| Generalized Chi-Square | 588.00 |
| Gener. Chi-Square / DF | 1.00 |
| Covariance Parameter Estimates | ||||
|---|---|---|---|---|
| Cov Parm | Subject | Group | Estimate | Standard Error |
| Var(1) | id | Species Setosa | 12.4249 | 2.5102 |
| Var(2) | id | Species Setosa | 14.3690 | 2.9030 |
| Var(3) | id | Species Setosa | 3.0159 | 0.6093 |
| Var(4) | id | Species Setosa | 1.1106 | 0.2244 |
| Corr(2,1) | id | Species Setosa | 0.7425 | 0.06409 |
| Corr(3,1) | id | Species Setosa | 0.2672 | 0.1327 |
| Corr(3,2) | id | Species Setosa | 0.1777 | 0.1383 |
| Corr(4,1) | id | Species Setosa | 0.2781 | 0.1318 |
| Corr(4,2) | id | Species Setosa | 0.2328 | 0.1351 |
| Corr(4,3) | id | Species Setosa | 0.3316 | 0.1271 |
| Var(1) | id | Species Versicolor | 26.6433 | 5.3828 |
| Var(2) | id | Species Versicolor | 9.8469 | 1.9894 |
| Var(3) | id | Species Versicolor | 22.0816 | 4.4612 |
| Var(4) | id | Species Versicolor | 3.9106 | 0.7901 |
| Corr(2,1) | id | Species Versicolor | 0.5259 | 0.1033 |
| Corr(3,1) | id | Species Versicolor | 0.7540 | 0.06163 |
| Corr(3,2) | id | Species Versicolor | 0.5605 | 0.09797 |
| Corr(4,1) | id | Species Versicolor | 0.5465 | 0.1002 |
| Corr(4,2) | id | Species Versicolor | 0.6640 | 0.07987 |
| Corr(4,3) | id | Species Versicolor | 0.7867 | 0.05445 |
| Var(1) | id | Species Virginica | 40.4343 | 8.1690 |
| Var(2) | id | Species Virginica | 10.4004 | 2.1012 |
| Var(3) | id | Species Virginica | 30.4588 | 6.1536 |
| Var(4) | id | Species Virginica | 7.5433 | 1.5240 |
| Corr(2,1) | id | Species Virginica | 0.4572 | 0.1130 |
| Corr(3,1) | id | Species Virginica | 0.8642 | 0.03616 |
| Corr(3,2) | id | Species Virginica | 0.4010 | 0.1199 |
| Corr(4,1) | id | Species Virginica | 0.2811 | 0.1316 |
| Corr(4,2) | id | Species Virginica | 0.5377 | 0.1015 |
| Corr(4,3) | id | Species Virginica | 0.3221 | 0.1280 |
| Tests of Covariance Parameters Based on the Restricted Likelihood |
|||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Label | DF | -2 Res Log Like | ChiSq | Pr > ChiSq | Estimates H0 | Note | |||||||||||||||||||||||||||||
| Est1 | Est2 | Est3 | Est4 | Est5 | Est6 | Est7 | Est8 | Est9 | Est10 | Est11 | Est12 | Est13 | Est14 | Est15 | Est16 | Est17 | Est18 | Est19 | Est20 | Est21 | Est22 | Est23 | Est24 | Est25 | Est26 | Est27 | Est28 | Est29 | Est30 | ||||||
| Equal Covariance Matrices | 20 | 2959.55 | 146.66 | <.0001 | 26.5 | 11.5 | 18.52 | 4.19 | 0.530 | 0.756 | 0.378 | 0.365 | 0.471 | 0.484 | 26.5 | 11.54 | 18.5 | 4.19 | 0.530 | 0.756 | 0.378 | 0.365 | 0.471 | 0.484 | 26.5 | 11.5 | 18.5 | 4.19 | 0.530 | 0.756 | 0.378 | 0.365 | 0.471 | 0.484 | DF |
| Equal Correlation Matrices | 12 | 2876.38 | 63.49 | <.0001 | 16.5 | 14.9 | 4.84 | 1.44 | 0.561 | 0.683 | 0.402 | 0.384 | 0.498 | 0.522 | 24.4 | 9.16 | 17.4 | 3.00 | 0.561 | 0.683 | 0.402 | 0.384 | 0.498 | 0.522 | 35.1 | 10.8 | 27.4 | 8.14 | 0.561 | 0.683 | 0.402 | 0.384 | 0.498 | 0.522 | DF |
The result of the homogeneity test is identical to that in Output 38.9.3. The hypothesis of equality of the correlation matrices is also rejected with a chi-square value of 
and a p-value of 
. Notice, however, that the chi-square statistic is smaller than in the test of homogeneity due to the smaller number of restrictions imposed on the full model. The estimate of the common correlation matrix in the restricted model is
 |
|
|
Copyright © 2009 by SAS Institute Inc., Cary, NC, USA. All rights reserved.