SAS/IML Robust Regression Examples
Example 9.4: Brainlog Data
The following data, consisting of the body weights (in kilograms) and brain weights (in grams) of N=28 animals, are reported by Jerison (1973) and can be found also in Rousseeuw and Leroy
(1987, p. 57). Instead of the original data, this example uses the logarithms of the measurements of the two variables.
title "*** Brainlog Data: Do MVE ***";
aa={ 1.303338E-001 9.084851E-001 ,
2.6674530 2.6263400 ,
1.5602650 2.0773680 ,
1.4418520 2.0606980 ,
1.703332E-002 7.403627E-001 ,
4.0681860 1.6989700 ,
3.4060290 3.6630410 ,
2.2720740 2.6222140 ,
2.7168380 2.8162410 ,
1.0000000 2.0606980 ,
5.185139E-001 1.4082400 ,
2.7234560 2.8325090 ,
2.3159700 2.6085260 ,
1.7923920 3.1205740 ,
3.8230830 3.7567880 ,
3.9731280 1.8450980 ,
8.325089E-001 2.2528530 ,
1.5440680 1.7481880 ,
-9.208187E-001 .0000000 ,
-1.6382720 -3.979400E-001 ,
3.979400E-001 1.0827850 ,
1.7442930 2.2430380 ,
2.0000000 2.1959000 ,
1.7173380 2.6434530 ,
4.9395190 2.1889280 ,
-5.528420E-001 2.787536E-001 ,
-9.136401E-001 4.771213E-001 ,
2.2833010 2.2552720 };
By default, the MVE subroutine (like the MINVOL subroutine) uses only 1500 randomly selected subsets rather than all subsets.
The following specification of the options vector requires that all 3276 subsets of 3 cases out of 28 cases are generated and
evaluated:
title2 "***MVE for BrainLog Data***";
title3 "*** Use All Subsets***";
optn = j(8,1,.);
optn[1]= 3; /* ipri */
optn[2]= 1; /* pcov: print COV */
optn[3]= 1; /* pcor: print CORR */
optn[6]= -1; /* nrep: all subsets */
call mve(sc,xmve,dist,optn,aa);
Specifying optn[1]=3, optn[2]=1, and optn[3]=1 requests that all output be printed. Therefore, the
first part of the output shows the classical scatter and correlation matrix.
Output 9.4.1: Some Simple Statistics
|
Minimum Volume Ellipsoid (MVE) Estimation
|
|
Consider Ellipsoids Containing 15 Cases.
|
|
Classical Covariance Matrix
|
|
|
VAR1
|
VAR2
|
|
VAR1
|
2.6816512357
|
1.3300846932
|
|
VAR2
|
1.3300846932
|
1.0857537552
|
|
Classical Correlation Matrix
|
|
|
VAR1
|
VAR2
|
|
VAR1
|
1
|
0.7794934643
|
|
VAR2
|
0.7794934643
|
1
|
|
Classical Mean
|
|
VAR1
|
1.6378572186
|
|
VAR2
|
1.9219465964
|
|
The second part of the output shows the results of the combinatoric optimization (complete subset sampling).
Output 9.4.2: Iteration History for MVE
|
There are 3276 subsets of 3 cases out of 28 cases.
|
|
The algorithm will draw 1500 random subsets of 3 cases.
|
|
Random Subsampling for MVE
|
|
Subset
|
Singular
|
Best
Criterion
|
Percent
|
|
375
|
0
|
0.453147
|
25
|
|
750
|
0
|
0.453147
|
50
|
|
1125
|
0
|
0.453147
|
75
|
|
1500
|
0
|
0.453147
|
100
|
|
Minimum Criterion= 0.4531471437
|
|
Among 1500 subsets 0 are singular.
|
|
Observations of Best Subset
|
|
1
|
22
|
2
|
Initial MVE Location
Estimates
|
|
VAR1
|
1.5140266
|
|
VAR2
|
1.9259543667
|
|
Initial MVE Scatter Matrix
|
|
|
VAR1
|
VAR2
|
|
VAR1
|
7.5699818309
|
5.2533273517
|
|
VAR2
|
5.2533273517
|
3.7329226108
|
|
The third part of the output shows the optimization results after local improvement.
Output 9.4.3: Table of MVE Results
|
Final MVE Estimates (Using Local Improvement)
|
|
Number of Points with Nonzero Weight=23
|
|
Robust MVE Location Estimates
|
|
VAR1
|
1.3154029661
|
|
VAR2
|
1.8568731174
|
|
Robust MVE Scatter Matrix
|
|
|
VAR1
|
VAR2
|
|
VAR1
|
2.139986054
|
1.6068556533
|
|
VAR2
|
1.6068556533
|
1.2520384784
|
Eigenvalues of Robust
Scatter Matrix
|
|
VAR1
|
3.363074897
|
|
VAR2
|
0.0289496354
|
|
Robust Correlation Matrix
|
|
|
VAR1
|
VAR2
|
|
VAR1
|
1
|
0.9816633012
|
|
VAR2
|
0.9816633012
|
1
|
|
The final output presents a table containing the classical Mahalanobis distances, the robust distances, and the weights
identifying the outlier observations.
Output 9.4.4: Mahalanobis and Robust Distances
|
Classical Distances and Robust (Rousseeuw) Distances
|
|
Unsquared Mahalanobis Distance and
|
|
Unsquared Rousseeuw Distance of Each Observation
|
|
N
|
Mahalanobis Distances
|
Robust Distances
|
Weight
|
|
1
|
1.006591
|
0.855347
|
1.000000
|
|
2
|
0.695261
|
1.477050
|
1.000000
|
|
3
|
0.300831
|
0.239828
|
1.000000
|
|
4
|
0.380817
|
0.517719
|
1.000000
|
|
5
|
1.146485
|
1.108362
|
1.000000
|
|
6
|
2.644176
|
10.599341
|
0
|
|
7
|
1.708334
|
1.808455
|
1.000000
|
|
8
|
0.706522
|
0.690099
|
1.000000
|
|
9
|
0.858404
|
1.052423
|
1.000000
|
|
10
|
0.798698
|
2.077131
|
1.000000
|
|
11
|
0.686485
|
0.888545
|
1.000000
|
|
12
|
0.874349
|
1.035824
|
1.000000
|
|
13
|
0.677791
|
0.683978
|
1.000000
|
|
14
|
1.721526
|
4.257963
|
0
|
|
15
|
1.761947
|
1.716065
|
1.000000
|
|
16
|
2.369473
|
9.584992
|
0
|
|
17
|
1.222253
|
3.571700
|
0
|
|
18
|
0.203178
|
1.323783
|
1.000000
|
|
19
|
1.855201
|
1.741064
|
1.000000
|
|
20
|
2.266268
|
2.026528
|
1.000000
|
|
21
|
0.831416
|
0.743545
|
1.000000
|
|
22
|
0.416158
|
0.419923
|
1.000000
|
|
23
|
0.264182
|
0.944610
|
1.000000
|
|
24
|
1.046120
|
2.289334
|
1.000000
|
|
25
|
2.911101
|
11.471953
|
0
|
|
26
|
1.586458
|
1.518721
|
1.000000
|
|
27
|
1.582124
|
2.054593
|
1.000000
|
|
28
|
0.394664
|
1.675651
|
1.000000
|
|
Distribution of Robust Distances
|
|
MinRes
|
1st Qu.
|
Median
|
Mean
|
3rd Qu.
|
MaxRes
|
|
0.2398280857
|
0.8719458303
|
1.4978856181
|
2.4419474031
|
2.0658619775
|
11.471952618
|
|
Cutoff Value = 2.7162030315
|
The cutoff value is the square root of the 0.975 quantile of the chi square distribution
with 2 degrees of freedom.
|
There are 5 points with large robust distances receiving zero weights. These may include
boundary cases. Only points whose robust distances are subs tantially larger than the
cutoff value should be considered outliers.
|
|
Again, you can call the subroutine scatmve(), which is included in the sample library in file robustmc.sas, to
plot the classical and robust confidence ellipsoids:
optn = j(8,1,.); optn[6]= -1;
vnam = { "Log Body Wgt","Log Brain Wgt" };
filn = "brl";
titl = "BrainLog Data: Use All Subsets";
call scatmve(2,optn,.9,aa,vnam,titl,1,filn);
The output follows.
Output 9.4.5: BrainLog Data: Classical and Robust Ellipsoid
Statistics and Operations Research | SAS/IML Software
