SAS博客列表

Many modern statistical techniques incorporate randomness: simulation, bootstrapping, random forests, and so forth. To use the technique, you need to specify a seed value, which determines pseudorandom numbers that are used in the algorithm. Consequently, the seed value also determines the results of the algorithm. In theory, if you know the seed value and the internal details of the pseudorandom algorithm, then the stream is completely determined, and the results of an algorithm are reproducible. For example, if I publish code for a simulation or bootstrap method in SAS, you can reproduce my computations as long as my program specifies the seed value for every part of the program that uses random numbers.

Occasionally, programmers post questions on discussion forums about how to compare the results that were obtained by using different software. A typical question might be, "My colleague used R to run a simulation. He gets a confidence interval of [14.585, 21.42]. When I run the simulation in SAS, I get [14.52, 21.19]. How can I get the same results as my colleague?" In most cases, you can't get the same answer. More importantly, it doesn't matter: these two results are equally valid estimates. When you run an algorithm that uses pseudorandom numbers, you don't get THE answer, you get AN answer. There are many possible answers you can get. Each time you change the random number stream, you will get a different answer, but all answers are equally valid.

This article examines how the random number seed affects the results of a bootstrap analysis by using an example of a small data set that has N=10 observations. I will run the same bootstrap analysis for 50 different seeds and compare the results of a confidence interval. The same ideas apply to simulation studies, machine learning algorithms, and other algorithms that use randomness to obtain an answer.

A bootstrap estimate of a confidence interval

The following data were used by Davison, Hinkley, and Schechtman (1986). They ran an analysis to examine the bootstrap distribution of the sample mean in a small nonnormal set of data. The goal of the following analysis is to obtain a bootstrap estimate of a 90% confidence interval for the population mean. Because we are going to run the same analysis many times with different random number seeds, I will put the bootstrap steps in a SAS macro. The following statements run two bootstrap analyses on the same data. The only difference is the value of the random-number seeds.

data Sample;
input x @@;
datalines;
9.6 10.4 13.0 15.0 16.6 17.2 17.3 21.8 24.0 33.8
;
 
/* macro to perform one bootstrap analysis that uses NumSamples resamples
   and uses a particular seed. Results are in the data set BootEst */
%macro Do1Boot(seed, NumSamples);
   /* 1. generate NumSamples bootstrap resamples from SEED value */
   proc surveyselect data=Sample NOPRINT seed=&SEED
        method=urs              /* URS = resample with replacement */
        samprate=1              /* each bootstrap sample has N obs */
        OUTHITS                 /* do not use a frequency var */
        reps=&NumSamples(repname=SampleID)        /* generate NumSamples samples */
        out=BootSamp;
   run;
 
   /* 2. compute sample mean for each resample; write to data set */
   proc means data=BootSamp NOPRINT;
      by SampleID;
      var x;
      output out=BootStats mean=Mean;  /* bootstrap distribution */
   run;
   /* 3. find boostrap estimates of mean, stdErr, and 90% CI */
   proc means data=BootStats NOPRINT;
      var Mean;
      output n=N mean=Mean stddev=StdErr p5=LCL90 p95=UCL90
             out=BootEst(drop=_TYPE_ _FREQ_);
   run;
   /* add the seed value to these statistics */
   data BootEst;
      Seed = &seed;
      set BootEst;
   run;
%mend;
 
title "A Boostrap Analysis for Seed=1";
%Do1Boot(1, 1000);
proc print data=BootEst noobs; run;
 
title "A Boostrap Analysis for Seed=1234";
%Do1Boot(1234, 1000);
proc print data=BootEst noobs; run;

The output shows that the bootstrap estimates are similar, but not exactly the same. This is expected. Each analysis uses a different seed, which means that the bootstrap resamples are different between the two runs. Accordingly, the means of the resamples are different, and so the bootstrap distributions are different. Consequently, the descriptive statistics that summarize the bootstrap distributions will be different. To make the differences easier to see, I used only 1,000 samples in each bootstrap analysis. If you increase the number of resamples to, say, 10,000, the differences between the two analyses will be smaller.

The fact that the random number seed "changes" the result makes some people nervous. Which value should you report to a journal, to a customer, or to your boss? Which is "the right answer"?

Both answers are correct. A bootstrap analysis relies on randomness. The results from the analysis depend on the pseudorandom numbers that are used.

Do the results depend on the seed value?

You should think of each bootstrap estimate as being accurate to some number of digits, where the number of digits depends on the number of resamples, B. For some estimates, such as the bootstrap standard error, you can calculate that the accuracy is proportional to 1/sqrt(B). Thus, larger values for B tend to shrink the difference between two bootstrap analyses that use different random number streams.

One way to visualize the randomness in a bootstrap estimate is to run many bootstrap analyses that are identical except for the random-number seed. The following SAS macro runs 50 bootstrap analyses and uses the seed values 1,2,3,..., 50:

/* Run nTrials bootstrap analyses. The i_th trials uses SEED=i. The estimates are written 
   to the data set AllBoot by concatenating the results from each analysis. */
%macro DoNBoot(NumSamples, nTrials);
   options nonotes;
   proc datasets noprint nowarn;   /* delete data set */
      delete AllBoot;
   quit;
   %do seed = 1 %to &nTrials;
       %Do1Boot(&seed, &numSamples);   /* run with seed=1,2,3,...,nTrials */
       proc append base=AllBoot data=BootEst force; 
       run;
   %end;
   options notes;
%mend;
 
/* Run 50 bootstrap analyses: numSamples=1000, nTrials=50, seeds=1,2,3,...,50 */
%DoNBoot(1000, 50); 
 
/* find boostrap estimates of mean, stdErr, and 90% CI */
title "Variation in 50 Traditional Boostrap Estimates";
proc means data=AllBoot N Mean StdDev Min Max;
   var Mean StdErr LCL90 UCL90;
run;

The table shows how four estimates vary across the 50 bootstrap analyses:

• The bootstrap estimate of the mean varies in the interval [17.76, 18.00]. The average of the 50 estimates is 17.88.
• The bootstrap estimate of the standard error varies in the interval [2.06, 2.25]. The average of the 50 estimates is 17.88.
• The bootstrap estimate of the lower 90% confidence interval varies in the interval [14.41, 14.80]. The average of the 50 estimates is 14.56.
• The bootstrap estimate of the upper 90% confidence interval varies in the interval [21.36, 22.09]. The average of the 50 estimates is 21.63.

The table shows that the variation among the 50 estimates is small (and can be made smaller by increasing B). If you and your colleague run an identical bootstrap analysis but use different seed values, you are likely to obtain results that are similar to each other. Even more comforting is that the inferences that you make from the analyses should be the same. For example, suppose you are testing the null hypotheses that these data came from a population in which the mean is μ=20. Let's look at the 90% confidence intervals for the 50 analyses:

title "Traditional Bootstrap Estimates";
title2 "Varying Seed Values";
proc sgplot data=AllBoot;
   refline &ObsStat / axis=Y label="Sample Mean";
   highlow x=Seed low=LCL90 high=UCL90 / legendlabel="90% CL";
   scatter x=Seed y=Mean / legendlabel="Bootstrap Estimate";
   yaxis label="x" ;
run;

The graph shows that the value μ=20 is inside all 50 of the confidence intervals. Therefore, the decision to accept or reject the null hypothesis does not depend on the seed value when μ=20.

On the other hand, if you are testing whether the population mean is μ=21.5, then some analyses would reject the null hypothesis whereas others would fail to reject it. What to do in that situation? I would acknowledge that there is uncertainty in the inference and conclude that it cannot be reliably decided based on these data and these analyses. It would be easy to generate additional bootstrap samples (increase B). You could also advocate for collecting additional data, if possible.

Summary

In summary, the result of a bootstrap analysis (and ANY statistical technique that uses randomness) depends on the stream of pseudorandom numbers that are used. If you change the seed or use different statistical software, you are likely to obtain a different result, but the two results should be similar in many cases. You can repeat the analysis several times with different seeds to understand how the results vary according to the random-number seed.

If some seed values enable you to reject a hypothesis and others do not, you should acknowledge that fact. It is not ethical to use the seed value that gives the result that you want to see! (See "How to lie with a simulation" for an example.)

If you are nervous about the result of an algorithm that uses randomness, I suggest that you rerun the analysis several times with different seed values. If the results are consistent across seeds, then you can have confidence in the analysis. If you get vastly different results, then you are justified in being nervous.

The post How much does a bootstrap estimate depend on the random number stream? appeared first on The DO Loop.

A reader asked whether it is possible to find a bootstrap sample that has some desirable properties. I am using the term "bootstrap sample" to refer to the result of randomly resampling with replacement from a data set. Specifically, he wanted to find a bootstrap sample that has a specific value of the mean and standard deviation. If not, can we find a bootstrap sample whose mean and standard deviation are close to the target value? It's an interesting question. Each bootstrap sample is composed of the original observations, but because it is a sample with replacement, some observations can appear more often than others.

For very small values of N, you can solve the reader's problem by generating all bootstrap samples. In a previous article, I used the data {-1, -0.1, -0.6, 0, 0.5, 0.8, 4}, which has the sample mean 0.51 for N=7. I mentioned that one possible bootstrap sample is {4, 4, 4, 4, 0.5, 4, 4}, which has the mean 3.5. This bootstrap sample contains only large values in the original data. The mean of that sample is much greater than the mean of the original data, and the standard deviation is much less.

Of course, generating all bootstrap samples is impractical for even moderate values of N. Furthermore, bootstrap samples like this are improbable if the probability of choosing each observation is 1/N. But what if you assign a higher probability to select some observations and a lesser probability to others? For example, the sample that has six 4s is not improbable if the probability of choosing the 4 is 0.9. If you sample with unequal probability (and still with replacement) the statistics in the samples might be very different from the statistics in the original data.

This kind of sampling is no longer classical bootstrap sampling. It could be considered a kind of importance sampling. Whatever you call it, each sample you generate is a resampling from the original data, but you can control the kinds of samples that you generate by adjusting the probability of selecting each observation.

Sampling with non-uniform probability

Recall that there are four main ways to sample. You can sample with or without replacement; you can sample with equal or unequal probability. In this article, we use sampling with replacement with unequal probability. SAS supports several ways to sample with replacement and with unequal probability. This article uses PROC SURVEYSELECT with the METHOD=PPS_WR option. The acronym PPS_WR stands for "probability proportional to size (PPS), with replacement."

To illustrate sampling with unequal probability, consider the following data with N=1000 that are simulated from a uniform distribution on [0, 20]:

%let N = 1000;    /* sample size */
data Have;
call streaminit(54321);
do i = 1 to &N;
   x = round( rand("uniform", 0, 20), 0.01);  /* round to nearest 0.01 */
   output;
end;
run;
 
/* compute samples mean and std dev of original data */
proc means data=Have mean std min max;
  var x;
run;

The mean and standard deviation of the data are 10.04 and 5.80, respectively. If you were to perform the usual bootstrap sampling, which selects each observation with equal probability, then most samples would have a mean in the interval [9.5, 10.5] and a standard deviation in the interval [5.6, 6.0].

Suppose you are seeking a resample that has unusual properties, such as Mean=9 and StdDev=3. It would be highly unusual to obtain these statistics by using the usual method of resampling observations with equal probability. However, you can assign a high probability to observations whose values are near the value 9. The following DATA step uses the normal density function for the N(9,3) distribution to assign the probability of selecting an observation:

/* Target: mean=9; SD = 3 */
%let mean = 9;
%let StdDev = 3;
/* Use the normal distribution to assign the probability of choosing an observation */
data Have2;
set Have;
Prob = pdf("Normal", x, &mean, &StdDev);
run;

You can use the METHOD=PPS_WR option to resample from the data by using this probability of selecting each observation. The following call to PROC SURVEYSELECT generates a resample of size N from the data but gives preference to observations that are closer to x=9.

/* Sample with replacement from the data, where the probability 
   of choosing an obs is given by the lProb variable */
proc surveyselect noprint data=Have2 out=Resample
     method=PPS_WR 
     seed=12345 N=&N;       /* use OUTHITS option if you want N obs */
   size Prob;               /* specify the probability variable */
run;
 
proc means data=Resample mean std;  /* a different seed will give a different result */
  var x;
  freq NumberHits;
run;

The output from PROC MEANS shows that for these uniformly distributed data, you can obtain a resample whose mean and standard deviations are close to the target values. Of course, using a different random number seed will result in a different resample and different statistics.

You can create a graph to visualize the observations that were chosen for this resample. The following graph plots the NumberHits variable, which records the number of times that each observation appears in the sample:

title "Number of Times Each Observation Appears in Resample";
proc sgplot data=Resample;
   scatter x=x y=NumberHits / markerattrs=(symbol=CircleFilled) transparency=0.75;
   xaxis values=(0 to 20 by 2);
   yaxis values=(1 to 8) valueshint;
run;

You can see that no observations for which x ≥ 18 were selected. Observations for which x ≈ 9 are selected more often than observations for which x is far from 9.

Random samples whose statistics are close to target values

If you change the random number seed, you will get a different resample that has different statistics. If you repeat this process many times, you can choose the resample whose statistics are closest to the target values. You can add the REPS= option to PROC SURVEYSELECT to generate additional resamples. This creates a variable (Replicate) that you can use as a BY-group variable in PROC MEANS. For each (Mean, StdDev) pair, you can compute a distance to the target value. You can then sort the statistics and identify the resamples whose statistics are closest to the target, as follows:

proc surveyselect noprint data=Have2 out=Resample100
     method=PPS_WR 
     Reps=100
     seed=12345 N=&N;       /* use OUTHITS option if you want N obs */
   size Prob;               /* specify the probability variable */
run;
 
/* Compute mean and StdDev of each resample */
proc means data=Resample100 mean std noprint;
  by Replicate;
  var x;
  freq NumberHits;
  output out=BootMeans mean=Mean stddev=StdDev;
run;
 
/* compute distance from statistics to target values */
data Diff;
set BootMeans;     
DiffMean = &Mean - Mean;
DiffStd  = &StdDev - StdDev;
EuclidDist = euclid(DiffMean, DiffStd);  /* distance between observation and a target value */
run;
 
proc sort data=Diff;
   by EuclidDist;    /* sort by how close statistics are to target values */
run;
 
proc print data=Diff(obs=5);
   var Replicate Mean StdDev DiffMean DiffStd EuclidDist;
run;

The top five resamples (from 100) are shown. For these resamples, the mean is close to 9 and the standard deviation is close to 3, which are the target values. If you measure "closeness" by using a Euclidean distance, then Replicate=45 is the resample whose statistics are closest to the target values.

Why the example works

It is possible to turn this example into a general-purpose method for finding a resample that has statistics that are close to target values. However, it requires some additional work. I will outline why this example works and what is required to generalize the example.

Let f(x) be the probability density function of the data. If you sample uniformly from the distribution, then $E(X) = \int x f(x)\, dx$ is the expected value of the distribution. The classical bootstrap method generates samples for which the expected sample mean is E(X). However, if you sample according to some other probability distribution, g(x), then the expected value of the sample is $\int x f(x) g(x)\, dx$. By strategically choosing g(x), you can control the expected value. The same is true for higher-order moments, such as the variance.

For this example, the data has a uniform density, so it was easy to choose a value for g(x) so that the expected mean and variance matches target values. For data that has a non-uniform density, the process is more complicated because you must consider the convolution of f(x) and g(x).

Even for this example, you might not be able to match every set of target values. For example, you probably can't find a resample that has a large mean and also a large variance.

Summary

This article considers whether it is possible to find a resample that has a desired set of statistics. A classical bootstrap sample selects observations with replacement and with uniform probability. With high probability, the bootstrap sample has statistical properties that are somewhat similar to the original data. However, if you sample with nonuniform probability, you can obtain resamples whose statistical properties are quite different from the original data. The article uses the METHOD=PPS_WR option in PROC SURVEYSELECT to create an example in the special case where the data are uniformly distributed.

The post Choose samples with specified statistical properties appeared first on The DO Loop.