Section 10: The Unequal Variance ANOVA Model

1. Checking Equality of Variance in the ANOVA Model
2. Testing Means with the Unequal Variance Model
3. Homogeneity of Variance Tests
4. Weighted Least Squares
5. How can PROC MIXED handle this unequal variance situation?
6. Multiple Comparisons in the Unequal Variance ANOVA Model
7. Dispersion Effects Model


One of the three assumptions of an Analysis of Variance (ANOVA) model 
listed in Section 1 (found at http://www.uoregon.edu/~robinh/glm01_assumpt.txt )
 
  "The residuals have equal variances across groups".

This section examines the several methods to check this assumption and a 
proposed solution if indeed variances across groups should not be 
considered equal.
 

1. Checking Equality of Variance in the ANOVA Model

For an ANOVA model with one "fixed" factor and at least two observations 
for each level of that factor, the MEANS statement from PROC GLM allows 
you to request several types of homogeneity of variance tests (HOV). These 
tests help you determine if one of the primary assumptions of ANOVA models 
-- equal variances across groups -- can be assumed true.

Heterogeneity of variance tests listed here are defined only for one-way 
ANOVAs. With k factors (k>1) it is possible to re-code the data from K-way 
ANOVA so that there is a single variable that acts as a surrogate for all 
K effects and then compute heterogeneity of variance tests on this recoded 
data. However, the greater the complexity of the design, the further you 
may diverge from the underlying assumptions of the test.

To illustrate with two fixed factors, (A and B), you can recode them into 
one variable that uniquely defines all possible combinations of levels of 
the two factors. First, sort the factors by their respective levels, with 
factor B the inner variable changing levels for each level of A. If factor 
A has 2 levels and factor B has 3 levels combine them into one 6-level 
factor which defines all 6 possible combinations:

A   B  class_var
1   1     1
1   2     2
1   3     3
2   1     4
2   2     5
2   3     6

If the factors are numerically coded as above, a simple math 
transformation is evident:

class_var = 3*(A-1) + B

where 3 is the number of levels of the inner factor B. (This type of 
recoding factors has other applications, esp. when plotting data, so 
recognizing how to do it quickly can be a big help.)

The PROC GLM commands to test for homogeneity of variance in this 
situation are the following:

PROC GLM;
CLASS class_var;
MODEL y= class_var;
MEANS class_var / HOVTEST= < option > ;
RUN;

HOVTEST= is an option entered on the MEANS statement. The classification 
factor must be listed on the CLASS, MODEL, and MEANS statements. If you 
have a two-way ANOVA (or with more factors), the HOVTEST option on the 
MEANS statement will not produce any test.

The available choices for homogeneity of variance tests are listed below. 
They request HOV tests for a one-way ANOVA where the groups are defined by 
the levels of the factor entered on the MEANS statement. Choices for 
HOVTEST= include:

 HOVTEST=bartlett
 HOVTEST=levene < ( TYPE= abs | square )>
 HOVTEST=bf
 HOVTEST=obrien < ( W=number )>

Specify the particular test you want for the <option>; if you do not 
specify a test, Levene's test with TYPE=square is the default.

The bartlett option specifies Bartlett's test which is an adaptation of 
the likelihood ratio procedure derived by Neyman and Pearson. It uses only 
sample variances and counts and is known to be *very* sensitive to 
non-normality. George Box pointed this out in 1953 when he made his famous 
statement:

  "To make a preliminary test on variances is rather like putting to sea
   in a row boat to find out whether conditions are sufficiently calm
   for an ocean liner to leave port!" - [Box, "Non-normality and tests
   on variances", 1953, Biometrika 40, pp. 318-335]

The levene option specifies Levene's test (either with absolute or squared 
deviations from the group means) which is generally the recommended method 
to examine whether the variances of data defined by two or more groups are 
equal. It can end up with an inflated Type I error rate. It essentially 
calculates the absolute value or squared value of each observation's 
residual from its group mean and then performs an analysis of variance on 
these positive deviations. Levene's test is less sensitive than other 
tests to non-normality (e.g., which is required for Bartlett's test).

If you choose Levene's test you may enter the TYPE= option in parentheses 
to specify whether to apply the test to the absolute residuals (TYPE=abs) 
or the squared residuals (TYPE=square), the default.

The bf option specifies Brown and Forsythe's variation of Levene's test. 
It takes the Levene approach, but works with the absolute differences from 
the group medians. This choice can be time-consuming for large datasets, 
but generally that is not an issue. The Brown-Forsythe method is perhaps 
the best choice if you have a one-way ANOVA problem.

The obrien option specifies O'Brien's test, which is basically a 
modification of HOVTEST=levene(TYPE=square). You can apply the W= option 
in parentheses to tune the variable to match the suspected kurtosis of the 
underlying distribution. By default, W=0.5.


2. Testing Means with the Unequal Variance Model

For the testing for equality of variances for an ANOVA, first consider the 
following small dataset which are measurements collected from 24 subjects 
randomly divided into one of three groups. Group 1 receives the standard 
treatment, group 2 receives a treatment similar to the standard, and group 
3 gets a new and different treatment. The objective is to determine if the 
new treatment from group 3 is better than the treatments received by 
groups 1 and 2.

DATA new; DROP y1 y2 y3;
LABEL group='Treatment';
INPUT group y1 y2 y3;
group=1; y=y1; OUTPUT;
group=2; y=y2; OUTPUT;
group=3; y=y3; OUTPUT;
CARDS;
1 3.791 3.122 4.761
2 5.174 2.514 3.884
3 3.019 3.850 6.840
4 3.218 2.564 9.150
5 3.079 4.486 2.776
6 4.054 3.199 5.398
7 3.131 3.406 6.405
8 2.822 3.986 2.115
;

PROC SORT DATA=new ;
BY group;
RUN;

PROC TABULATE DATA=new NOseps;
CLASS group;
VAR y;
TABLE group, y*(n*f=3.0 mean*f=7.2 var*f=7.3) / rts=12;
RUN;

The associated summary statistics collected across the 8 subjects in each 
treatment group:

--------------------------------
|          |         y         |
|          |-------------------|
|          | N | Mean  |  Var  |
|----------+---+-------+-------|
|Treatment |   |       |       |
|1         |  8|   3.54|  0.610|
|2         |  8|   3.39|  0.475|
|3         |  8|   5.17|  5.299|
--------------------------------

The table of summary statistics indicates that groups 1 and 2 have similar 
means and variances that are nearly equal (at least with a rather small 
sample of 8 subjects per group) whereas the mean for group 3 is somewhat 
larger and its variance appears to be much "larger". The issue is to 
determine if the assumption of equal variances no longer holds, and then 
to run a one-way ANOVA to determine if the means are equal.


3. Homogeneity of Variance Tests

For the example dataset introduced above, Levene's test for the 
homogeneity of variance tests is applied with PROC GLM with absolute 
deviations from the group means:

PROC GLM DATA=new;
CLASS group;
MODEL y = group / ss3;
MEANS group / HOVTEST=levene(TYPE=abs) welch;
OUTPUT OUT=rsd(keep=group y pred resid) P=pred R=resid;
RUN;

Levene's Test for Homogeneity of y Variance
ANOVA of Absolute Deviations from Group Means

                      Sum of        Mean
Source        DF     Squares      Square    F Value      Pr > F

group          2      7.8270      3.9135       5.86    0.009499
Error         21     14.0236      0.6678

The p-value from Lavene's test indicates sufficient evidence to indicate 
the variances across the three treatments are unequal.

Inference about two means in the presence of unequal variance is commonly 
known as the Behrens-Fisher problem. Welch's test is one choice for the 
situation when you have k populations with unequal means. When k=2 it is 
known as the Welch t-test.  In this example the welch option entered on 
the MEANS statement tests the hypothesis that the means of the 3 groups 
are equal when there is evidence the variances of the groups are unequal.

Welch's ANOVA for y

Source          DF    F Value    Pr > F

group       2.0000       2.08    0.1656
Error      12.7640

The p-value for group (0.1656) indicates means are not significantly 
different according to Welch's test.


4. Weighted Least Squares

Another approach with unequal variances is to apply a weighted least 
squares analysis. This test consists of first running an unweighted ANOVA 
on the data from the three groups, determine the variance of the residuals 
from each group, compute the inverse of these variances, and then merge 
this value back into the dataset to be applied as a weight. In the example 
above the residuals from the unweighted model were saved when Levene's 
test was computed. The remaining SAS steps are as follows:

* compute the variances within each group;

PROC MEANS DATA =rsd NOprint;
BY group; VAR resid;
OUTPUT OUT=vr(drop= _: ) VAR=vr N=n;
run;

* compute the weight for each group which is 1 over the group variance;

DATA vr; SET vr; wgt = 1/vr; RUN;

PROC PRINT DATA=vr NOobs; run;

group       vr      n      wgt

  1      0.60964    8    1.6403
  2      0.47518    8    2.1045
  3      5.29882    8    0.1887

These weights are now merged into the respective records from the main 
dataset.

DATA new;
MERGE new vr(keep=group wgt);
BY group;

PROC PRINT DATA=new; run;

obs    group      y        wgt

  1      1      3.791    1.6403
  2      1      5.174    1.6403
etc...
  9      2      3.122    2.1045
 10      2      2.514    2.1045
etc..
 17      3      4.761    0.1887
 18      3      3.884    0.1887etc..

* run a weighted ANOVA with the weights as the inverse of the group 
  variances;

PROC GLM DATA=new;
CLASS group;
WEIGHT wgt;
MODEL y = group / ss3;
RUN; quit;

The one-way weighted ANOVA from PROC GLM contains the following output:

Weight: wgt
                                Sum of
Source                DF       Squares   Mean Square  F Value  Pr > F

Model                  2    4.36649870    2.18324935     2.18  0.1376
Error                 21   21.00000000    1.00000000
Corrected Total       23   25.36649870

The pvalue for Model (0.1376) indicates the means are not different 
according to these weighted ANOVA calculations.


5. How can PROC MIXED handle this unequal variance situation?

PROC MIXED provides the capability to work directly with unequal 
variances, as it is able to apply different estimates of the variances for 
levels of a grouping factor specified as one of the CLASS variables. As a 
result, all observations having the same level of the variable entered on 
the GROUP= option of the REPEATED statement as shown below will have the 
same estimated variance.

ODS SELECT COVPARMS FITSTATISTICS TESTS3;

PROC MIXED DATA=new;
CLASS group;
MODEL y= group / solution ddfm=satterthwaite;
                * apply Sattherwaite ddfm option with unequal variances;
REPEATED / group=group; * specify the unequal variance model;
RUN;

Covariance Parameter Estimates

Cov Parm     Group      Estimate

Residual     group 1      0.6095
Residual     group 2      0.4750
Residual     group 3      5.2979

Fit Statistics

-2 Res Log Likelihood            68.8
AIC (smaller is better)          74.8

Type 3 Tests of Fixed Effects

              Num     Den
Effect         DF      DF    F Value      Pr > F

group           2    9.72       2.18    0.16465

This table contains a pvalue for group (p=0.16465) which is nearly the 
same as computed with the Welch test above (p=0.1656). (Without the 
ddfm=sattherthwaite option on the MODEL statement, the pvalue is the same 
as found with the weighted ANOVA with PROC GLM, p=0.1376). The results 
here were computed with far less effort. [Note: With unequal variances, 
one should add the ddfm=sattherthwaite option to the MODEL statement and 
the group= option to the REPEATED statement to the commands listed within 
the PROC MIXED step.] Another way to fit this unequal variance model with 
PROC MIXED is through the "Dispersion Effects Model" shown below.
You can also set equal variances for subsets of classification levels.  
That is, you can group the data for unequal variances by assigning values 
to a variable based on your knowledge of the variability of the data; you 
do not necessarily need to assume the variances across all levels of a 
classification variable are unequal. For the example data, the summary 
statistics indicate we may be able to assign a different variance to group 
3 while keeping the same pooled variance estimate for groups 1 and 2:

DATA new; SET new; grp=(group=3); RUN;

This dataset adds a dummy variable called grp (assigned to have the value 
1 for group 3, 0 otherwise) and will be entered into the statements for 
PROC MIXED to determine if one should compute a separate variance estimate 
for group 3 than compute a common variance for all three groups.

The following PROC MIXED step fits an unequal-variance ANOVA model for the 
grouping of variances as described above. Namely, you can test the 
equality of group variances by entering the option GROUP=grp on the 
REPEATED statement (the name grp must also appear on the CLASS statement).

REPEATED / GROUP=grp;

This choice defines the effect which specifies heterogeneity in the 
covariance structure of R. If all three variances are assumed unequal, 
GROUP=grp on the REPEATED statement; this will also be tested later.

To determine which model is better, compare the fit statistics of the 
various models such as AIC, BIC, etc. and select the model with a 
"substantially" smaller value.

The output from both the LSMEANS (with computed differences) and ESTIMATE 
statements are obtained for comparisons of means for groups 1, 2, and 3:

PROC MIXED DATA=new;
CLASS group grp;      * grp=0 for group=1 & 2;
                      * grp=1 for group=3;
MODEL y= group / solution ddfm=satterthwaite;
REPEATED / group=grp; * specify the unequal variance model;
LSMEANS group / diff ;
ESTIMATE 'group 1 vs 2' group 1 -1  0;
ESTIMATE 'group 1 vs 3' group 1  0 -1;
ESTIMATE 'group 2 vs 3' group 0  1 -1;
RUN;

Fit Statistics

-2 Res Log Likelihood            68.9
AIC (smaller is better)          72.9 **

The Fit Statistics are smaller than with this choice of the unequal 
variance model. Note for comparison note that the AIC value for the 
unweighted ANOVA model (results not shown here) is 83.7

Type 3 Tests of Fixed Effects

              Num     Den
Effect         DF      DF    F Value    Pr > F

group           2    9.75       2.16    0.1673


Least Squares Means
                          Standard
Treatment    Estimate       Error

   1           3.5360      0.2604
   2           3.3909      0.2604
   3           5.1661      0.8139 **

** The two standard errors for the lsmeans from group=1 and 2 are equal 
   and that the standard error is larger for group=3.

Here are the results from the Estimate statements which compare all three 
pairwise combinations of means.

Estimates
                            Standard
Label           Estimate       Error      DF    t Value    Pr > |t|

group 1 vs 2      0.1451      0.3682      14       0.39    0.6994
group 1 vs 3     -1.6301      0.8545    8.46      -1.91    0.0909
group 2 vs 3     -1.7752      0.8545    8.46      -2.08    0.0695

Differences of Least Squares Means (same results as from Estimate 
statements)

                                        Standard
Effect  Treatment  Treatment  Estimate     Error    DF  t Value  Pr > |t|

group   1          2            0.1451    0.3682    14     0.39  0.6994
group   1          3           -1.6301    0.8545  8.46    -1.91  0.0909
group   2          3           -1.7752    0.8545  8.46    -2.08  0.0695


The Standard Error of the difference for the comparison of group=1 and 
group=2 is much smaller than the comparisons for group=1,2 with group=3. 
The lsmeans for groups 1 and 2 appear to be approximately equal; however 
both lsmeans are less than the lsmean for group 3 with a much larger 
standard error. This is one result of group 1 and 2 assigned to have the 
same variance while the variance of group 3 is allowed to differ.

If the REPEATED statement is removed (that is, choose an ANOVA model which 
assumes equal variances), the estimated lsmeans are exactly the same yet 
with equal standard errors.

Least Squares Means

                    Standard
group    Estimate      Error

  1        3.5359     0.5157
  2        3.3908     0.5157
  3        5.1663     0.5157

This results in different p-values for the same comparisons because the 
standard errors of the differences are now all equal.

Differences of Least Squares Means

                              Standard
group    group    Estimate       Error    DF    t Value    Pr > |t|

 1        2         0.1451      0.7293    21       0.20    0.844
 1        3        -1.6304      0.7293    21      -2.24    0.036
 2        3        -1.7755      0.7293    21      -2.43    0.024

With the unequal variance model, the pvalues for the comparisons of groups 
1 and 2 with 3 are larger than the same pvalues with the equal variance 
model. In fact for those who take 0.05 as a typical cutoff for 
significance, the unequal variance model produces non-significant results 
whereas the equal variance model is more ambiguous and perhaps misleading.

Note how PROC GLM produces the same results as PROC MIXED with the equal 
variance assumption:

PROC GLM DATA=new;
CLASS group;     
MODEL y= group / ss3;
LSMEANS group / pdiff tdiff ;
ESTIMATE 'group 1 vs 2' group 1 -1  0;
ESTIMATE 'group 1 vs 3' group 1  0 -1;
ESTIMATE 'group 2 vs 3' group 0  1 -1;RUN; quit;

                         LSMEAN
group     y LSMEAN       Number

 1       3.53590603        1
 2       3.39078786        2
 3       5.16627977        3


Least Squares Means for Effect group
t for H0: LSMean(i)=LSMean(j) / Pr > |t|
          Dependent Variable: y

i/j              1             2             3

   1                    0.198985      -2.23555
                          0.8442        0.0364
   2      -0.19898                    -2.43454
            0.8442                      0.0239
   3      2.235554      2.434539


            0.0364        0.0239

The output from the 3 ESTIMATE statements is one way to easily display the 
standard error of the difference between two lsmeans:

                             Standard
Parameter         Estimate      Error    t Value    Pr > |t| 
            
group 1 vs 2       0.14512     0.7293       0.20    0.8449       
group 1 vs 3      -1.63037     0.7293      -2.24    0.0364       
group 2 vs 3      -1.77549     0.7293      -2.43    0.0239     


The Table of Fit Statistics produced for the unequal variance analysis is:

-2 Res Log Likelihood            68.9
AIC (smaller is better)          72.9
AICC (smaller is better)         73.6
BIC (smaller is better)          75.3

The Fit Statistics for the output produced with the equal group variances 
are all considerably larger, indicating the unequal variance model may be 
appropriate.

-2 Res Log Likelihood            81.7
AIC (smaller is better)          83.7
AICC (smaller is better)         83.9
BIC (smaller is better)          84.7

Since the Fit Statistics are smaller with the unequal variance assumption, 
it confirms the results from Levene's test presented above.

What if we had assumed all three variances are unequal? This model can be 
tested by adding GROUP=group to the REPEATED statement in the PROC MIXED 
step listed above. The Fit Statistics for this model is somewhat larger 
than the model with two residual variances, thus indicating the 2 variance 
model is preferred.

-2 Res Log Likelihood            68.8
AIC (smaller is better)          74.8
AICC (smaller is better)         76.2
BIC (smaller is better)          78.4


6. Multiple Comparisons in the Unequal Variance ANOVA Model

Tests for multiple comparisons of means such as Tukey (among several 
others) have not been considered in this section since they are only valid 
under the assumption of equal variances across groups. If this assumption 
is violated, that is the heteroscedastic model as described here is 
needed, one way to compute a multiple comparison of means test that 
adjusts for unequal variances is to apply the ADJUST=simulate option on 
the LSMEANS statement as described in Section 10.2 of "Multiple 
Comparisons and Multiple Tests" by Westfall, et. al.

ODS OUTPUT DIFFS=dfs(drop=effect alpha);
ODS EXCLUDE diffs;

PROC MIXED DATA=new;
CLASS group grp;      * grp=0 for group=1 & 2;
                      * grp=1 for group=3;
MODEL y= group / solution ddfm=satterthwaite;
REPEATED / group=grp; * specify the unequal variance model;
LSMEANS group / diff ADJUST=simulate(seed=92953) cl;
RUN;

proc print data=dfs NOobs; run;


group  _group  Estimate   StdErr    DF  tValue   Probt   Adjp   Adjustment

  1      2       0.1451   0.3682    14    0.39  0.6994   0.9133  Simulate
  1      3      -1.6301   0.8545  8.46   -1.91  0.0909   0.1742  Simulate
  2      3      -1.7752   0.8545  8.46   -2.08  0.0695   0.1342  Simulate

    Lower       Upper    AdjLower    AdjUpper

  -0.6447      0.9349     -0.8542      1.1444
  -3.5820      0.3218     -3.9490      0.6888
  -3.7271      0.1766     -4.0941      0.5436

Under this method of multiple comparisons, the adjusted p-values (column 
Adjp) become even less significant compared with the actual p-values 
(column Probt).


7. Dispersion Effects Model

Another Method to Model Unequal Variances

The same unequal variance component estimates can be found by modeling 
group as a "Dispersion Effect", or, a factor that influences both the mean 
and the variance of the response.  This is done with the local=exp(trt) 
option for the REPEATED statement as shown below.

ODS OUTPUT covparms=covparms;

PROC MIXED DATA=new NOitPrint;
CLASS group;
MODEL y =group / solution ddfm=satterthwaite;
REPEATED / local=exp(group); /*different variances estimated for the */
                             /*values of group variable */
run;

The resulting output (not shown) is the same as shown earlier with the 
GROUP=group option. For one, the type3 test of the fixed effects variable 
group is the same:

Type 3 Tests of Fixed Effects

              Num     Den
Effect         DF      DF    F Value     Pr > F

group           2    9.72       2.18    0.16465

However, the covariance parameters appear to be quite different. How can 
they be converted to variance estimates?

  CovParm      Estimate

EXP group 1     -0.6377
EXP group 2     -0.8870
Residual         1.1533

The first two rows contain the CovEsts values and the third is the 
Residual. For this Dispersion Effects model the actual variances are 
computed as:

Variance = Residual*Diag(U*CovEsts)

where in this case

U= { 1  0
     0  1
    -1 -1}

To form the individual group variance estimates, first transpose the 
estimates:

PROC TRANSPOSE DATA=covparms OUT=ParmTr(drop= _name_);
VAR estimate; ID covparm; run;

* CovParms Transposed;

PROC PRINT DATA=ParmTr NOobs; run;

    EXP_        EXP_
 group_1     group_2    Residual

 -0.6377     -0.8870      1.1533

* compute variance estimates in a DATA step;

DATA Vars;
set ParmTr;
Variance1=Residual*exp(EXP_group_1);
Variance2=Residual*exp(EXP_group_2);
Variance3=Residual*exp(-1*EXP_group_1-1*EXP_group_2);
run;

proc print data=Vars NOobs ; VAR var: ;
run;

Notice the variance estimates with local=exp(group) [Dispersion Effects 
Model] are the same values as computed with the GROUP=group option.

Variance1    Variance2    Variance3

   0.6095       0.4750       5.2979


Why would anyone want to get these estimates of the three different 
variances this harder way rather than just enter the GROUP= option? 
Sometimes, the factor that affects the response variance is continuous or 
"analytic". And we may want to estimate what the variance might be between 
the specific levels of our available data. These parameter estimates can 
be computed as a function of the design values, so one can interpolate 
between (in this case) our actual numeric coded 1,2,3 levels (Note: be 
sure to only interpolate between and not extrapolate beyond - which is 
always dangerous!)

More details on this unequal variance estimation process with PROC MIXED 
can be found in "SAS System for Mixed Models" (2006).

References

Littell, R.C., Milliken, G.A., Stroup, W.W., and Wolfinger, R.D. (2006). 
SAS System for Mixed Models, 2nd ed., SAS Institute, Cary, NC.

Milliken, George A. and Dallas E. Johnson, "Analysis of Messy Data: 
Analysis of Covariance" Chapman & Hall: Boca Raton, 2002.

Westfall, Peter H., Tobias, Randall D., Rom, Dror, Wolfinger, Russell D., 
and Hochberg, Yosef, "Multiple Comparisons and Multiple Tests Using the 
SAS System", Cary, NC: SAS Institute Inc., 1999, 416 pp.