Updating old SAS programs: A case study with PROC LOESS

SAS programmers love to brag that the SAS will still run a program they wrote twenty or forty years. This is both a blessing and a curse. It's a blessing because it frees the statistical programmer from needing to revisit and rewrite code that was written long ago. It's a curse for the same reasons! Although revisiting an old SAS program requires time, it can be beneficial to revise an old program by incorporating newer features in SAS. The newer code can be shorter and more efficient.

I was reminded of this fact recently when I was looking at the SAS documentation of the LOESS procedure and noticed that two of the graphs in the documentation example use the Graph Template Language (GTL) instead of PROC SGPLOT. These examples were created for SAS 9.2 (in 2008). As part of the SAS 9.2 release, many documentation examples were updated to use the newer "SG" procedures instead of traditional graphing procedures (like GPLOT, G3D, and GCONTOUR). In SAS 9.2, PROC SGPLOT supported many basic features, but some advanced features were available only in GTL.

The current SAS release is SAS 9.4M8 (in 2023) with SAS 9.4M9 scheduled for 2025. In this article, I use modern features in PROC SGPLOT to remake two graphs in the PROC LOESS documentation. The old versions still run and look great, but the newer versions are simpler to produce and enable you to make modifications that improve the visualization of the data and the model.

Visualization of sulfate deposits

I will remake the graphs in a PROC LOESS example that is a spatial analysis of sulfate deposits at 179 locations in the contiguous US states. The doc example shows how to fit a loess smoother to these data, thus producing a model that predicts the mean sulfate deposit at any spatial location. The variables are the Longitude and Latitude of a location and the SO4 (sulfate, in g/m2) measurement at the location. The data are provided in the Appendix of this article. If you want to run the code in this article, you must first create the S04 data set.

After creating the S04 data set, the example uses the GTL to overlay two scatter plots. The first scatter plot uses the SYMBOL=CIRCLEFILLED option to create filled markers and uses the MARKERCOLORGRADIENT option to create a color scale that associates a color ramp with values of the S04 variable. The second scatter plot uses the same coordinates, but uses the SYMBOL=CIRCLE option. Thus, the result of the two SCATTERPLOT statements is a plot in which markers have a color-filled interior and a dark outline.

Today, this plot can be created by using a single statement in the SGPLOT procedure. First, you use the COLORRESPONSE=SO4 option to use the values of the SO4 variable to color the markers, then you use the FILLEDOUTLINEDMARKERS option to put a dark outline around the markers. Thus, the first graph can be remade by using the following call to PROC SGPLOT:

title "Sulfate Deposits in the US (1990)";
proc sgplot data=SO4;
scatter x=longitude y=latitude / colorresponse=SO4 markerattrs=(symbol=circleFilled size=12)
colormodel=ThreeColorRamp filledoutlinedmarkers;
run;

By looking at the location of the markers, you can see the outline of the contiguous United States. The left side of the graph shows measurements in California. The lower right is Florida. The upper right shows the New England states. You can see that many of the high values of SO4 are in or east of the industrial Midwest, including Ohio, Michigan, and Pennsylvania. Areas downwind of industrial emissions can show higher sulfate deposition due to the influence of acid rain.

Spatial smoothing of sulfate deposits

As you might expect, the main part of the example shows how to use PROC LOESS to create a loess smoother and to visualize the model. The example shows two separate calls to PROC LOESS. The first fits a model and creates a contour plot to visualize the model. (That call can be improved by using the newer OBS=OUTLINEGRADIENT option to add filled-and-outlined markers to the contour plot.) The second example shows how to create a scoring data set that consists of a regular grid of Longitude-Latitude locations. You can then use the SCORE statement in PROC LOESS to evaluate the regression model on the regular grid, in preparation for additional visualization. The following SAS code combines those two separate doc steps. It creates a regular grid on the range of the Longitude and Latitude coordinates, then scores the model at those locations.

/* Create a scoring data set on a regular grid.
Step 1: find range of Longitude and Latitude */
proc means data=SO4 min max;
var Longitude Latitude;
output out=MeansOut (keep=min_long min_lat max_long max_lat)
min(Longitude)=min_long min(Latitude)=min_lat
max(Longitude)=max_long max(Latitude)=max_lat;
run;

/* Step 2: construct regular grid from ranges */
data PredPoints;
num_intervals = 80;
set MeansOut;
d_long = (max_long - min_long) / num_intervals;
d_lat = (max_lat - min_lat) / num_intervals;
do Latitude = min_lat to max_lat by d_lat;
do Longitude = min_long to max_long by d_long;
output;
end;
end;
keep Latitude Longitude;
run;

/* run the loess fit and score the model on a uniform grid */
proc loess data=SO4 plots(only)=(ContourFit(OBS=OUTLINEGRADIENT));
model SO4=Longitude Latitude / degree=2 interp=cubic;
score data=PredPoints;
ods Output ScoreResults=ScoreOut; /* write predicted values to a data set */
run;

The contour plot is similar to the one shown in the documentation, but this one uses the OUTLINEGRADIENT option to improve the visualization.

I have two minor complaints about this visualization:

The plot doesn't display a gradient color ramp that shows the association of color and the values of SO4. The contour levels provide this information for the predicted values, but the graph does not enable you to determine the largest SO4 value for the markers, which is more than 4 g/m2.

The model predicts negative grams of sulfates off the coast of Florida and south of California. This is a well-known issue that can occur with any statistical model that extrapolates outside the data region. I have previously written about this issue and shown how you can modify the contour plot and the response surface plot to address the issue.

Switching from a contour plot to a heat map

The remainder of the PROC LOESS example uses GTL to create a response surface. PROC SGPLOT does not support response surfaces, so the GTL is the best way to create a response surface. However, PROC SGPLOT does support creating heat maps, which are often effective alternatives to response surfaces.

I'm not a big fan of response surfaces. I find that it is difficult to visualize the surface and the data at the same time. The hills and valleys of the surface can obstruct data markers. If you use an opaque surface (which is usually the case), no markers beneath the surface are visible. Heat maps (or contour plots) do not suffer from these issues. Since PROC LOESS creates a contour plot automatically, let's try to overlay a heat map and a scatter plot by using PROC SGPLOT.

In the heat map, we can get rid of any locations at which the model predicts a negative value. Let's change all negative predicted values to a missing value. We can then use the HEATMAPPARM statement in PROC SGPLOT to visualize the predictions from the model:

data ScoreOut_noneg;
set ScoreOut;
if p_SO4 < 0 then p_SO4=.;
run;

title "Predicted Mean Sulfate Deposits";
proc sgplot data=ScoreOut_noneg;
heatmapparm x=Longitude y=Latitude colorresponse=p_SO4;
run;

You can see grey areas, which represent missing values, in the lower-left and lower-right of the graph. This changes the lowest value for the color ramp. You can overlay the observed and predicted values by merging the data with the predicted values, as follows:

data all;
set ScoreOut_nomiss SO4(rename=(Longitude=x Latitude=y));
run;

title "Predicted and Observed Sulfate Deposits";
proc sgplot data=all;
heatmapparm x=Longitude y=Latitude colorresponse=p_SO4 / colormodel=ThreeColorRamp;
scatter x=x y=y / colorresponse=SO4 markerattrs=(symbol=circleFilled)
colormodel=ThreeColorRamp filledoutlinedmarkers;
run;

This combined heat map and scatter plot contains a lot of information. I think it is a better choice than a surface plot. I like that the observed and predicted values use a common color model.

Summary

This article looks at a documentation example that was written for SAS 9.2. The example can be simplified and improved by using features of SAS 9.4 (M5 or above). The examples from SAS 9.2 still run and provide an effective visualization, but the GTL in the example can be eliminated in favor of newer statements and options in PROC SGPLOT.

While writing this article, I was reminded of my former colleague, Warren Repole, who used to give conference presentations titled, "Don't Be a SAS Dinosaur," in which he would show SAS programmers newer features in SAS software that would make their programs easier to write and maintain. In some organizations, the decision about whether to update old SAS programs is a complicated one. On the one hand, if your old SAS program still runs, why mess with it? On the other hand, if your program was written long ago, there might be newer features in SAS 9.4 that can simplify your SAS program or even improve its effectiveness at analyzing or visualizing data or statistical models.

Appendix: The sulfate data

data SO4;
input Latitude Longitude SO4 @@;
/* Longitudes decrease from west to east in the western hemisphere,
so use negative values when graphing longitudes in the western hemisphere. */
Longitude = -Longitude;
format Latitude f4.0;
format Longitude f4.0;
format SO4 f4.1;
label S04 = "Sulfate (g / m^2)";
datalines;
32.45833 87.24222 1.403 34.28778 85.96889 2.103
33.07139 109.86472 0.299 36.07167 112.15500 0.304
31.95056 112.80000 0.263 33.60500 92.09722 1.950
34.17944 93.09861 2.168 36.08389 92.58694 1.578
36.10056 94.17333 1.708 39.00472 123.08472 0.096
36.56694 118.77722 0.259 41.76583 122.47833 0.065
34.80611 119.01139 0.053 38.53528 121.77500 0.135
37.44139 105.86528 0.247 38.11778 103.31611 0.326
39.99389 105.48000 0.687 39.40306 107.34111 0.225
39.42722 107.37972 0.339 37.19806 108.49028 0.559
40.50750 107.70194 0.250 40.36417 105.58194 0.307
40.53472 106.78000 0.564 37.75139 107.68528 0.557
39.10111 105.09194 0.371 40.80639 104.75472 0.286
29.97472 82.19806 1.279 28.54278 80.64444 1.564
25.39000 80.68000 0.912 30.54806 84.60083 1.243
27.38000 82.28389 0.991 32.14111 81.97139 1.225
33.17778 84.40611 1.580 31.47306 83.53306 0.851
43.46139 113.55472 0.180 43.20528 116.74917 0.103
44.29778 116.06361 0.161 40.05333 88.37194 2.940
41.84139 88.85111 2.090 41.70111 87.99528 3.171
37.71000 89.26889 2.523 38.71000 88.74917 2.317
37.43556 88.67194 3.077 40.93333 90.72306 2.006
40.84000 85.46389 2.725 38.74083 87.48556 3.158
41.63250 87.08778 3.443 40.47528 86.99222 2.775
42.90972 91.47000 1.154 40.96306 93.39250 1.423
37.65111 94.80361 1.863 39.10222 96.60917 0.898
38.67167 100.91639 0.536 37.70472 85.04889 2.693
37.07778 82.99361 2.195 38.11833 83.54694 2.762
36.79056 88.06722 2.377 29.92972 91.71528 1.276
30.81139 90.18083 1.393 44.37389 68.26056 2.268
46.86889 68.01472 1.551 44.10750 70.72889 1.631
45.48917 69.66528 1.369 39.40889 76.99528 2.535
38.91306 76.15250 2.477 41.97583 70.02472 1.619
42.39250 72.34444 2.156 42.38389 71.21472 2.417
45.56083 84.67833 1.701 46.37417 84.74139 1.539
47.10472 88.55139 1.048 42.41028 85.39278 3.107
44.22389 85.81806 2.258 47.53111 93.46861 0.550
47.94639 91.49611 0.563 46.24944 94.49722 0.591
44.23722 95.30056 0.604 32.30667 90.31833 1.614
32.33472 89.16583 1.135 34.00250 89.80000 1.503
38.75361 92.19889 1.814 36.91083 90.31861 2.435
45.56861 107.43750 0.217 48.51028 113.99583 0.387
48.49917 109.79750 0.100 46.48500 112.06472 0.209
41.15306 96.49278 0.743 41.05917 100.74639 0.391
36.13583 115.42556 0.139 41.28528 115.85222 0.075
38.79917 119.25667 0.053 39.00500 114.21583 0.273
43.94306 71.70333 2.391 40.31500 74.85472 2.593
33.22028 108.23472 0.377 35.78167 106.26750 0.315
32.90944 105.47056 0.355 36.04083 106.97139 0.376
36.77889 103.98139 0.326 42.73389 76.65972 3.249
42.29944 79.39639 3.344 43.97306 74.22306 2.322
44.39333 73.85944 2.111 41.35083 74.04861 3.306
43.52611 75.94722 3.948 42.10639 77.53583 2.231
41.99361 74.50361 3.022 36.13250 77.17139 1.857
35.06056 83.43056 2.393 35.69694 80.62250 2.082
35.02583 78.27833 1.729 34.97083 79.52833 1.959
35.72833 78.68028 1.780 47.60139 103.26417 0.354
48.78250 97.75417 0.306 47.12556 99.23694 0.273
39.53139 84.72417 3.828 40.35528 83.06611 3.401
39.79278 81.53111 3.961 40.78222 81.92000 3.349
36.80528 98.20056 0.603 34.98000 97.52139 0.994
36.59083 101.61750 0.444 44.38694 123.62306 0.629
44.63472 123.19000 0.329 45.44917 122.15333 0.716
43.12167 121.05778 0.050 44.21333 122.25333 0.423
43.89944 117.42694 0.071 45.22444 118.51139 0.109
40.78833 77.94583 3.275 41.59778 78.76750 4.336
40.65750 77.93972 3.352 41.32750 74.82028 3.081
33.53944 80.43500 1.456 44.35500 98.29083 0.372
43.94917 101.85833 0.224 35.96139 84.28722 3.579
35.18250 87.19639 2.148 35.66444 83.59028 2.474
35.46778 89.15861 1.811 33.95778 102.77611 0.376
28.46667 97.70694 0.886 29.66139 96.25944 0.934
30.26139 100.55500 0.938 32.37861 94.71167 2.229
31.56056 94.86083 1.472 33.27333 99.21528 0.890
33.39167 97.63972 1.585 37.61861 112.17278 0.237
41.65833 111.89694 0.271 38.99833 110.16528 0.143
41.35750 111.04861 0.172 42.87611 73.16333 2.412
44.52833 72.86889 2.549 38.04056 78.54306 2.478
37.33139 80.55750 1.650 38.52250 78.43583 2.360
47.86000 123.93194 1.144 48.54056 121.44528 0.837
46.83528 122.28667 0.635 46.76056 117.18472 0.255
37.98000 80.95000 2.396 39.08972 79.66222 3.291
45.79639 88.39944 1.054 45.05333 88.37278 1.457
44.66444 89.65222 1.044 43.70194 90.56861 1.309
46.05278 89.65306 1.132 42.57917 88.50056 1.809
45.82278 91.87444 0.984 41.34028 106.19083 0.335
42.73389 108.85000 0.236 42.49472 108.82917 0.313
42.92889 109.78667 0.182 43.22278 109.99111 0.161
43.87333 104.19222 0.306 44.91722 110.42028 0.210
45.07611 72.67556 2.646
;