Statistical programming in SAS with an emphasis on SAS/IML programs

In statistics, distances between observations are used to form clusters, to identify outliers, and to estimate distributions. Distances are used in spatial statistics and in other application areas.

There are many ways to define the distance between observations. I have previously written an article that explains Mahalanobis distance, which is used often in multivariate analysis, and I have showed how to compute the Mahalanobis distance in SAS. Today's article is simpler: how do you compute the usual Euclidean distance in SAS?

Recall that the squared Euclidean distance between the point p = (p₁, p₂, ..., p_n) and the point q = (q₁, q₂, ..., q_n) is the sum of the squares of the differences between the components: Dist²(p,q) = Σ_i (p_i – q_i)². The Euclidean distance is then the square root of Dist²(p,q). This article shows three ways to compute the Euclidean distance in SAS:

1. By using the DISTANCE procedure in SAS/STAT software.
2. By using the DISTANCE function in SAS/IML software.
3. By writing your own SAS/IML function.

The following DATA step creates 27 observations that are arranged on an integer lattice in three dimensions. Each row of the data set contains an observation in three variables. The 27 points are (0,0,0), (0,0,1), (0,0,2), (0,1,0), ..., (2,2,2).

data Obs;
do x=0 to 2;
   do y=0 to 2;
      do z = 0 to 2;
         output;
      end;
   end;
end;
run;

Compute distance by using SAS/STAT software

PROC DISTANCE can compute many kinds of distance, and can also standardize the data variables, which is useful when your variable represent different quantities (such as height, weight, and age). In the simple case of Euclidean distance without any standardization, specify the METHOD=EUCLID option and the NOSTD option on the PROC DISTANCE statement, as follows:

proc distance data=Obs out=Dist method=Euclid nostd;
   var interval(x y z);
run;
 
proc print data=Dist(obs=4);
   format Dist: 8.6;
   var Dist1-Dist4;
run;

The output data set has 27 variables and 27 observations. The preceding table shows the first four observations and the first four variables. You can see that the output data set is the lower-triangular portion of the distance matrix. The ith row gives the distance between the ith observation and the jth observation for j ≤ i. For example, the distance between the fourth observation (0,1,0) and the second observation (0,0,1) is sqrt(0² + 1² + 1²)= sqrt(2) = 1.414.

If you prefer to output the full, dense, symmetric matrix of distances, use the SHAPE=SQUARE option on the PROC DISTANCE statement.

Computing distance in SAS/IML Software

In SAS/IML software, you can use the DISTANCE function in SAS/IML to compute a variety of distance matrices. The DISTANCE function was introduced in SAS/IML 12.1 (SAS 9.3M2).

By default, the DISTANCE function computes the Euclidean distance, and the output is always a square matrix. Therefore, the following statements compute the Euclidean pairwise distances between the 27 points in the Obs data set:

proc iml;
use Obs;
read all var _NUM_ into X;
close Obs;
 
D = distance(X);
print (D[1:4, 1:4])[format=8.6 c=("Dist1":"Dist4")];

The output shows that the values in the upper-left portion of the distance matrix are the same as were computed by PROC DISTANCE.

A SAS/IML Module for computing Euclidean distance

Sometimes you do not want to compute the complete matrix of pairwise distances between observations. Sometimes you only need the distances between observations that belong to different groups. For example, you might have one set of locations that represent warehouses and another set that represents the locations of retail stores. You might be interested in computing the distances from each store to the warehouses, so that you can efficiently ship goods to each store from the closest warehouse.

The following SAS/IML module computes the distances between observations in the matrix X and the matrix Y. The (i,j)th element of the result is the distance between the ith row of X and the jth row of Y.

/* compute Euclidean distance between points in x and points in y.
   x is a n x d matrix, where each row is a point in d dimensions.
   y is a m x d matrix.
   The function returns the n x q matrix of distances, D, such that
   D[i,j] is the distance between x[i,] and y[j,]. */
start PairwiseDist(x, y);
   if ncol(x)^=ncol(y) then return (.);       /* different dimensions */
   n = nrow(x);  m = nrow(y);
   idx = T(repeat(1:n, m));                   /* index matrix for x   */
   jdx = shape(repeat(1:m, n), n);            /* index matrix for y   */
   diff = X[idx,] - Y[jdx,];
   return( shape( sqrt(diff[,##]), n ) );     /* sqrt(sum of squares) */
finish;
 
p = { 1  0, 
      1  1,
     -1 -1};
q = { 0  0,
     -1  0};
PD = PairwiseDist(p,q); 
print PD;

If you are still running SAS/IML 9.3 or earlier, you can use the PairwiseDist function to define your own function that computes the Euclidean distance between rows of a matrix:

/* compute Euclidean distance between points in x.
   x is a p x d matrix, where each row is a point in d dimensions.
   Use the DISTANCE function in SAS/IML 12.1 and later releases. */
start EuclideanDistance(x);  
   y=x;
   return( PairwiseDist(x,y) );
finish;
 
D = EuclideanDistance(X);
print (D[1:4, 1:4])[format=8.6 c=("Dist1":"Dist4")];

The output is the same as for the built-in DISTANCE function, and is not printed.

Someone recently asked a question on the SAS Support Communities about estimating parameters in ridge regression. I answered the question by pointing to a matrix formula in the SAS documentation. One of the advantages of the SAS/IML language is that you can implement matrix formulas in a natural way. The SAS/IML expressions can often be written almost verbatim from the formula. This article uses a ridge regression formula from the PROC REG documentation to illustrate this feature.

When to use ridge regression?

Ridge regression is a variant to least squares regression that is sometimes used when several explanatory variables are highly correlated. The "usual" ordinary least squares (OLS) regression produces unbiased estimates for the regression coefficients (in fact, the Best Linear Unbiased Estimates). However, when the explanatory variables are correlated, the OLS parameter estimates have large variance. It might be desirable to use a different regression technique, such as ridge regression, in order to obtain parameter estimates that have smaller variance. The trade-off is that the estimates for the ridge regression method are biased.

If X is the centered and scaled data matrix, then the crossproduct matrix X`X is nearly singular when columns of X are highly correlated. Mathematically, ridge regression adds a multiple (called the ridge parameter, k) of the identity matrix to the X`X matrix. The nonsingular matrix (X`X + kI) is then used to compute the parameter estimates.

Each value of k results in a different set of parameter estimates. There have been many papers written about how to choose the best value of k. In this article, I merely want to show how to compute the parameters for ridge regression when the ridge parameter is given.

How to compute ridge regression in SAS

In SAS software, you can compute ridge regression by using the REG procedure. The following example from the PROC REG documentation is used to illustrate ridge regression. The RIDGE= option specifies the value(s) of the ridge parameter, k. The OUTEST= option is used to create an output data set that contains the parameter estimates for each value of k.

/* Ridge regression example from PROC REG documentation */
data acetyl;
input x1-x4 @@;
x1x2 = x1 * x2;
x1x1 = x1 * x1;
datalines;
1300  7.5 0.012 49   1300  9   0.012  50.2 1300 11 0.0115 50.5
1300 13.5 0.013 48.5 1300 17   0.0135 47.5 1300 23 0.012  44.5
1200  5.3 0.04  28   1200  7.5 0.038  31.5 1200 11 0.032  34.5
1200 13.5 0.026 35   1200 17   0.034  38   1200 23 0.041  38.5
1100  5.3 0.084 15   1100  7.5 0.098  17   1100 11 0.092  20.5
1100 17   0.086 29.5
;
 
proc reg data=acetyl outest=b ridge=0.02 noprint;
   model x4=x1 x2 x3 x1x2 x1x1;
run;
 
proc print data=b(where=(_TYPE_="RIDGE")) noobs;
   var _RIDGE_ Intercept x1 x2 x3 x1x2 x1x1;
run;

The parameter estimates for the ridge regression are shown for the ridge parameter k = 0.02.

Implementing a matrix formula for ridge regression by using SAS/IML software

The question that was asked on the SAS Discussion Forum was about where to find the matrix formula for estimating the ridge regression coefficients. The documentation for the REG procedure includes a section that provides the formula that PROC REG uses to compute the ridge regression coefficients. The accompanying text says:

Let X be the matrix of the independent variables after centering [and scaling] the data, and let Y be a vector corresponding to the [centered] dependent variable. Let D be a diagonal matrix with diagonal elements as in X`X. The ridge regression estimate corresponding to the ridge constant k can be computed as D<sup>-1/2</sup>(Z`Z + kI)^-1Z`Y.

The terms in brackets do not appear in the original documentation, but I included them for clarity. Since this is a matrix formula, let's use the SAS/IML language to implement the formula. The following SAS/IML program uses the formula to compute the ridged parameter estimates:

/* Ridge regression coefficients computed in SAS/IML */
proc iml;
use acetyl;        
read ALL var {x1 x2 x3 x1x2 x1x1} into X; 
read all var {x4} into y;
close acetyl;        
 
/* ridge regression */
xBar = mean(X);       /* save row vector of means */
yBar = mean(y);       /* save mean of Y */
X = X - xBar;         /* center X and y; do not add intercept */
y = y - yBar;                      
D = vecdiag(X`*X);
Z = X / sqrt(D`);     /* Z`Z = corr(X) */
k = 0.02;             /* choose a ridge parameter */
b =  (1/sqrt(D)) # inv(Z`*Z + k*I(ncol(X))) * (Z`*y); /* formula in REG doc */
print (b`)[colname={"x1" "x2" "x3" "x1x2" "x1x1"}];

So that the formula and the SAS/IML statements match each other, I have written the computation in a "natural" way, rather than in a more efficient way. See the article "Do you really need to compute that matrix inverse?" In any case, you can see that the parameter estimates that are compute from the formula agree with the estimates that are computed by PROC REG.

However, notice that the PROC IML computations do not include an intercept term. This is because the independent and dependent variables were all centered, so the intercept in these coordinates is exactly zero. To obtain the intercept in the original (uncentered) coordinates, you can use a simple formula that recovers the intercept estimate:

   /* The ridge regression uses centered x and y. Recover the intercept:
      y-yBar = b1(x1-x1Bar) + ... + bn(xn-xnBar)
   so      y = yBar-Sum(bi*xibar) + b1 x1 + ... + bn xn
   Consequently, b0 = yBar-Sum(bi*xibar). */
   b0 = ybar - xbar * b;
   print b0[label="Intercept"];

Got a matrix formula? Use SAS/IML software

As this article shows, the SAS/IML language enables you to implement matrix formulas in a natural way. Although this article implements a formula that is already built into a SAS procedure, the same ideas apply for formulas and algorithms that are not available in any SAS procedure.

So the next time that you need to implement a matrix formula that appears in a textbook, in documentation, or in a journal article, reach for the SAS/IML language!