Example 42.3 Scatter Plot Smoothing by Selecting Spline Functions

This example shows how you can use model selection to perform scatter plot smoothing. It illustrates how you can use the experimental EFFECT statement to generate a large collection of B-spline basis functions from which a subset is selected to fit scatter plot data.

The data for this example come from a set of benchmarks developed by Donoho and Johnstone (1994) that have become popular in the statistics literature. The particular benchmark used is the "Bumps" functions to which random noise has been added to create the test data. The following DATA step, extracted from Sarle (2001), creates the data. The constants are chosen so that the noise-free data have a standard deviation of 7. The standard deviation of the noise is , yielding bumpsNoise with a signal-to-noise ratio of ().

   %let random=12345;
   
   data DoJoBumps;
      keep x bumps bumpsWithNoise;
   
      pi = arcos(-1);
   
      do n=1 to 2048;
         x=(2*n-1)/4096;
         link compute;
         bumpsWithNoise=bumps+rannor(&random)*sqrt(5);
         output;
      end;
   stop;
   
   compute:
      array t(11) _temporary_ (.1 .13 .15 .23 .25 .4 .44 .65 .76 .78 .81);
      array b(11) _temporary_ (   4    5    3   4   5 4.2 2.1 4.3  3.1  5.1  4.2);
      array w(11) _temporary_ (.005 .005 .006 .01 .01 .03 .01 .01 .005 .008 .005);
   
      bumps=0;
      do i=1 to 11;
         bumps=bumps+b[i]*(1+abs((x-t[i])/w[i]))**-4;
      end;
      bumps=bumps*10.528514619;
   return;
   run;

The following statements use the SGPLOT procedure to produce the plot in Output 42.3.1. The plot shows the bumps function superimposed on the function with added noise.

   proc sgplot data=DoJoBumps;
     yaxis display=(nolabel);
     series  x=x y=bumpsWithNoise/lineattrs=(color=black);
     series  x=x y=bumps/lineattrs=(color=red);
   run;

Output 42.3.1 Donoho-Johnstone Bumps Function

Suppose you want to smooth the noisy data to recover the underlying function. This problem is studied by Sarle (2001), who shows how neural nets can be used to perform the smoothing. The following statements use the LOESS statement in the SGPLOT procedure to show a loess fit superimposed on the noisy data (Output 42.3.2). (See Chapter 50, The LOESS Procedure, for information about the loess method.)

   proc sgplot data=DoJoBumps; 
      yaxis display=(nolabel);
      series x=x y=bumps;
      loess  x=x y=bumpsWithNoise / lineattrs=(color=red) nomarkers;
   run;

The algorithm selects a smoothing parameter that is small enough to enable bumps to be resolved. Because there is a single smoothing parameter that controls the number of points for all local fits, the loess method undersmooths the function in the intervals between the bumps.

Output 42.3.2 Loess Fit

Another approach to doing nonparametric fitting is to approximate the unknown underlying function as a linear combination of a set of basis functions. Once you specify the basis functions, then you can use least squares regression to obtain the coefficients of the linear combination. A problem with this approach is that for most data, you do not know a priori what set of basis functions to use. You need to supply a sufficiently rich set to enable the features in the data to be approximated. However, if you use too rich a set of functions, then this approach yields a fit that undersmooths the data and captures spurious features in the noise.

The penalized B-spline method (Eilers and Marx 1996) uses a basis of B-splines (see the section Constructed Effects and the EFFECT Statement of Chapter 18, Shared Concepts and Topics ) corresponding to a large number of equally spaced knots as the set of approximating functions. To control the potential overfitting, their algorithm modifies the least squares objective function to include a penalty term that grows with the complexity of the fit.

The following statements use the PBSPLINE statement in the SGPLOT procedure to show a penalized B-spline fit superimposed on the noisy data (Output 42.3.3). See Chapter 90, The TRANSREG Procedure, for details about the implementation of the penalized B-spline method.

   proc sgplot data=DoJoBumps; 
      yaxis display=(nolabel);
      series    x=x y=bumps;
      pbspline  x=x y=bumpsWithNoise / 
                  lineattrs=(color=red) nomarkers;
   run;

As in the case of loess fitting, you see undersmoothing in the intervals between the bumps because there is only a single smoothing parameter that controls the overall smoothness of the fit.

Output 42.3.3 Penalized B-spline Fit

An alternative to using a smoothness penalty to control the overfitting is to use variable selection to obtain an appropriate subset of the basis functions. In order to be able to represent features in the data that occur at multiple scales, it is useful to select from B-spline functions defined on just a few knots to capture large scale features of the data as well as B-spline functions defined on many knots to capture fine details of the data. The following statements show how you can use PROC GLMSELECT to implement this strategy:

   proc glmselect data=dojoBumps; 
      effect spl = spline(x / knotmethod=multiscale(endscale=8) 
                                split details);
      model bumpsWithNoise=spl;
      output out=out1 p=pBumps;
   run;
   
   proc sgplot data=out1; 
      yaxis display=(nolabel);
      series    x=x y=bumps;
      series    x=x y=pBumps / lineattrs=(color=red);
   run;
   

The KNOTMETHOD=MULTISCALE suboption of the EFFECT spl = SPLINE statement provides a convenient way to generate B-spline basis functions at multiple scales. The ENDSCALE=8 option requests that the finest scale use B-splines defined on equally spaced knots in the interval . Because the cubic B-splines are nonzero over five adjacent knots, at the finest scale, the support of each B-spline basis function is an interval of length about (5/256), enabling the bumps in the underlying data to be resolved. The default value is ENDSCALE=7. At this scale you will still be able to capture the bumps, but with less sharp resolution. For these data, using a value of ENDSCALE= greater than eight provides unneeded resolution, making it more likely that basis functions that fit spurious features in the noise are selected.

Output 42.3.4 shows the model information table. Since no options are specified in the MODEL statement, PROC GLMSELECT uses the stepwise method with selection and stopping based on the SBC criterion.