Distribution System Simulation Action Set

This section contains PROC CAS code.

Note: Input data must be accessible in your CAS session, either as a CAS table or as a transient-scope table. A CAS table has a two-level name: the first level is your CAS engine libref, and the second level is the table name. You refer to this table in the CAS procedure by specifying only the second level. For more information about two-level names, see Chapter 2, Shared Concepts. A transient-scope table is called directly from the action and exists in memory for the duration of the action. For more information about accessing data, see SAS Viya: System Programming Guide. For more information about PROC CAS and programming in CASL, see SAS Cloud Analytic Services: CASL Programmer’s Guide and SAS Cloud Analytic Services: CASL Reference.

You are an industrial engineer who is interested in how different equipment affects a certain characteristic of a manufactured product. You know that changing the equipment might shift the expected mean value of the characteristic, so analysis of variance (ANOVA) is an appropriate technique. However, you also know that the characteristic of interest is not normally distributed. Calculations from past data indicate that the characteristic’s values are skewed and heavy-tailed, with a skewness coefficient of about 1.5 and a kurtosis of 7. You are concerned that ANOVA F tests might not be valid for these kinds of data. You can use the simsk action to simulate data similar to what you might observe in practice, and then you can use the simulated data to study the behavior of ANOVA tests.

The following statements use the simsk action to simulate 100 values from a distribution that has mean 0, standard deviation 1, skewness 1.5, and kurtosis 7:

proc cas;
   simSystem.simsk /
      seed=12345
      nsimulations=100
      skewness=1.5
      kurtosis=7
      outtable={ name="Sim" };
run;

Because the system parameter is not specified, the default Pearson system of distributions is used. The seed parameter specifies a value for the random number generator in order to make the results reproducible. Because no values for the mean and standard deviation were specified, the chosen distribution has the default values of 0 and 1, respectively.

The simsk action determines the distribution of the Pearson system that has the specified skewness and kurtosis, and it displays the family and the parameters of this distribution. It also displays the theoretical and observed moments of this distribution for the sampled values, as shown in Figure 3.

As shown in the table of Pearson distribution parameters, the selected distribution is in the Type VI family of the Pearson system, which also happens to be equivalent to a shifted and scaled F distribution with 6.777 and 38 degrees of freedom. The sample skewness and especially the sample kurtosis in the moments table are not all that close to the specified values, but this is because of the relatively small sample size.

To address the original question of how well standard ANOVA works with this distribution, you can first simply run the regression.glm action on these data, with five randomly generated groups, as follows:

proc cas;
   ds2.runds2 program="
      data sim / overwrite=yes;
         declare double a y;
         method run();
            set sim;
            a = mod( iObs - 1, 5 ) + 1;
            y = Variate;
         end;
      enddata;";
   run;
   regression.glm /
      table={ name="Sim" }
      class={ "a" }
      model={ depVar="y" effects={ "a" } }
      display={ names="Anova" };
   run;
quit;

Because all the data have the same mean, you expect to get a nonsignificant F test for the group effect. Figure 4 shows that this is precisely what you get.

However, one F test does not really say much about the general characteristics of ANOVA on such data. A more extensive way to explore the question is to generate many replicates of such skewed and kurtotic data and to see how often the F test is rejected. For a valid statistical technique, the F test on these 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 the simsk action with the nReplicates parameter to generate 1,000 samples of size 100 from the distribution. The parameters of the distribution are the same as in Figure 3, but the sample moments are now averaged over the 1,000 samples, as shown in Figure 5.

proc cas;
   simSystem.simsk /
      seed=12345
      skewness=1.5
      kurtosis=7
      nsimulations=100
      nreplicates=1000
      outtable={name="Sim" replace=True};
run;

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 population moments, with a bias that can be rather high even for moderately large samples if the kurtosis is large.

In the following code, the regression.glm action step shown previously is performed on each of these 1,000 samples of size 100, and the runds2 and simple.summary actions are used to compute how often the F test is rejected at the 10%, 5%, and 1% levels:

proc cas;
   ds2.runds2 program="
      data sim / overwrite=yes;
         declare double a y;
         method run();
            set sim;
            a = mod( iObs - 1, 5 ) + 1;
            y = Variate;
         end;
      enddata;";
   run;
   regression.glm /
      table={ name="Sim" groupby={name="Rep"} }
      class={ "a" }
      model={ depVar="y" effects={ "a" } }
      display={ excludeall=True }
      outputTables={names={Anova="AOV"}}
      ;
   run;
   ds2.runds2 program="
      data AOV / overwrite=yes;
         declare double Sig10 Sig05 Sig01;
         method run();
            set AOV;
            if ( Source = 'Model' );
            Sig10 = ( ProbF < 0.10 );
            Sig05 = ( ProbF < 0.05 );
            Sig01 = ( ProbF < 0.01 );
         end;
      enddata;";
   run;
   simple.summary /
      table={name="AOV"},
      inputs={"Sig10", "Sig05", "Sig01"},
      casOut={name="AOV_summary", replace=True}
      ;
   run;
quit;
proc print data=mycas.AOV_summary noobs;
   var _Column_ _Mean_;
run;

The results, shown in Figure 6, indicate that the simulated significances for ANOVA on these data are reasonably close to the nominal levels. The implication is that you can be confident in using ANOVA with these nonnormal data.