Statistical Graphics Using ODS

Polynomial splines have interior knots (knots that are inside the range of the minimum and maximum x). B-splines have interior knots and exterior knots (knots that are outside the range of the minimum and maximum x). For example, if x ranges from –2 to 2, a reasonable knot vector for three interior knots is , where is a small value such as . These knots are all but evenly spaced. The addition and subtraction of ensure that the first three exterior knots are less than the minimum and the last three exterior knots are greater than the maximum. Statistical procedures such as PROC TRANSREG and the procedures that support the EFFECT statement enable you to explicitly control both interior and exterior knots. Simpler options enable you to specify the number of knots and let the procedure determine their values. In many cases, fit functions based on slightly different knots are only slightly different. However, when you fit a model to one set of data and then score another set of data, you must ensure that you use the same knots for both purposes. The following procedure steps all fit the same model and use the knot vector :

data x2;
   do x = -2 to 2 by 0.1;
      y = x + sin(x) + normal(17);
      output;
   end;
run;

proc transreg data=x2 details ss2;
   model identity(y) = bspline(x / nknots=3 evenly=3);
run;

proc transreg data=x2 details ss2;
   model identity(y) = spline(x / nknots=3 evenly=3);
run;

%let k = -4.000000000001 -3.000000000001 -2.000000000001 -1 0 1
          2.000000000001  3.000000000001  4.000000000001;
proc transreg data=x2 details ss2;
   model identity(y) = bspline(x / knots=&k);
run;

proc orthoreg data=x2;
   effect spl = spline(x / knotmethod=equal(3));
   model  y = spl / noint;
run;

proc orthoreg data=x2;
   effect spl = spline(x / knotmethod=listwithboundary(&k));
   model  y = spl / noint;
run;

Although the models are the same, how the effects are grouped can vary. PROC TRANSREG with the DETAILS option displays the knots shown in Output 24.6.25. PROC TRANSREG with the BSPLINE expansion fits a model that has no intercept and seven 1-df B-spline effects. PROC TRANSREG with the SPLINE transformation fits a model that has an intercept and one 6-df B-spline effect. PROC ORTHOREG fits a model that has no intercept and seven 1-df B-spline effects.

PROC TRANSREG, because it has more than a 30-year history, does not create evenly spaced knots by default. By default, in the data set X2, the knot list is . (Exterior knots are repeated.) Both repeated and all-but-equally-spaced exterior knots work, but the latter method better matches later software developments.

The following four steps show four equivalent methods of putting the B-spline basis into a SAS data set:

proc transreg data=x2 details ss2;
   model identity(y) = bspline(x / nknots=3 evenly=3);
   output out=b1(keep=x_:) replace;
run;

proc glimmix data=x2 outdesign=b2(drop=x y
   rename=(%macro ren; %do i = 1 %to 7; _x&i=x_%eval(&i-1) %end; %mend; %ren));
   effect spl = spline(x / knotmethod=equal(3));
   model y = spl / noint;
run;

%let k = -4.000000000001 -3.000000000001 -2.000000000001 -1 0 1
          2.000000000001  3.000000000001  4.000000000001;

proc iml;
   use x2(keep=x); read all into x;
   b = bspline(x, 3, {&k});
   vname = 'x_0' : rowcatc('x_' || char(ncol(b)-1));
   create b3 from b[colname=vname]; append from b;
quit;

data b4(keep=x:);
   %let d   = 3;                        /* degree                  */
   %let nkn = %eval(&d * 2 + 3);        /* total number of knots   */
   %let nb  = %eval(&d + 1 + 3);        /* number of cols in basis */
   array k[&nkn] _temporary_ (&k);      /* knots                   */
   array b[&nb] x_0 - x_%eval(&nb - 1); /* basis                   */
   array w[%eval(2 * &d)];              /* work                    */
   set x2;
   do i = 1 to &nb; b[i] = 0; end;

   * find the index of first knot greater than current data value;
   do ki = 1 to &nkn while(k[ki] le x); end;
   kki = ki - &d - 1;

   * make the basis;
   b[1 + kki] = 1;
   do j = 1 to &d;
      w[&d + j] = k[ki + j - 1] - x;
      w[j] = x - k[ki - j];
      s = 0;
      do i = 1 to j;
         t = w[&d + i] + w[j + 1 - i];
         if t ne 0.0 then t = b[i + kki] / t;
         b[i + kki] = s + w[&d + i] * t;
         s = w[j + 1 - i] * t;
      end;
      b[j + 1 + kki] = s;
   end;
run;

The %REN macro is used to rename the PROC GLIMMIX variables, because different procedures have different rules for naming variables. Each data set has seven columns (intercept, degree 3, and three knots, although no column corresponds directly to these terms). The B-spline basis has at most order 4 (degree + 1) nonzero values in each row. This sparseness can be exploited by modeling code, which is one of the many appealing characteristics of B-splines. The following steps compare and display the bases:

options nolabel;
proc compare error note briefsummary criterion=1e-12
   data=b1 compare=b2 method=relative(1);
run;

proc compare error note briefsummary criterion=1e-12
   data=b1 compare=b3 method=relative(1);
run;

proc compare error note briefsummary criterion=1e-12
   data=b1 compare=b4(drop=x) method=relative(1);
run;
options label;

proc sgplot data=b4;
   %macro s; %do i = 0 %to &nb-1; series y=x_&i x=x; %end; %mend; %s
   yaxis label='B-Spline Basis Functions';
run;

The PROC COMPARE output (not displayed here) shows that the matrices are identical. The B-spline basis is plotted in Output 24.6.26.

Output 24.6.26: B-Spline Basis

The next steps show that the polynomial-spline basis and the B-spline basis are linear transformations of each other:

proc transreg data=x2 details ss2;
   model identity(y) = pspline(x / name=(p) nknots=3);
   output out=p1(keep=p_:) replace;
run;

data both; merge p1 b1; run;

proc print data=both noobs;
   format _numeric_ bestd6.5 p_1 p_2 best6.;
run;

proc cancorr data=both;
  var x:;
  with p:;
run;

The polynomial spline (variables that begin with p_) and B-spline (variables that begin with x_) bases are displayed in Output 24.6.27. The variable p_1 matches the original X variable. All canonical correlations in Output 24.6.28 are 1.0. This shows that every vector in the space that is defined by the columns in the B-spline basis can be constructed from a linear combination of the columns of the polynomial-spline basis.