Nevertheless, many SAS customers use macro constructs as part of their daily SAS programming tasks, and that practice often continues when they write SAS/IML programmers. A customer recently asked a question about the macro language that required knowledge of the way that macro variables are handled within a SAS/IML loop. This post shares my response.

Here's the crux of the customer's question. Run the following SAS/IML program and see if you can understand why it behaves as it does:

proc iml;
i = 7;
call symputx("j", i);    /* 1. Put value of i into macro variable j */
y1 = &j;                 /* 2. Assign y1 the value of &j            */
print y1;                /* success! */
 
y = j(1,4,.);
do i = 1 to ncol(y);     /* 3. Start processing the DO block of statements */
   call symputx("j", i); /* 4. Put value of i into macro variable j */
   y[i] = &j;            /* 5. Hmmmm, what does this do inside the loop? */
end;
print y;                 /* Not what you might expect? */

As you can see from the output, the first use of the macro variable (outside the DO loop), works as expected. But the second does not. The customer wanted to know why the elements of y are not set to 1, 2, 3, 4 within the loop.

The key point to remember about macro variables is that SAS code never sees them. Macro variables are evaluated by the macro preprocessor at parse time, not at run time. The SAS/IML code never sees &j, only the constant value that the preprocessor substitutes for &j.

It is also important to remember that PROC IML is an interactive procedure. (The "I" in IML stands for interactive!) Each statement or block of statements is parsed as it is encountered, as opposed to the DATA step, which parses the entire program before beginning execution.

Let's examine the program step-by-step to understand why the first construct works but the second does not. The following steps refer to the numbers in the program comments:

1. The value of the SAS/IML scalar i is copied (as text) into the macro variable j.
2. The statement is encountered. The value of the macro variable j is substituted by the macro preprocesser. Then the statement is executed. The SAS/IML variable y1 is assigned to the value 7.
3. A DO loop is encountered by the SAS/IML parser. The parser finds the matching END statement and proceeds to parse the entire body of the loop in order to check for syntax errors. This parsing phase occurs exactly one time. Because the block of statements contain a macro variable, the macro preprocessor substitutes the value of the macro variable j, which is 7.
4. For each iteration, the value of the SAS/IML scalar i is copied (as text) into the macro variable j.
5. For each iteration, the ith element of the y vector is assigned the value 7. In particular, this statement does not contain a reference to the macro varible j.

To the casual reader of the program, it looks like &j will have a different value during each step of the iteration. But but it doesn't. The expression &j is resolved at parse time. SAS/IML parses the entire body of the DO loop once, before any execution occurs, and at parse time the expression &j is 7.

There is a way to get what the customer wants. The SYMGET function retrieves the value of a macro variable at run time. Therefore the following statements fill the vector y with the values 1 through 4:

do i = 1 to ncol(y);
   call symputx("j", i);
   y[i] = num(symget("j"));  /* get macro value at run time */
end;
print y;                     /* Yes! This is what we want! */

For me, this blog post emphasizes three facts:

• Always remember that macro substitution is done by a preprocessor, which operates at parse time.
• The SAS/IML language parses an entire block of statements (between the DO and END statements) one time before executing the block.
• Mixing macro code and SAS/IML statements can be confusing and hard to debug. When you have the option, use SAS/IML language features instead of relying on macro language constructs.

For example, the following SAS/IML statement defines a vector for which the value 2.2 is repeated two times and the value 3.3 is repeated three times. The resulting vector has five elements:

proc iml;
x = {[2] 2.2  [3] 3.3};

This vector is constructed as a matrix literal. If instead you have the values and frequencies in separate vectors, then use the ExpandFreq module in my previous post:

values = {2.2 3.3};
freq   = {2   3};
load module=(ExpandFreq); /* define or load ExpandFreg module here */
y = ExpandFreq(values, freq);

This is probably a good time to remind everyone about the SAS/IML Support Community. You can post your SAS/IML questions there 24 hours a day. That is always better than sending me direct email. There are lots of experienced SAS/IML experts out there, so use the SAS/IML Support Community to tap into that knowledge.

In a previous blog post, I described how to use a spread plot to compare the distributions of several variables. Each spread plot is a graph of centered data values plotted against the estimated cumulative probability. Thus, spread plots are similar to a (rotated) plot of the empirical cumulative distribution function. Users of the SAS regression procedures will recognize the spread plots as one of the plots that are created automatically by procedures such as PROC REG. The spread plots of the fitted and residual values appear in the middle column of the third row of the regression diagnostics panel.

In the SAS documentation, the residual-fit spread plot is also called an "RF plot." This article describes how to interpret the R-F spread plot.

The residual-fit spread plot in SAS output

When I first saw the R-F spread plot in the PROC REG diagnostics panel, there were two things that I found confusing:

• The title of the left plot is "Fit–Mean." I read the title as "fit hyphen mean," and I didn't know what that meant. However, the correct way to read the title is "fit minus mean," which is equivalent to "centered fit."
• The label for the horizontal axis is "Proportion Less." I didn't know what that meant, either. I now know that scatter plot shows empirical quantiles versus their plotting positions. Recall that the pth empirical quantile is the data value that is greater than (or equal to) a proportion p of the data values. Therefore, if a point on the scatter plot has coordinates (p_i, q_i), it means that the vertical coordinate is the ith quantile, and approximately p_i of the other data values are less than that proportion.

History of the residual-fit spread plot

The spread plot is a graph of the centered data versus the corresponding plotting position. Essentially, it is a plot of the sorted data against the corresponding rank, except that using the plotting position instead of ranks makes it possible to compare variables that have different numbers of nonmissing observations. Also, using centered data instead of raw values enables you to compare the spread of variables that have different means.

William S. Cleveland featured the residual-fit spread plot in his book Visualizing Data (1993). He describes how to create a quantile plot on pp. 17–20, and describes quantile plots for fitted and residual values on p. 35–38. Then he says (p. 40):

It is informative to study how influential the [explanatory] variable is in explaining the variation in the [response variable]. The fitted values and the residuals are two sets of values each of which has a distribution. If the spread of the fitted-value distribution is large compared with the spread of the residual distribution, then the [explanatory] variable is influential. If it is small, the [explanatory] variable is not as influential.... Since it is the spreads of the distributions that are of interest, the fitted values minus their overall mean are graphed.... This residual-fit spread plot, or r-f spread plot, shows [whether] the spreads of the residuals and fit values are comparable.

Cleveland goes on to use the R-F spread plot about 20 times in multiple examples.

The residual-fit spread plot as a regression diagnostic

Following Cleveland's examples, the residual-fit spread plot can be used to assess the fit of a regression as follows:

• Compare the spread of the fit to the spread of the residuals. This is the main idea. If the left side of the plot (the centered fitted values) is taller than the right side (the residual values), then you conclude that the spread of the residual values is small relative to the spread of the fitted values.
• Examine the distribution of the residual values. The quantile plot of the residual values contains all of the information that a box plot does—and more. If the distribution does not appear to be normally distributed, the model might not fit the data.
• Are there extreme values for the distribution of the residual values? These indicate outliers: observations for which the observed value is far from the fitted value.
• Are there extreme values for the distribution of the fitted values? These might indicate influential observations that have high leverage in the model. They need to be investigated.

Some examples on simulated data

The best way to practice interpreting the R-F spread plots are to view some examples for which the true model is known. The following DATA step simulates two response variables:

data RegData(drop=i);
call streaminit(12345);
do i = 1 to 100;
   x = rand("Normal");
   y1 = 2 + 4*x    + rand("Normal", 0, 0.25); /* small error */
   y2 = 2 + 4*x**2 + rand("Normal", 0, 1);    /* not linear in x */
   output;
end;
run;

For a real regression analysis, I would look at several diagnostic plots, but in the subsequent examples I will only interpret the residual-fit spread plots. I use the DIAGNOSTICS(UNPACK) option on the PLOTS= option to extract the R-F spread plot from the diagnostics panel.

Example 1: Examining the residual variation in a model

The y1 variable has a small error term. The following statements display the R-F spread plot:

title "y = 2 + 4*x + eps, eps ~ N(0,0.25)";
ods select RFPlot;
proc reg data = RegData plots=diagnostics(unpack);
   model y1 = x;
quit;

Notice that the left plot (the centered fitted values) is "taller" than the right plot (the residual values), which indicates that the residual values have a smaller spread. In terms of the model, the x variable accounts for a significant portion of the variation in the model, with only a little residual variation.

You can change the standard deviation of the error term to 5 and rerun the program. The new R-F spread plot (not shown), shows that the spread of the residual values is larger than the spread of the fitted values. The interpretation would be that considerable variation remains after accounting for the effect of the x variable.

Example 2: A misspecified model

In the previous example, the model was correctly specified. In this second example, the true model is quadratic in the x variable, but the fitted model is linear in x.

title "y = 2 + 4*x**3 + eps, eps ~ N(0,0.25)";
title2 "Model is y = x";
ods select RFPlot;
proc reg data = RegData plots(only)=diagnostics(unpack);
   model y2 = x;
quit;

In the R-F spread plot for the (misspecified) model, the right-hand plot is taller than the left-hand plot. This shows that there is a lot of variation that is not explained by the model. Furthermore, the residual distribution does not appear to be normally distributed. The right tail of the residual distribution is long, and the distribution is skewed. If I saw a plot like this in real data, I would look at other plots (such as the plot of residuals versus the predicted values) to see if the residuals exhibit a pattern that can be modeled.

Closing Comments

Several SAS regression procedures produce a regression diagnostics panel automatically. Each graph reveals information about the regression model and whether it fits the data well. This article has described how to interpret a residual-fit plot, which is located in the last row of the diagnostics panel. The residual-fit spread plot, which was featured prominently in Cleveland's book, Visualizing Data, is one tool in the arsenal of regression diagnostic plots.

The usual way to compare data distributions is to use histograms. One technique is to display a panel of histograms, which are known as comparative histograms. I have used this approach to compare salaries between two categories of workers. The comparative histogram is produced automatically by PROC UNIVARIATE when the analysis includes a classification variable.

A related approach is to use transparency to overlay two histograms on the same axis. This method works well for two distributions. However, for three or more distributions, an overlay of histograms can be difficult to read.

A problem of scale

If one of the variables has a range that is an order of magnitude greater than the range of another variable, the comparative histogram can lose its effectiveness. To illustrate this, consider the following comparative histogram of three widely varying quantities:

The variable in the top histogram has a range that is 10 times the range of the variable in the lower histogram. Consequently, all data in the bottom panel is inside of a single histogram bin. This plot does not reveal anything about the distribution of the third variable.

A different way to compare distributions is to plot a panel of the empirical cumulative distribution functions (CDF). The following call to PROC UNIVARIATE creates these "comparative CDF plots," as well as the comparative histograms, for simulated data:

/* simulate three variables with different distributions */
data Wide;
call streaminit(123);
do i = 1 to 100;
   Uniform = rand("Uniform");
   Normal = rand("Normal");
   Gamma = rand("Gamma", 4);
   output;
end;
run;
 
/* transpose from wide to long data format to create a CLASS variable */
proc transpose data=Wide name=ID out=Long(drop=i);  by i;  run;
data Long; 
   set Long(rename=(Col1=x)); 
   label ID=;
run;
 
/* create comparative histograms and CDF plots */
ods select Histogram CDFPlot;
proc univariate data=Long;
   class ID;
   histogram x / nrows=3;
   cdfplot x / nrows=3;
run;

The comparative CDF plot shows the empirical distributions. The horizontal axis indicates the range of the data: [0, 11] for the Gamma variable, [-3,3] for the Normal variable, and [-0.5, 0.5] for the Uniform variable. However, the plot is far from perfect. Although the CDF plot reveals the approximate center and scale of the data, I would find it difficult to conclude from the comparative CDF plot that the Gamma variable is skewed whereas the Normal variable is symmetric.

The spread plot for distributions

In the classic books Graphical Methods for Data Analysis (Chambers, et al., 1983) and Visualizing Data (Cleveland, 1993), the authors recommend using scatter plots of quantiles to visualize and compare distributions.

Consider rotating the CDF panel 90 degrees in the counter-clockwise direction. If you also center the data so that each variable has zero mean, then the result called a spread plot. The spread plot, as its name implies, is used to compare the spread (range) of distributions, as well as to give some indication of the tail behavior (long tails, outliers, and so forth) in the data.

I have previously blogged about how to compute empirical quantiles. The following SAS/IML program centers each variable and estimates the empirical cumulative probabilities. The SGPANEL procedure is then used to visualize the results:

proc iml;
varNames = {"Gamma" "Normal" "Uniform"};
use Wide;  read all var varNames into x;  close Wide;
 
/* assume nonmissing values */
x = x - mean(x);          /* 1. Center the variables */
N = nrow(x);
F = j(nrow(x), ncol(x));
do i = 1 to ncol(x);      /* 2. Compute quantile of each variable */
   F[,i] = (rank(X[,i])-0.375)/(N+0.25); /* Blom (1958) fractional rank */
end;
ID = repeat(varNames, N); /* 3. ID variables */
create Spread var {ID x F};  append;  close Spread;
quit;
 
title "Comparative Spread Plot";
proc sgpanel data=Spread;
   panelby ID / rows=1 novarname;
   scatter x=f y=x;
   refline 0 / axis=y;
   colaxis label="Proportion of Data";
   rowaxis display=(nolabel);
run;

There does not appear to be a standard term for the quantity on the horizontal axis, which I have labeled "proportion of data." I call it the estimated cumulative probability, but it is also called the estimated proportion, the fractional rank, the estimate for cumulative proportion, or simply plotting positions. In Chambers et al. (1983) it is called the fraction of data. In Cleveland (1993) it is called the f-value, which is short for fractional value. In some SAS output it is labeled "proportion less [than the value shown]." (There is also many different ways to compute the plotting positions. Chambers and Cleveland each use the simpler fractional rank given by r_i–0.5)/(N+1), where r_i is the rank of the ith observations. I have used Blom's 1958 formula, which is used heavily in SAS software.)

You could create a similar plot by using PROC RANK and the DATA step. The nice thing about the spread plot is that the range of the data is evident. The range for the centered Gamma variable is [-3,8]. Furthermore, notice that the mean value is greater than the median value, which is a standard rule of thumb for continuous skewed distributions. In contrast, the Normal and Uniform variables are both visually symmetric. The Normal variable clearly has two moderate tails, whereas the Uniform variable appears to be a bounded distribution.

A plot that uses quantiles to compare distributions is more powerful than the technique of comparing histograms. However, the quantile plot requires more skill to interpret.

Although the spread plot is not a well-known display, it appears prominently in the output of several popular SAS procedures. My next blog post will discuss how to interpret a spread plot when it appears in the output of SAS regression procedures.

Of course, the value 1.2 was completely arbitrary. You could just as well ask for the probability of detecting a difference of size 0.5, 1.0, or 2.0. You can repeat the simulation for any value of the difference and thereby construct a power curve. That's what I do in this article, which is based on material in Chapter 5 of Simulating Data with SAS.

Here is the result of running the previous simulation for differences of size 0, 0.1, 0.2, ..., 2.0. The points along the curve are power estimates from the simulation, along with 95% confidence intervals. The size of the error bars reflects the number of simulations used for each point estimate, which is 5,000 for this simulation. The underlying curve is the result of running PROC POWER, which can solve this problem exactly. For more complicated statistical tests, an exact power curve is not available and simulation is the only way to estimate power.

The simulation takes about 30 seconds to run on a desktop PC from 2010. You can download the complete SAS program that runs the simulation. The program is virtually identical to the program I discussed last week, except that there is now a loop over the parameter that controls the mean difference (do Delta = 0 to 2 by 0.1), and I have included the Delta parameter in two BY statements.

A few comments about the power curve:

• In this simulation, the two populations have unit standard deviations. The graph shows that when the true difference between the means is 0.5, there is a 20% chance of detecting the difference with such a small sample size (n₁=n₂=10). The chance rises to 56% when the difference is 1. To have an 80% chance, the means must differ by about 1.35.
• When there is no difference between the population means, the graph shows that there is a 5% chance that the t test will reject the null hypothesis (which is that the means are equal). This is consistent with the fact that the test is performed at the α = 5% significance level.
• Notice that the confidence intervals for the power are smaller when the curve is flatter and larger when the curve is steeper. You can exploit this to make your simulations run faster by using fewer simulation trials in certain parameter ranges. See the Chapter 6, "Strategies for Efficient and Effective Simulation," of Simulating Data with SAS.

Download the program and experiment with changing the parameters. Change NumSamples to 1000 to make your experiments run faster. What does the curve look like when (n₁=n₂=20)?

There are two ways to pass values from PROC IML to other procedures. One solution stores the values in a macro variable at run time. The other solution calls the SAS procedure from within the SAS/IML program and passes the SAS/IML vector as a parameter.

To give a concrete example, suppose that you have a set of numbers in a SAS/IML vector, and you want to use the numbers as a list of percentile values in PROC UNIVARIATE, as follows:

/* Goal: specify a list of numbers to the PCTLPTS= option */
proc univariate data=sashelp.cars noprint;
   var MPG_City;
   output out=pctl pctlpre=p pctlpts=2.5 10 50 90 97.5; /* numbers go here */
run;

In this example, we assume that the numbers on the PCTLPTS= option are the result of some SAS/IML computation.

The macro variable solution

The first solution converts the vector into a string of values that are separated by blanks. That string can then be placed in a macro variable. The macro variable can be referenced in the SAS procedure. The following steps create a macro variable from values in a SAS/IML row vector:

1. Convert the numeric vector into a character vector. You can use the PUTN function in Base SAS or the SAS/IML CHAR function to convert a numeric vector.
2. Use the ROWCAT function to concatenate the values into a single blank-delimited string of values. (If your values are in a column vector, transpose the vector before calling the ROWCAT function.)
3. Use the SYMPUT function to store the string of values in a macro variable.

You can then call PROC UNIVARIATE and reference the macro variable. The program follows:

proc iml;  
pctl = {2.5 10 50 90 97.5};               /* these are the values  */
s = rowcat( char(pctl) );
call symputx("PctList", s);                     /* create macro variable */
quit;
 
proc univariate data=sashelp.cars noprint;
   var MPG_City;
   output out=pctl pctlpre=p pctlpts=&PctList; /* use SAS/IML values here */
run;

I have used the same trick to pass values from SAS/IML vectors into Base SAS functions that expect a comma-delimited list. Use the macro variable solution when you are finished with PROC IML and are ready to use other procedures.

An easier way: Use the SUBMIT statement

There is even an easier way that passes the values of a SAS/IML vector directly into a SAS procedure without creating a macro variable and without leaving PROC IML: use the SUBMIT and ENDSUBMIT statements.

proc iml;  
pctl = {2.5 10 50 90 97.5};
 
submit  pctl;
   proc univariate data=sashelp.cars noprint;
      var MPG_City;
      output out=pctl pctlpre=p pctlpts=&pctl;  /* list from SAS/IML vector */
   run;
endsubmit;

If you include the name of a SAS/IML vector on the SUBMIT statement, the contents of that vector are substituted for the expression (in this case, &pctl) before the SUBMIT block is sent to SAS. Notice that there is not a macro variable called pctl. Although the &pctl expression looks like a macro variable substitution, it is actually accomplished through text substitution at run time, not through a macro preprocessor.

Use the SUBMIT method when you will need to compute further quantities in PROC IML. To learn more, watch my video about calling SAS functions from the SAS/IML language.

The ability to detect differences in the group means depends on the sample sizes in the study, but it is safe to say that the t test is unlikely to detect small differences in the students' mean scores, and it is more likely to detect larger differences.

In general, the power of a test is the probability that the test will reject the null hypothesis when a specified alternative is true. For the t test, power means the probability that the test can detect a mean difference of a specified magnitude. In general, this can be a difficult computation because it requires knowing the sampling distribution of the test statistic under the alternative hypothesis. For simple tests (such as the two-sample t test), the sampling distribution is known, but for more complicated statistical tests the power computation might be available only by using simulation methods.

This article describes how to use simulation to estimate the power of the t test. For this simple test, we can check the simulation by using the exact answer, as provided by the POWER procedure in SAS. However, the simulation will illustrate general ideas that you can use to estimate the power of more complicated statistical tests. This article is based on material in Chapter 5 of Simulating Data with SAS.

Formulation of the problem

In this simulation, I assume that each population is normally distributed. For simplicity, assume the population for Group 1 is N(0, 1) and the population for Group 2 is N(δ, 1), where δ > 0 is the difference between the population means.

The null hypothesis for the t test is that δ = 0. Given two samples, the t test will either reject the null hypothesis at the α = 0.05 significance level or it won't. In the simulation approach, you simulate many samples, and estimate the probability of rejecting the null hypothesis by using the empirical proportion of simulated samples that rejected the null hypothesis. (This uses the law of large numbers.)

The following steps estimate the power for the two-sample pooled t test:

1. Simulate data from the model for each group's population. These are the samples. The populations are chosen so that the true difference between the population means is δ > 0. (The null hypothesis is false.)
2. Run the TTEST procedure on the samples. For each sample, record whether the t test rejects the null hypothesis.
3. Count how many times the t test rejects the null hypothesis. This proportion is an estimate for the power of the test.

Simulating the data

The following SAS statements define parameters for the simulation and use the DATA step to simulate 5,000 simulated trials. All of the data are in a single data set, and the SampleID variable identifies which observations belong to which simulated trial. In this simulation, each group has 10 observations and the true difference between the population means is 1.2, which is a little larger than the standard deviations of the populations.

%let n1 = 10;             /* group sizes*/
%let n2 = 10;
%let NumSamples = 5000;   /* number of simulated samples */  
%let Delta = 1.2;         /* true size of mean difference in population */
 
data PowerSim(drop=i);
call streaminit(321);
do SampleID = 1 to &NumSamples;
   do i = 1 to &n1;
      c=1;  x = rand("Normal", 0, 1);      /* Group 1: x ~ N(0,1) */
      output;
   end;
   do i = 1 to &n2;
      c=2;  x = rand("Normal", &Delta, 1); /* Group 2: x ~ N(Delta, 1) */
      output;
   end;
end;
run;

Analyzing the Simulated Data

As I have written previously, use BY-group processing to carry out efficient simulation and analysis in SAS. Also, be sure to suppress the display of tables and graphs during the analysis by using the the %ODSoff and macro. The following SAS statements define the %ODSOff and %ODSOn macros, and analyze all data for the simulated trials:

%macro ODSOff(); /* Call prior to BY-group processing */
   ods graphics off;  ods exclude all;  ods noresults;
%mend;
 
%macro ODSOn(); /* Call after BY-group processing */
   ods graphics on;  ods exclude none;  ods results;
%mend;
 
/* 2. Compute (pooled) t test for each sample */
%ODSOff 
proc ttest data=PowerSim; 
   by SampleID; 
   class c; 
   var x; 
   ods output ttests=TTests(where=(method="Pooled")); 
run; 
%ODSOn

The TTEST procedure creates an output data set that contains 5,000 rows, one for each simulated trial. The data set includes a variable named Probt, which gives the result of the t test on each trial.

Count the rejections to estimate power

The last step in the simulation is to estimate the power, which is the probability of rejecting the null hypothesis. The following SAS statements create an indicator variable that has the value 1 if the t test rejected the null hypothesis at the 0.05 significance level, and the value 0 otherwise. (Alternatively, you could define a SAS format.) You can then use PROC FREQ to count the proportion of trials for which the t test rejected the null hypothesis.

/* Construct indicator var for obs that reject H0 at 0.05 significance */ 
data Results; 
   set TTests; 
   RejectH0 = (Probt <= 0.05); 
run; 
 
/* 3. Compute proportion: (# that reject H0)/NumSamples and CI */ 
proc freq data=Results; 
   tables RejectH0 / nocum binomial(level='1');
run;

The FREQ procedure indicates that the power of the two-sample t test is about 72%. The 95% confidence interval for that estimate is [0.708, 0.733]. This estimate is for the scenario of samples of sizes 10, where one sample is drawn from N(0,1) and the other is drawn from N(1.2, 1).

As mentioned earlier, you can use PROC POWER to find the exact power for the two-sample t test, as follows:

proc power;
  twosamplemeans  power = .           /* missing ==> "compute this" */
    meandiff=1.2 stddev=1 ntotal=20;  /* 20 obs in the two samples  */
run;

The estimate from the simulation is very close to the true power of the t test.

Next week I will extend this simulation to estimate points along a power curve.

Here is how the fit plot looks for fewer than 5,000 observations:

/* simulate data from a linear regression model */
data RegData;
call streaminit(1);
do i = 1 to 6000;
   x = rand("Uniform");
   y = 1 + 2*x + 3*rand("Normal");
   if i < 20 then y = y + 8*rand("Normal"); /* add a few outliers */
   output;
end;
 
ods select FitPlot(persist);
proc reg data=RegData(obs=4000) plots(only)=FitPlot;
   model y = x;
quit;

With fewer than 5,000 observations, I get the usual fit plot that consists of a scatter plot overlaid with a curve of predicted value, a band for the confidence interval for the mean, and dashed lines that indicate the confidence intervals for individual predictions. (The confidence limits are barely visible. Click the graph to enlarge it.) However, watch what happens when I use more than 5,000 observations:

proc reg data=RegData plots(only)=FitPlot;
   model y = x;
quit;
ods select all;

The scatter plot is gone, replaced by a heat map that shows the density of the data. The predicted values are still present, although the graphical style used to draw it is different, which results in a red line. The confidence intervals are gone.

Overall, this is a nice feature and I think that the change is a good idea. The reason for the change is easy to understand: scatter plots suffer from overplotting when there are many points, so it is more useful to visualize the density of the observations than the individual observations. Furthermore, although scatter plots are very fast to construct, when there are many points a heat map (which bins observations) is faster to compute and render than a scatter plot.

The plot does not currently include confidence intervals, but there is no reason why these can't be added in a future release. However, the confidence interval for the mean predictions will usually be tiny for large data sets—already it is barely visible in the plot of 4,000 points.

Controlling the appearance of the heat map

Prior to SAS/STAT 12.1, the REG procedure created a fit plot as a scatter plot for small data sets (less than 5,000 points). For larger sample sizes, the procedure suppressed the fit plot. The behavior was controlled by using MAXPOINTS= option on the PLOTS= option on the PROC REG statement.

In SAS/STAT 12.1, the MAXPOINTS= option accepts two arguments, and the default values are MAXPOINTS=5000 150000. The first argument specifies the data size for which heat maps are used instead of scatter plots. The documentation of the MAXPOINTS= option states that "when the number of points exceeds [the first number] but does not exceed [the second number] divided by the number of independent variables, heat maps are displayed instead of scatter plots for the fit and residual plots." In other words, if you have a regression with k explanatory variables, heat maps are used when the number of observations is between 5,000 and 150,000/k. Of course, you can use the MAXPOINTS= option to change either or both of those values.

Any comments on this new behavior in PROC REG?

About 50 SAS/STAT procedures support the NOPRINT option in the PROC statement. When you specify the NOPRINT option, ODS is temporarily disabled while the procedure runs. This prevents SAS from displaying tables and graphs that would otherwise be produced for each BY group. For a simulation that computes statistics for thousands of BY groups, suppressing the display of tables results in a substantial savings of time.

Newer SAS procedures do not always support a NOPRINT statement. However, you can still suppress the ODS output. The following macros encapsulate statements that turn the ODS system off and on. I call the %ODSOff macro before I start the BY-group analysis; I call the %ODSOn macro after the analysis completes.

%macro ODSOff(); /* Call prior to BY-group processing */
ods graphics off;
ods exclude all;
ods noresults;
%mend;
 
%macro ODSOn(); /* Call after BY-group processing */
ods graphics on;
ods exclude none;
ods results;
%mend;

For example, if I were using PROC ROBUSTREG to analyze many samples of simulated data, I might use the following pseudo-code:

%ODSOff
proc robustreg data=MySimData;
   BY SampleID;
   model y = x;
   ods output ParameterEstimates = OutputStats;  /* <== insert name of ODS table */
run;
%ODSOn

Even though ODS is suppressed to the display destinations (such as LISTING and HTML), you can capture the statistics that result from each analysis by using an ODS OUTPUT statement, which saves an ODS table to a SAS data set. Other ways to save statistics include using an OUTPUT statement, an OUT= or OUTEST= data set, and so forth.

Be aware that some SAS procedures (such as PROC MIXED) write a NOTE to the SAS log as part of their normal operation. The NOTE might say something like "NOTE: Convergence criteria met." For these procedures, you will also want to turn off notes, lest they fill the SAS log:

%ODSOff
options nonotes;  /* use NONOTES to suppress notes to the log */
proc mixed ...;
model y = ...;
run;
options notes;   /* turn NOTES back on */
%ODSOn

The material in this blog post is taken from my book Simulating Data with SAS, which contains many more tips and techniques for the efficient simulation of data.

Someone asked me what kind of speed-ups can be expected by vectorizing a SAS/IML program. In general, the answer depends on the size of your data and the particular operations that you are performing, but it is straightforward to run a series of examples that compare the performance at each step of the vectorization of the previous example.

In this post, I generate a square matrix of size N for N=100, 200, 400, 600, 800, and 1000. For each size, I time how long it takes to compute an operation that is equivalent to X`*X. The operations are as follows:

• Level-0: The original program, which consists of three nested loops and scalar operations.
• Level-1: A program that contains two nested loops and dot products of vectors.
• Level-2: A program that contains one loops and matrix-vector multiplication.
• Level-3: A program that contains one matrix-matrix multiplication.

You can download the program that performs the experiment.

On my desktop computer, the results are summarized by the following plots. The first plot shows all of the times. However, the graph is dominated by the slow times that are associated with the original program that had three nested loops and only scalar operations. The middle plot omits the times for the original program. In this view you can see that the level-1 operation is about 25 times faster than the level-0 operations for N=1000. Again, however, the graph is dominated by the slowest method, and so the last plot shows only the timings for the level-2 and level-3 operations. The level-2 operations are about 16 times faster than the level-3 operations for N=1000. The level-3 operations are almost three times faster than level-2.

In summary, the slowest method (scalar operations in three nested loops) is about 1,000 times slower than the equivalent level-3 operation, which does not require any loops. This is good motivation: the time you invest to fully vectorize your code can pay dividends by running 1,000 times faster.

If you are adept at interpreting logarithmic scales, you can also see the results plotted on a Log10 axis.

No matter how you visualize the results, the conclusion is the same: if you want your program to scale to handle large data, vectorize the program to take advantage of the SAS/IML matrix operations.