How to compute p-values for a bootstrap distribution

22

I was recently asked the following question:

I am using bootstrap simulations to compute critical values for a statistical test. Suppose I have test statistic for which I want a p-value. How do I compute this?

The answer to this question doesn't require knowing anything about bootstrap methods. An equivalent formulation (for a one-sided p-value) is "How do I count the number of values in a vector that are greater than a given value?" For p-values, you assume that the vector contains random values sampled from the null distribution.

You can find a fully-worked example of a bootstrap computation for a paired t test on pages 11–14 of my SAS Global Forum paper, "Rediscovering SAS/IML Software: Modern Data Analysis for the Practicing Statistician." The empirical p-value is computed on page 14.

Here's one way think about this problem. Suppose a vector, s, contains random values from the null distribution. In a bootstrap situation, this means that s1, s2, ..., sN are the bootstrapped statistics, where si is the statistic computed on the ith bootstrap sample, and where each bootstrap sample is sampled from the null distribution (that is, according to the null hypothesis). Let s0 be the value of the test statistic. Then a one-sided empirical p-value for s0 is computed as follows:

  • The simplest computation is to apply the definition of a p-value. To do this, count the number of values (statistics) that are greater than or equal to the observed value, and divide by the number of values. In code, pval = sum(s >= s0)/N;
  • The previous formula has a bias due to finite sampling. Some authors suggest the modification pval = (1+sum(s >= s0))/(N+1); For example, see Davison and Hinkley (1997), Bootstrap Methods and their Application, p. 141. Obviously, the two formulas are essentially the same when the number of values, N, is large.

See also my article on computing empirical estimates from the data.

Incidentally, if you'd like to run the bootstrap computations yourself, you can download the airlines data that I used in my SAS Global Forum paper.

Share

About Author

Rick Wicklin

Distinguished Researcher in Computational Statistics

Rick Wicklin, PhD, is a distinguished researcher in computational statistics at SAS and is a principal developer of SAS/IML software. His areas of expertise include computational statistics, simulation, statistical graphics, and modern methods in statistical data analysis. Rick is author of the books Statistical Programming with SAS/IML Software and Simulating Data with SAS.

22 Comments

  1. Jitendra Kumar on

    i want to find Bootstrap-p and Bootstrap -t.
    so i want a programme to find that values.

    Thank you.

  2. I think the intent of the question is how to obtain p-values for a test when we have the original data (in which we don't know whether the effect is null or not), and do not already have data from the null distribution. Could you modify the bootstrap procedure to get p-values that are consistent with the 95% CIs from the bootstrap? (Eg. if the CI straddles 0 the p-value is >0.05.) Randomly permuting the treatment labels works for p-values but doesn't give CIs.

    • Rick Wicklin

      You ask two questions. The first concerns hypothesis tests. The second concerns CIs.
      1) I am not assuming that the data follows the null distribution. To apply the bootstrap,you have to choose a resampling scheme. When testing a hypothesis, you should resample AS IF the hypothesis is true. You can use the resulting bootstrap distribution to estimate the p-value.
      2) Yes, if you have a parameter estimate from the data you can use a bootstrap technique to form an approximate CI for the parameter. Unfortunately, the technique might depend on the estimate. For example, bootstrapping a median requires a different technique than bootstrapping a mean.

  3. Christian Hinze on

    i was wondering how many instances of the null hypothesis i would need for a pvalue of say 0.05.

    i could for instance only compute 10 instances of my null hypthesis and apply the above formula (pval = sum(s >= s0)/N). It is maybe 0. So, how many instance (N) do I need for a 0.05 or 0.01 pvalue?

    Thanks a lot!

    • Rick Wicklin

      That is an excellent question. Unfortunately, there is no simple answer. If the sampling distribution is normal, you might only need hundreds of simulations; if the distribution is severely skewed or long-tailed, more are needed. In my book Simulating Data with SAS I discuss this problem and make several additional recommendations.

      • M Amir Saeed on

        respected sir can you send me this book Simulating Data with SAS. i cant purchase it. please send me the bootstrap simulation topic.

  4. Pingback: Resampling and permutation tests in SAS - The DO Loop

  5. Thanks a lot for this post. just to clarify, the formula you have provided: pval = sum(s >= s0)/N is for a 2-sided p-value and if I wanted a one-sided p-value I would instead use pval = sum(s > s0)/N or pval = sum(s < s0)/N?

  6. KAUSHIK BHATTACHARJE on

    Hi Rick,
    I am trying to calculate the p-values for all the values of my original sample(m=53) .My bootstrapped estimates are (N=)5000. for that I am writing the following code :

    data Ave; set stacked (keep= Mean); run;

    %let Mydata1=Ave; * Contain Bootstrapped estimate .Dim is 5000x1.Variable name is Mean;
    %let MyData2 = pcm; * Contain actual sample.Dim is 53x1. Variable name is pcm;

    proc iml;

    use &Mydata1; read all var {Mean} into y; close &MyData1;

    use &MyData2; read all var {pcm} into x; close &MyData2;

    print x; * Just to see;

    N=nrow(y) ;

    m=nrow(x);

    print m N;

    P=j(m, 1,1); *mx1 vector of 1s;

    Do i = 1 to m ;

    a=x[i,1] ; print a;

    idx = loc(y> a); print idx; p = nrow(idx`); * how many?;

    print p N; p_val= p/N; print p_val;

    P[i]= p_val ;

    end;

    quit;
    after execution I am getting

    ERROR: Matrix idx has not been set to a value.

    Can you help me to fix the error?

    Kaushik

      • KAUSHIK BHATTACHARJE on

        ok & thanks .
        I will do that .
        However I still think I should ask you here beacuse it contains& related to the DO Loop(my be I am wrong): The point is , as I understand from your blogs & books that the LOC function is dynamic in nature and there is no need to initialize it.It is working fine idx= loc(y>8) say.When I am putting that into a loop i.e. idx=loc(y>x[i,1]) , it is giving error message : Matrix idx has not been set to a value.
        Can you give any clue as to where I am wrong?

        • Rick Wicklin

          Yes. If the condition inside the LOC function is not satisfied, then the LOC function returns an empty matrix. When you try to use idx as a subscript, the error occurs. You can read about this situation and how to handle it correctly in the article "Beware the naked LOC."

          • KAUSHIK BHATTACHARJE on

            Thanks Rick. I am able to figure out now. Though I have been following your blogs last couple of years this is the first time I have interacted with you. I have learnt many a things from your blogs : not about SAS only sometimes about statistics too. Thanks again

  7. What if the observed beta could be bigger or smaller than the bootstrap value. How would we calculate a 2 sided p value?

    • The computation doesn't change. For the first definition, the bootstrap estimate of the probability can be 0. You should display the estimate by using the PVALUE format, such as
      print pval[format=pvalue6.];
      which will display the text "<.0001".
      For the second estimate, the p-value will always be nonzero.

  8. Anders Møller Greve on

    I have a question. I computed nonparametric bootstrap hazard ratios for the variable X (X=0 vs. X=1) using 250 replications (proc surveyselect), which are given below:

    CI95_lower CI95_median CI95_upper
    0.66051 0.90034 1.23374

    There seems to be no difference in rates of the investigated endpoint as a function of X. When I try to calculate the p-value for 1 being included (no difference between X=0 and X=1) in the bootstrap confidence interval, I get the p-values below:

    N lt1 gt1
    250 0.692 0.308

    Here I would choose P=0.31 for a protective effect of X on the investigated endpoint (hazard ratio <1 for X=1 vs. X=0). But what if I had no knowledge about the direction of the effect? How can I know, which of the one-sided p-values I should use?

    Thanks!

    • Rick Wicklin

      Yes, it can be confusing because there are so many symmetries. For example, do you model the event as X=1 or as X=0? Are you interested in the ratio P(Y|X=0)/P(Y|X=1) or the inverse ratio? In practice, we are guided by domains-specific knowledge such as assuming that the experimental treatment will have a certain effect compared to the control group. That domain knowledge helps you to write down the null and alternative hypotheses.

Leave A Reply

Back to Top