The MBC Procedure

This example uses the Old Faithful data on geyser eruptions (Cook and Weisberg 1982) to illustrate a mixture of multivariate Gaussian components. This data set contains time before eruption and duration of eruption for Old Faithful geyser in Yellowstone National Park in the United States. The following DATA step creates the data table mycas.getStarted, which consists of 272 observations on the duration of each eruption and the wait time until the next eruption. Both duration and wait time are measured in minutes.

data mycas.getStarted;
  input Duration Wait @@;
  datalines;
3.600 79  1.800 54  3.333 74  2.283 62  4.533 85  2.883 55  4.700 88  3.600 85
1.950 51  4.350 85  1.833 54  3.917 84  4.200 78  1.750 47  4.700 83  2.167 52

   ... more lines ...   

1.983 43  2.250 60  4.750 75  4.117 81  2.150 46  4.417 90  1.817 46  4.467 74
;

These statements assume that your CAS engine libref is named mycas, but you can substitute any appropriately defined CAS engine libref.

The following statements produce a scatter plot of the duration and wait times:

proc sgplot data=mycas.getStarted;
  scatter x=duration y=wait / markerattrs=(symbol=circlefilled size=4);
run; quit;

It is clear from the scatter plot in Figure 12.1 that there are two clusters of eruptions. The first consists of short, frequent eruptions, and the second consists of long, infrequent eruptions. However, the shape of these clusters is not spherical. The MBC procedure can help you explore the shapes.

Figure 12.1: Scatter Plot of Old Faithful Data

The following statements fit several different mixture models to these data:

proc mbc data=mycas.getStarted nclusters=(1 to 5)
                               noise=(YES NO)
                               covstruct=(EEV EII VVV)
                               seed=1418410433;
  var duration wait;
  output out=mycas.scores copyvars=(duration wait);
run;

The NCLUSTERS= option list fits separate models by using one to five Gaussian clusters. The NOISE= option list fits these models with and without noise. The COVSTRUCT= option list fits each of these models by using three different covariance structures. Finally, the SEED= option specifies a starting random seed to generate initial cluster weights. In total, PROC MBC fits 30 different models, corresponding to the unique combinations of each element in the NCLUSTERS=, NOISE=, and COVSTRUCT= option lists.

The VAR statement identifies the two variables, Duration and Wait, that define the points to be clustered.

The OUTPUT statement specifies a CAS data table to contain the posterior cluster weights for each observation. The COPYVARS= option includes the two analysis variables, Duration and Wait, in the output data table.

After fitting all 30 models, PROC MBC chooses the model that has the best Bayesian information criterion (BIC) value. BIC is a quantity composed of the log likelihood and a penalty term. Smaller BIC values indicate models with better fit. The "Fit Summary" table in Figure 12.2 shows model information for the 10 models that have the best (smallest) BIC values, sorted so that the best model appears first. The selected model has two nonspherical clusters, as suggested by the plot in Figure 12.1.

Figure 12.3 shows the "Model Information" table. This table summarizes the characteristics of the selected model and of the modeling technique that was used.

The "Parameter Estimates" table in Figure 12.4 contains the estimates of the centers of the two clusters and the estimates of the covariance for the multivariate Gaussian that is represented by each cluster.

If you compare the mean estimates to the plot in Figure 12.1, you can see that the estimates correspond well to the data.

Figure 12.1 suggests that one cluster is slightly heavier than the other. The "Mixing Estimates" table in Figure 12.5 supports this impression.

The following statements use the output table mycas.scores to produce a plot of the maximum posterior weights for each observation:

data scores;
  set mycas.scores;
  maxwt = max(of next:);
run;

proc sgplot data=scores noautolegend;
  * ---- shaded corresponding to maximum posterior weight --------- ;
  scatter x=duration y=wait   / markerattrs=(symbol=circlefilled)
                                colorresponse=maxwt
                                colormodel=TwoColorRamp;
  scatter x=duration y=wait   / markerattrs=(symbol=circle color=black);
run;
quit;

Figure 12.6 plots the observations with shading to indicate the strength of the affinity that each point has for the cluster that it is more strongly associated with. Only one point has light shading, indicating that it is in a region where the larger posterior cluster weight is not as large as it is in other regions. This indicates a region where the posterior weights for the two clusters are closer in value.

Figure 12.6: Old Faithful Data and Selected Model