Partial Least Squares 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.

The example in this section illustrates basic features of the pls action in the pls action set. The data are reported in Umetrics (1995); the original source is Lindberg, Persson, and Wold (1983). Suppose you are researching pollution in the Baltic Sea and you want to use the fluorescence spectra of seawater samples to determine the amounts of three compounds present in those samples: lignin sulfonate (ls: pulp industry pollution), humic acids (ha: natural forest products), and optical whitener from laundry detergent (dt). Spectrometric calibration is a type of problem in which partial least squares can be very effective. The predictors are the spectra emission intensities at different frequencies in a sample spectrum, and the responses are the amounts of various chemicals in the sample.

For the purpose of calibrating the model, samples that have known compositions are used. The calibration data consist of 16 samples of known concentrations of ls, ha, and dt, with spectra based on 27 frequencies (or, equivalently, wavelengths). In order to demonstrate the use of test set validation, the data contain the variable Role, which is used to assign observations to the training and testing roles. In this case, the training role has nine samples and the testing role has seven samples.

The following DATA step creates the mycas.Sample data table, which provides the calibration data, in your CAS session. This DATA step assume that your CAS engine libref is named mycas, but you can substitute any appropriately defined CAS engine libref.

data mycas.Sample;
   input obsnam $ v1-v27 ls ha dt Role $5. @@;
   datalines;
EM1   2766 2610 3306 3630 3600 3438 3213 3051 2907 2844 2796
      2787 2760 2754 2670 2520 2310 2100 1917 1755 1602 1467
      1353 1260 1167 1101 1017         3.0110  0.0000   0.00 TRAIN
EM2   1492 1419 1369 1158  958  887  905  929  920  887  800
       710  617  535  451  368  296  241  190  157  128  106
        89   70   65   56   50         0.0000  0.4005   0.00 TEST
EM3   2450 2379 2400 2055 1689 1355 1109  908  750  673  644
       640  630  618  571  512  440  368  305  247  196  156
       120   98   80   61   50         0.0000  0.0000  90.63 TRAIN
EM4   2751 2883 3492 3570 3282 2937 2634 2370 2187 2070 2007
      1974 1950 1890 1824 1680 1527 1350 1206 1080  984  888
       810  732  669  630  582         1.4820  0.1580  40.00 TEST
EM5   2652 2691 3225 3285 3033 2784 2520 2340 2235 2148 2094
      2049 2007 1917 1800 1650 1464 1299 1140 1020  909  810
       726  657  594  549  507         1.1160  0.4104  30.45 TEST
EM6   3993 4722 6147 6720 6531 5970 5382 4842 4470 4200 4077
      4008 3948 3864 3663 3390 3090 2787 2481 2241 2028 1830
      1680 1533 1440 1314 1227         3.3970  0.3032  50.82 TRAIN
EM7   4032 4350 5430 5763 5490 4974 4452 3990 3690 3474 3357
      3300 3213 3147 3000 2772 2490 2220 1980 1779 1599 1440
      1320 1200 1119 1032  957         2.4280  0.2981  70.59 TRAIN
EM8   4530 5190 6910 7580 7510 6930 6150 5490 4990 4670 4490
      4370 4300 4210 4000 3770 3420 3060 2760 2490 2230 2060
      1860 1700 1590 1490 1380         4.0240  0.1153  89.39 TRAIN
EM9   4077 4410 5460 5857 5607 5097 4605 4170 3864 3708 3588
      3537 3480 3330 3192 2910 2610 2325 2064 1830 1638 1476
      1350 1236 1122 1044  963         2.2750  0.5040  81.75 TEST
EM10  3450 3432 3969 4020 3678 3237 2814 2487 2205 2061 2001
      1965 1947 1890 1776 1635 1452 1278 1128  981  867  753
       663  600  552  507  468         0.9588  0.1450 101.10 TRAIN
EM11  4989 5301 6807 7425 7155 6525 5784 5166 4695 4380 4197
      4131 4077 3972 3777 3531 3168 2835 2517 2244 2004 1809
      1620 1470 1359 1266 1167         3.1900  0.2530 120.00 TRAIN
EM12  5340 5790 7590 8390 8310 7670 6890 6190 5700 5380 5200
      5110 5040 4900 4700 4390 3970 3540 3170 2810 2490 2240
      2060 1870 1700 1590 1470         4.1320  0.5691 117.70 TEST
EM13  3162 3477 4365 4650 4470 4107 3717 3432 3228 3093 3009
      2964 2916 2838 2694 2490 2253 2013 1788 1599 1431 1305
      1194 1077  990  927  855         2.1600  0.4360  27.59 TRAIN
EM14  4380 4695 6018 6510 6342 5760 5151 4596 4200 3948 3807
      3720 3672 3567 3438 3171 2880 2571 2280 2046 1857 1680
      1548 1413 1314 1200 1119         3.0940  0.2471  61.71 TRAIN
EM15  4587 4200 5040 5289 4965 4449 3939 3507 3174 2970 2850
      2814 2748 2670 2529 2328 2088 1851 1641 1431 1284 1134
      1020  918  840  756  714         1.6040  0.2856 108.80 TEST
EM16  4017 4725 6090 6570 6354 5895 5346 4911 4611 4422 4314
      4287 4224 4110 3915 3600 3240 2913 2598 2325 2088 1917
      1734 1587 1452 1356 1257         3.1620  0.7012  60.00 TEST
;

To isolate a few underlying spectral factors that provide a good predictive model, you can fit a PLS model to the 16 samples by using the pls action in the pls action set in the following PROC CAS statements:

proc cas;
   action pls.pls /
      table='SAMPLE',
      method='PLS',
      model={depVars={'ls', 'ha', 'dt'},
             effects={'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                      'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                      'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                      'v23', 'v24', 'v25', 'v26', 'v27'}};
run;

By default, the pls action extracts at most 15 factors. The default output from this analysis is presented in Output 18.1.1 and Output 18.1.2.

Output 18.1.1 displays the "Model Information," "Dimensions," and "Number of Observations" tables.

The "Model Information" table identifies the data source and shows that the factor extraction method is partial least squares regression and that the nonlinear iterative partial least squares (NIPALS) algorithm (which is the default) is used to compute extracted PLS factors.

The "Dimensions" table shows the number of response variables, the number of effects, the number of predictor parameters, and the number of factors to extract.

The "Number of Observations" table shows that all 16 of the sample observations in the input data are used in the analysis because all the samples contain complete data.

Output 18.1.2 lists the amount of variation, both individual and cumulative, that is accounted for by each of the 15 factors. All the variation in both the predictors and the responses is accounted for by only 15 factors because there are only 16 sample observations. More important, almost all the variation is accounted for by even fewer factors—one or two for the predictors and three to eight for the responses.

A PLS model is not complete until you choose the number of factors. You can choose the number of factors by using test set validation, in which the data table is divided into two groups called the training data and the test data. You fit the model to the training data, and then you check the capability of the model to predict responses for the test data. The predicted residual sum of squares (PRESS) statistic is based on the residuals that are generated by this process.

To select the number of extracted factors by test set validation, you use the partitionByVar or partitionByFrac parameter to specify how to logically divide observations in the input data table into two subsets for model training and testing. For example, you can designate a variable in the input data table and a set of formatted values of that variable to determine the role of each observation, as in the following PROC CAS statements. The ODS TRACE command displays the full paths of the output tables in the SAS log.

ods trace on;
proc cas;
   action pls.pls /
      table='SAMPLE',
      method='PLS',
      model={depVars={'ls', 'ha', 'dt'},
             effects={'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                      'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                      'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                      'v23', 'v24', 'v25', 'v26', 'v27'}},
      partitionByVar={name='Role', train='TRAIN', test='TEST'};
run;

The resulting output is shown in Output 18.1.3 through Output 18.1.5.

In Output 18.1.3, the "Model Information" table indicates that test set validation is used. The "Dimensions" table shows that the maximum number of factors to extract is nine. The "Number of Observations" table shows that nine sample observations are assigned to training roles and seven are assigned to testing roles.

Output 18.1.4 provides details about the results from test set validation. These results show that the absolute minimum PRESS is achieved with five extracted factors. Notice, however, that this is not much smaller than the PRESS for three factors. By using the cvTest parameter, you can perform the statistical model comparison that is suggested by Van der Voet (1994) to test whether this difference is significant, as shown in the following PROC CAS statements:

proc cas;
   action pls.pls /
      table='SAMPLE',
      method='PLS',
      model={depVars={'ls', 'ha', 'dt'},
             effects={'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                      'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                      'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                      'v23', 'v24', 'v25', 'v26', 'v27'}},
      partitionByVar={name='Role', train='TRAIN', test='TEST'},
      cvTest={pValue=0.15, seed=12345};
run;

The model comparison test is based on a rerandomization of the data. By default, the seed for this randomization is based on the system clock, but it is specified here. The resulting output is presented in Output 18.1.6 through Output 18.1.8.

The "Model Information" table in Output 18.1.6 displays information about the options that are used in the model comparison test. In Output 18.1.7, the p-value in comparing the test-set-validated residuals from models that have five and three factors indicates that the difference between the two models is insignificant; therefore, the model with fewer factors is preferred. The variation summary in Output 18.1.8 shows that more than 99% of the predictor variation and more than 90% of the response variation are accounted for by the three factors.

s:loadtable{casLib="casuser", path="Sample.csv"}

s:loadactionset{actionset='pls'}
r=s:pls{
   table='Sample',
   method='PLS',
   model={depVars={'ls', 'ha', 'dt'},
          effects={'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                   'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                   'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                   'v23', 'v24', 'v25', 'v26', 'v27'}} }

print(r.ModelInfo)
print(r.Dimensions)
print(r.NObs)
print(r.PercentVariation)
print(r.Timing)

s:loadactionset{actionset='pls'}
r=s:pls{
   table='Sample',
   method='PLS',
   model={depVars={'ls', 'ha', 'dt'},
          effects={'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                   'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                   'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                   'v23', 'v24', 'v25', 'v26', 'v27'}},
   partitionByVar={name='Role', train='TRAIN', test='TEST'} }

print(r.ModelInfo)
print(r.Dimensions)
print(r.NObs)
print(r["CrossValidation.ResidualSummary"])
print(r["CrossValidation.CVResults"])
print(r.PercentVariation)
print(r.Timing)

s:loadactionset{actionset='pls'}
r=s:pls{
   table='Sample',
   method='PLS',
   model={depVars={'ls', 'ha', 'dt'},
          effects={'v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                   'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                   'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                   'v23', 'v24', 'v25', 'v26', 'v27'}},
   partitionByVar={name='Role', train='TRAIN', test='TEST'},
   cvTest={pValue=0.15, seed=12345} }

print(r.ModelInfo)
print(r["CrossValidation.ResidualSummary"])
print(r["CrossValidation.CVResults"])
print(r.PercentVariation)

r

s.upload_file('Sample.csv')

s.loadactionset(actionset='pls')
r=s.pls(
   table='Sample',
   method='PLS',
   model={'depVars':['ls', 'ha', 'dt'],
          'effects':['v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                     'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                     'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                     'v23', 'v24', 'v25', 'v26', 'v27']} )

print(r.ModelInfo)
print(r.Dimensions)
print(r.NObs)
print(r.PercentVariation)
print(r.Timing)

s.loadactionset(actionset='pls')
r=s.pls(
   table='Sample',
   method='PLS',
   model={'depVars':['ls', 'ha', 'dt'],
          'effects':['v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                     'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                     'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                     'v23', 'v24', 'v25', 'v26', 'v27']},
   partitionByVar={'name':'Role', 'train':'TRAIN', 'test':'TEST'} )

print(r.ModelInfo)
print(r.Dimensions)
print(r.NObs)
print(r["CrossValidation.ResidualSummary"])
print(r["CrossValidation.CVResults"])
print(r.PercentVariation)
print(r.Timing)

s.loadactionset(actionset='pls')
r=s.pls(
   table='Sample',
   method='PLS',
   model={'depVars':['ls', 'ha', 'dt'],
          'effects':['v1', 'v2', 'v3', 'v4', 'v5', 'v6', 'v7', 'v8',
                     'v9', 'v10', 'v11', 'v12', 'v13', 'v14', 'v15',
                     'v16', 'v17', 'v18', 'v19', 'v20', 'v21', 'v22',
                     'v23', 'v24', 'v25', 'v26', 'v27']},
   partitionByVar={'name':'Role', 'train':'TRAIN', 'test':'TEST'},
   cvTest={'pValue':0.15, 'seed':12345} )

print(r.ModelInfo)
print(r["CrossValidation.ResidualSummary"])
print(r["CrossValidation.CVResults"])
print(r.PercentVariation)

print(r)