Suppose you are interested in statistics related to movies released in the US during the time period 2005-2007. The Movies data set contains a sampling of movies that were also rated by the kids-in-mind.com Web site. Strictly speaking, these movies are not a random sample from the population of all US movies, since the data includes all major movies with nationwide release. (In fact, the data essentially equal the population of US movies; the data set is lacking only minor releases that had limited distribution.) In discussing bootstrap statistics, it is important to differentiate between a random sample and the underlying population. To make this distinction perfectly clear, the following DATA step extracts a random sample of approximately 25% of the movies into the Sample data set:

/* extract random sample of movies for 2005-2007 */
submit;
data Sample;
   set Sasuser.Movies;
   if (ranuni(1)<0.25);
run;
endsubmit;

Figure 14.1 shows the distribution of the budgets for movies in the Sample data set. The histogram is created by the following statements:

/* read in movies and draw histogram */
use Sample;
read all var {Budget} into x;
close Sample;

declare Histogram hist;
hist = Histogram.Create("Movies", x);
hist.SetAxisLabel(XAXIS,"Sample Movie Budgets");

Figure 14.1. Distribution of Movie Budgets in Sample

The data are not normally distributed. The sample mean of the budgets is $39.2 million. There are several movies with very large budgets, although the complete Movies data set has three movies with budgets greater than $200 million that are not included in the Sample data. (For comparison, Figure 7.4 shows the distribution for all of the budgets in the Movies data set.) The remainder of this chapter analyzes the Sample data.

Implementing a bootstrap estimate often includes the following steps:

For the movies in the Sample data set, you can use bootstrap methods to examine the sampling distribution of the mean (or any other statistic). Since the first step in the bootstrap method is to resample (with replacement) from the data, it is convenient to use the SampleWithReplace module that was developed in Chapter 13, "Sampling and Simulation."

The following statements use the SampleWithReplace module to carry out a simple bootstrap on the budgets for movies:

/* bootstrap method: the statistic to bootstrap is the mean */
/* 1. Compute the statistic on the original data */
Mean = x[:];                             /* sample mean             */

/* 2. Resample B times from the data (with replacement)
   to form B bootstrap samples. */
call randseed(12345);
load module=SampleWithReplace;           /* not required in IMLPlus */
B = 1000;
n = nrow(x);
EQUAL = .;
xBoot = SampleWithReplace(x, B||n, EQUAL);

/* 3. compute the statistic on each resample */
s = xBoot[,:];

As indicated previously, the main steps of the bootstrap algorithm are as follows:

The vector s in the previous program contains B statistics. The union of the statistics is the bootstrap distribution for the statistic. It is an estimate for the sampling distribution of the statistic. A histogram of s is shown in Figure 14.2 and is created with the following IMLPlus statements:

/* 4. If desired, graph the bootstrap distribution */
declare Histogram hBoot;
hBoot = Histogram.Create("Means", s);
hBoot.SetAxisLabel(XAXIS,"Means of Resamples");
hBoot.ReBin(40, 1);
attrib = BLACK || SOLID ||2;
run abline(hBoot, Mean, ., attrib);

Figure 14.2. Estimate of the Sampling Distribution of the Mean

The figure shows that most of the means for the resampled data are between $35 and $43 million per movie. A small number were as small as $27 million per movie, or as large as $54 million per movie. The sample statistic has the value 39.2, which is shown as a vertical line in Figure 14.2. This figure enables you to see the variation in the mean of the resampled data.

This section shows how you can use the bootstrap distribution shown in Figure 14.2 to compute bootstrap estimates for familiar quantities: the mean, the standard error of the mean (SEM), and the confidence intervals for the mean. The mean of the bootstrap distribution is an estimate of the mean of the original sample. The standard deviation of the bootstrap distribution is the bootstrap estimate of the standard error of s. Percentiles of the bootstrap distribution are the simplest estimates for confidence intervals. The following statements (which continue the program in the previous section) compute these quantitites:

/* 5. Compute standard errors and confidence intervals. */
MeanBoot = s[:];                   /* a. mean of bootstrap dist     */
StdErrBoot = sqrt(var(s));         /* b. std error                  */
alpha = 0.05;
prob = alpha/2 || 1-alpha/2;       /* lower/upper percentiles       */
call qntl(CIBoot, s, prob);        /* c. quantiles of sampling dist */
print MeanBoot StdErrBoot CIBoot;

Figure 14.3. Bootstrap Estimates for the Sampling Distribution of the Mean

The output is shown in Figure 14.3. The following list describes the main steps:

You can visualize these quantities by adding vertical lines to the histogram of the sampling distribution. The following statements draw lines on Figure 14.2. The resulting graph is shown in Figure 14.4.

attrib = BLACK || DASHED || 1 || PLOTFOREGROUND;
run abline(hBoot, MeanBoot//CIBoot, {.,.,.}, attrib);

Figure 14.4. Estimate of the Mean and Confidence Intervals for the Sampling Distribution

In Figure 14.4, the solid line is the original statistic on the data. The dashed lines are bootstrap estimates.

The power of the bootstrap is that you can create bootstrap estimates for any statistic, s, even those for which there are no closed-form expression for the standard error and confidence intervals. However, because the statistic in this example is the mean, you can use elementary statistical theory to compute classical statistics (based on the sample) and compare the bootstrap estimates to the classical estimates.

The sampling distribution of the mean is well understood. The central limit theorem says that the sampling distribution of the mean is approximately normally distributed with parameters (μ,σ). The sample mean, , is an estimate for μ. An estimate for σ is the quantity where s²_x is the sample variance and n is the size of the sample. Furthermore, the two-sided 100(1 − α)% confidence interval for the mean has upper and lower limits given by where t_{1−α/2,n−1} is the 1 − α/2 percentile of the Student's t distribution with n − 1 degrees of freedom. The following program computes these estimates:

/* Compute traditional estimates for the sampling distribution of
   the mean by computing statistics of the original data */
StdErr = sqrt(var(x)/n);           /* estimate SEM                  */
t = quantile("T", prob, n-1);      /* percentiles of t distribution */
CI = Mean + t * StdErr;            /* 95% confidence interval       */
print Mean StdErr CI;

Figure 14.5. Estimates for the Sampling Distribution of the Mean

The output is shown in Figure 14.5. The sample mean is 39.15. The SEM is estimated by 3.76 and the 95% confidence interval by [31.7,46.6].

A comparison of Figure 14.5 with Figure 14.3 shows that the classical estimates are similar to the bootstrap estimates. However, note that the classical confidence interval is symmetric about the mean (by definition), whereas the bootstrap confidence interval is usually not symmetric. The bootstrap estimate is based on percentiles of the statistic on the resampled data, not on any formula. This fact implies that bootstrap confidence intervals might be useful for computing confidence intervals for a statistic whose sampling distribution is not symmetric.

It is interesting to note that the mean budget in the Movies data set is $44.8 million, which falls within the 95% confidence limits for either method. The Movies data set contains essentially all (major) movies released in 2005-2007, so you can take $44.8 million to be the mean of the population.

Some SAS programmers use the DATA step for everything, including generating random samples from data. The DATA step can sample reasonably fast. But it is even faster to use the SURVEY-SELECT procedure for generating a random sample.

The SURVEYSELECT procedure can generate B random samples of the input data, and is therefore ideal for bootstrap resampling (Cassell 2007, 2010; Wicklin 2008). It implements many different sampling techniques, but this section focuses on random sampling with replacement. The SURVEYSELECT procedure refers to "sampling with replacement" as "unrestricted random sampling" (URS).

As an example, suppose your data consists of four values: A, B, C, and D. The following statements create a data set named BootIn that contains the data:

/* create small data set */
x = {A,B,C,D};
create BootIn var {"x"};
append;
close BootIn;

Suppose you want to generate B bootstrap resamples from the data in the BootIn data set. The following statements use the SUBMIT statement to call the SURVEYSELECT procedure and to pass in two parameters to the procedure. The procedure reads data from BootIn and writes the resamples into the BootSamp data set.

/* Use the SURVEYSELECT procedure to generate bootstrap resamples */
N = nrow(x);
B = 5;

submit N B;
proc surveyselect data=BootIn out=BootSamp noprint
     seed   = 12345               /* 1 */
     method = urs                 /* 2 */
     n      = &N                  /* 3 */
     rep    = &B                  /* 4 */
     OUTHITS;                     /* 5 */
run;
endsubmit;

The following list describes the options to the SURVEYSELECT procedure:

The output data set created by the SURVEYSELECT procedure is displayed by the following statements, and is shown in Figure 14.8:

submit;
proc print data=BootSamp;
   var Replicate x;
run;
endsubmit;

Figure 14.8. Bootstrap Resamples Using OUTHITS Option

It is not necessary to specify the OUTHITS option. In fact, it is more efficient to omit the OUTHITS option. If you omit the OUTHITS option, the NumberHits variable in the output data set contains the frequency that each observation is selected. Omitting the OUTHITS option and using NumberHits as a frequency variable when computing statistics is typically faster than using the OUTHITS option because there are fewer observations in the BootSamp data set to read. Most, although not all, SAS procedures have a FREQ statement for specifying a frequency variable. Figure 14.9 shows the output data set created by the SURVEYSELECT procedure when you omit the OUTHITS option:

Figure 14.9. Bootstrap Resamples without Using OUTHITS Option

The BootSamp data set contains a copy of all variables in the input data set and contains the variable Replicate, which identifies the bootstrap resample that is associated with each observation. The data set also contains the NumberHits variable, which can be used on a FREQ statement to specify the frequency of each observation in the resample.

The preceding section decribed how to use the SURVEYSELECT procedure to generate B bootstrap resamples. This section shows how to use the BY statement (and, if desired, the FREQ statement) of a SAS procedure to generate the bootstrap distribution of a statistic.

This section duplicates the example in Section 14.2.2 by creating a bootstrap distribution for the mean statistic. However, instead of computing the means of each resample in the SAS/IML language, this section uses the MEANS procedure to do the same analysis.

Assume that the data set Sample contains the Budget variable for a subset of movies, as described in Section 14.2.2. The following program uses the SURVEYSELECT procedure to generate bootstrap resamples and uses the MEANS procedure to generate the bootstrap distribution for the mean statistic:

/* use SURVEYSELECT to resample; use PROC to compute statistics */
DSName = "Sample";                          /* 1 */
VarName = "Budget";
B = 1000;

submit DSName VarName B;                    /* 2 */
/* generate bootstrap resamples */
proc surveyselect data=&DSName out=BootSamp noprint
     seed=12345 method=urs rep= &B
     rate = 1;                              /* 3 */
run;

/* use procedure to compute statistic on each resample */
proc means data=BootSamp noprint;           /* 4 */
   by Replicate;                            /* a */
   freq NumberHits;                         /* b */
   var &VarName;
   output out=BootDist mean=s;              /* c */
run;
endsubmit;

The following list describes the main steps of the program:

You can visualize the bootstrap distribution by using IMLPlus graphics, as shown in the following statements and in Figure 14.10:

/* create histogram of bootstrap distribution */
use BootDist;
read all var {s};
close BootDist;

declare Histogram hBoot;
hBoot = Histogram.Create("Means", s);
hBoot.SetAxisLabel(XAXIS,"Means of Resamples (PROC MEANS)");
hBoot.ReBin(40, 1);

Figure 14.10. Estimate of the Sampling Distribution of the Mean

You can compare Figure 14.10, which was computed by calling the SURVEYSELECT and MEANS procedures, with Figure 14.4, which was computed entirely in the SAS/IML language. The distributions are essentially the same. The small differences are due to the fact that the random numbers that were used to generate the resamples in SAS/IML are different from the random numbers used to generate the resamples in the SURVEYSELECT procedure.

You can create a scatter plot that consists of all of the B scree plots generated by the principal component analysis computed on the B resamples:

/* create union of scree plots for bootstrap resamples */
declare DataObject dobj;
dobj = DataObject.CreateFromServerDataSet("Work.BootEVals");
dobj.SetVarFormat("Proportion","BEST5.");

declare ScatterPlot p;
p = ScatterPlot.Create(dobj, "Number", "Proportion");
p.SetMarkerSize(3);

To that scatter plot, you can overlay the PVE for the original sample, in addition to confidence intervals for the underlying parameter in the population of movies. The following statements overlay the quantities; the result is shown in Figure 14.14.

/* overlay statistics of the original data */
use EigenValues;                 /* 1. read PVE for data */
read all var {Number Proportion};
close EigenValues;
p.DrawUseDataCoordinates();      /* define coordinate system         */
p.DrawMarker(Number, Proportion, MARKER_X, 7); /* plot proportions   */

use BootEVals;                   /* 2. read stats for resamples      */
read all var {Number Proportion};
close BootEVals;

NumPC = ncol(VarNames);          /* number of principal components   */
s = shape(Proportion, 0, NumPC); /* 3. reshape results               */

mean = s[:,];                    /* 4. compute mean of each column   */
alpha = 0.05;                    /* significance level for C.I.      */
prob = alpha/2 || 1-alpha/2;     /* lower/upper values for quantiles */
call qntl(CI, s, prob);          /* compute C.I. for each column     */
print mean, CI[rowname={LCI UCI} label="Confidence Intervals"];

p.DrawSetRegion(PLOTBACKGROUND);
p.DrawSetBrushColor(0xC8C8C8);   /* light gray */
dx = 0.1;                        /* half-width of rectangles         */
do i = 1 to NumPC;               /* 5. draw rectangles and mean line */
   p.DrawRectangle(i-dx, CI[1,i], i+dx, CI[2,i], true);
   p.DrawLine(i-dx, mean[i], i+dx, mean[i]);
end;

The program consists of the following main steps:

Figure 14.13. Principal Component Analysis

Figure 14.14. Scree Plot with Bootstrap Confidence Intervals

Notice that Figure 14.13 answers the question "how much uncertainty is in the PVE estimate for the sample data?" The answer is that there is quite a bit of uncertainty. The estimate for the proportion of variance explained by the first principal component is 0.416, but the confidence interval is fairly wide: [0.368,0.484]. The half-width of the confidence interval is more than 10% of the size of the estimate. Similar results apply to proportion of variance explained by the other principal components.

If you want to view the bootstrap distributions of the PVE for each principal component, you can create a panel of histograms. The i th column of the s matrix (computed in the previous section) contains the data for the ith histogram. The mean and ci matrices contain the mean and 95% confidence intervals for each PVE. The following statements create a data object from the s matrix. For each column, the program creates a histogram and overlays the mean and bootstrap percentile confidence interval. The resulting panel of histograms is shown in Figure 14.15.

/* create data object from bootstrap distribution;
   create histograms and overlay bootstrap mean and C.I. */
names = "Proportion1":"Proportion"+strip(char(NumPC));
declare DataObject dobj2;
dobj2 = DataObject.Create("PVE", names, s);

declare Histogram hist;
attrib = BLACK || DASHED || 1 || PLOTFOREGROUND;
do i = 1 to NumPC;
   hist = Histogram.Create(dobj2, names[i]);
   /* use a module distributed with SAS/IML Studio to compute
      the position of i_th plot in an array of 6 plots */
   run CalcPlotPosition(i, 6, left, top, width, height);
   hist.SetWindowPosition(left, top, width, height);

   /* plot mean and C.I. */
   run abline(hist, mean[i]//CI[,i], {.,.,.}, attrib);

   /* draw rectangle in background to indicate C.I. */
   hist.DrawUseDataCoordinates();
   hist.GetAxisViewRange(YAXIS, YMin, YMax);
   hist.DrawSetRegion(PLOTBACKGROUND);
   hist.DrawSetBrushColor(0xC8C8C8);
   hist.DrawRectangle(CI[1,i], 0, CI[2,i], YMax+10, true);
end;

Figure 14.15 shows the bootstrap distributions for the proportion of variance explained by each principal component. The distributions are slightly skewed; this is especially noticeable for Propor-tion2.

The figure also shows how the dynamically linked graphics in SAS/IML Studio can give additional insight into the fact that these five statistics are not independent, since the sum of the five proportions must equal unity. In the figure, large values of the Proportion1 variable are selected. The selected observations correspond to bootstrap resamples for which the value of the Proportion1 variable is large. For these bootstrap resamples, the Proportion2 variable tends to be smaller than the mean taken over all bootstrap resamples. Similar relationships hold for the other variables, and also for the bootstrap resamples for which the Proportion1 variable is small.

Figure 14.15. Panel of Histograms of Bootstrap Distributions