The SIMSYSTEM Procedure

This example demonstrates a basic application of the SIMSYSTEM procedure in which the MOMENTS statement is used to simulate samples from a small number of nonnormal distributions with specified skewness (denoted by ) and kurtosis (denoted by ); see Moments and Moment Ratios. The example concludes by explaining the use of the MOMENTGRID statement to simulate samples from a large number of nonnormal distributions whose skewness and kurtosis, also referred to as moment ratios, vary evenly over a grid.

An industrial engineer in a manufacturing facility is investigating the effect of five operators with different levels of training on a critical functional parameter of a product. The engineer obtains 100 measurements of the parameter, 20 for each operator. A one-way analysis of variance (ANOVA) seems like an appropriate method for testing the hypothesis of no training effect, because if there is an effect the level of training would shift the expected (mean) value of the parameter.

However, it cannot be assumed that the measurements are normally distributed. The distribution of measurements taken previously was found to be skewed and heavy-tailed, with a skewness of 1.5 and a kurtosis of 7. There is concern that ANOVA F tests might not be valid for this type of data. You can use PROC SIMSYSTEM to simulate data with this particular combination of skewness and kurtosis, and with combinations that depart from the observed combination (including a skewness of 0 and a kurtosis of 3, which correspond to the normal distribution). You can then use the simulated data to study the behavior of ANOVA tests.

To start, consider how to simulate data from a single distribution with the observed combination of skewness and kurtosis. The following statements use PROC SIMSYSTEM to simulate 100 values from the Pearson distribution that has mean 0, standard deviation 1, skewness 1.5, and kurtosis 7:

ods graphics on;
proc simsystem system=pearson n=100 seed=12345;
   moments skewness=1.5 kurtosis=7;
   output out=mycas.sim;
run;

The SYSTEM= option specifies the system of distributions to sample from. The N= option specifies the number of simulated values, and the SEED= option specifies the seed for the random number generator in order to make the results reproducible. The MOMENTS statement specifies the values of skewness and kurtosis. Note that skewness and kurtosis are defined somewhat differently for theoretical distribution in the statistics literature and in SAS; see Definitions and Interpretations of Skewness and Kurtosis for the definitions used by the SIMSYSTEM procedure. Because no values are specified for the mean and standard deviation, default values of 0 and 1 are used, respectively. The OUTPUT statement requests a CAS table named SIM that contains the simulated values.

The SIMSYSTEM procedure determines the unique distribution in the Pearson system that has the specified skewness, kurtosis, mean, and standard deviation. Distributions in the Pearson system are classified into families, also referred to as Pearson types, which are designated by Roman numerals; see Pearson System of Distributions. The procedure determines both the Pearson type and the mathematical parameters for the distribution, as displayed by the Parameters for Pearson Distributions table in Figure 1. The specified moments and the sample moments for the simulated values are displayed in the Distribution Moments table.

As shown in Figure 1, the Pearson distribution determined by the procedure belongs to the Type VI family, which corresponds to a shifted and scaled F distribution with 6.777 and 38 degrees of freedom; see Pearson Type VI Family. Because the sample size is relatively small, the sample skewness and especially the sample kurtosis in the Distribution Moments table are not all that close to the specified values.

In addition to these tables, if ODS Graphics is enabled, PROC SIMSYSTEM displays a graphical map of the Pearson system of distributions in which the specified distribution is marked. The procedure also displays the theoretical density for the distribution and a histogram for the 100 sampled values. These plots are shown in Figure 2 and Figure 3.

Figure 2: Pearson Distributions

Figure 3: Density of Selected Distribution and Histogram of Simulated Values

To carry out an F test with the simulated data, you can create a categorical variable (named A) that randomly groups the measurements by operator and run PROC REGSELECT, as follows:

data mycas.sim; set mycas.sim;
   a = mod(iObs-1,5)+1;
   y = variate;
run;

proc regselect data=mycas.sim;
   class a;
   model y = a;
run;

Because all the measurements come from a distribution with the same mean, you expect that the F test will fail to reject the hypothesis of no operator effect. Figure 4 shows that this is what you get.

However, one F test does not say anything about the general characteristics of ANOVA on such data. One way to explore this is to generate many replicates of such null data and see how often the F test is rejected. In addition, because you cannot be certain that the distribution of the functional parameter has a precise skewness of 1.5 and a kurtosis of 7, you need to study how well ANOVA works for a variety of skewness-kurtosis combinations.

For a valid analysis of variance, the F test on null data should be rejected about 5% of the time at the 0.05 level, and likewise 1% and 10% of the time at the 0.01 and 0.10 levels, respectively. The following code checks this by using PROC SIMSYSTEM with the NREP= option to generate 1,000 samples of size 100, and by using multiple skewness and kurtosis pairs in the MOMENTS statement to draw these samples from a selection of six Pearson distributions whose skewness ranges from 0 to 1.5 and whose kurtosis ranges from 3 to 7.

proc simsystem system=pearson n=100 seed=12345 nrep=1000;
   moments skewness = 0 0 0.75 0.75 0.75 1.5
           kurtosis = 3 5 3    5    7    7  ;
   output out=mycas.sim;
run;

The Parameters table in Figure 5 displays the Pearson family types and parameters determined by the procedure for these six distributions. The first distribution is the standard normal distribution, which is the Pearson distribution with a skewness of 0 and a kurtosis of 3. The second distribution is a Pearson Type VII distribution, also known as a t distribution, with 7 degrees of freedom; see Pearson Type VII Family. The third distribution is a Pearson Type I distribution, also known as a beta distribution, with shape parameters 1.342 and 3.769; see Pearson Type I Family. The fourth and fifth distributions are Pearson Type IV distributions, which do not correspond to any well-known family of distributions; see Pearson Type IV Family. The sixth distribution is the Pearson Type VI distribution, also known as an F distribution, that was determined at the beginning of this example; see Pearson Type VI Family.

For each of these distributions, the Distribution Moments table in Figure 5 displays the specified moments and the the sample moments averaged over the 1000 replicate samples.

Figure 6 locates the specified skewness-kurtosis combinations for the six distributions in a moment-ratio map of the Pearson system. Figure 7 shows the shapes of the six corresponding density functions.

Figure 6: Selected Pearson Distributions

Figure 7: Densities and Histograms of Selected Pearson Distributions

Notice how the Pearson system provides all of the distribution shapes necessary to match the entire set of skewness-kurtosis combinations specified in the MOMENTS statement. This could not be accomplished with a single family of distributions.

The average sample moments in Figure 5 are all close to the specified values, but there are noticeable discrepancies for the standard deviation, for the skewness, and especially for the kurtosis. This is because these sample moments are biased estimates of the specified population moments, with a bias that can be rather high even for moderately large samples if the kurtosis is large.

In the following statements, the PROC REGSELECT step shown earlier is performed on each of the 1,000 samples of size 100, for each combination of skewness and kurtosis values, and the DATA step and PROC MDSUMMARY are used to compute how often the F test is rejected at the 10%, 5%, and 1% levels. The ODS SELECT NONE statement suppresses printed output for all the ANOVA analyses.

data mycas.sim; set mycas.sim;
   a = mod(iObs-1,5)+1;
   y = variate;
run;

ods select none;
proc regselect data=mycas.sim;
   by Skewness Kurtosis Rep;
   class a;
   model y = a;
   displayout ANOVA=AOV;
run;
ods select all;

data mycas.AOV; set mycas.AOV;
   where (Source = "Model");
   Rej10 = (ProbF < 0.10);
   Rej05 = (ProbF < 0.05);
   Rej01 = (ProbF < 0.01);
proc mdsummary data=mycas.AOV;
   groupby Skewness Kurtosis;
   var Rej:;
   output out=mycas.AOV_summary;
proc transpose data=mycas.AOV_summary(rename=(_Column_=_Name_)) out=summary;
   by notsorted Skewness Kurtosis;
   var _Mean_;
proc sort data=summary;
   by Skewness Kurtosis;
proc print data=summary noobs;
   var Skewness Kurtosis Rej01 Rej05 Rej10;
   format Rej01 Rej05 Rej10 5.3;
run;

The results, shown in Figure 8, indicate that the simulated rejection rates are reasonably close to the nominal levels over the entire set of skewness-kurtosis combinations. The conclusion is that you can be confident in using ANOVA with these nonnormal data.

This example demonstrates how easy it is to simulate data from a variety of distributions with a single invocation of the SIMSYSTEM procedure. In practice, you will often want to simulate data from a larger number of distributions whose moment ratios lie on an evenly spaced grid of skewness-kurtosis combinations. The MOMENTGRID statement is useful for this purpose, as shown in the following statements.

proc simsystem system=pearson n=100 seed=12345 nrep=1000 plot(only)=mrmap;
   momentgrid skewness = 0 to 1.5 by 0.25
              kurtosis = 3 to 7   by 1;
   output out=mycas.sim;
run;

When you use the MOMENTGRID statement, each value of that you specify in the SKEWNESS= option is paired with all the values of that you specify in the KURTOSIS= option to form a grid of combinations of and .

The Pearson distributions determined by these combinations are mapped in Figure 8.

Figure 9: Selected Pearson Distributions

Although the preceding MOMENTGRID statement requests Pearson distributions for a grid of moment ratio combinations, the combination with skewness 1.5 and kurtosis 3 is not feasible because it does not satisfy the condition .

In general, you can use the MOMENTGRID statement to simulate data from distributions specified by a comprehensive grid of evenly spaced combinations of skewness and kurtosis. This enables you to explore the effects of a wide variety of distributional shapes and to form conclusions that are framed in terms of skewness and kurtosis.