I have a relatively complicated model-fitting procedure, the output of which is a set of coefficients $\beta$ (one for each feature in the training set) and a set of hyperparameters $\lambda$ (which control the amount of regularization).

I would like to compute p-values for each feature. At the moment, my procedure is this:

  1. Determine good values for the hyperparameters $\lambda=\lambda^*$ using cross-validation.
  2. Fit the model on the whole data set with these hyperparameters, obtaining $\beta=\beta^*$
  3. For each feature $i$
    • For $j=1$ to $n$
      • Shuffle the values of feature $i$ to produce a new data set $X^j$
      • Fit the model to $X^j$ with the hyperparameters $\lambda$ to obtain $\beta^j$
    • Count the proportion of times that $|\beta^j_i| > |\beta^*_i|$ and call this the p-value for feature $i$

This matches my intuitive understanding of what a p-value is (i.e. each of the reshufflings produces a data set where the null hypotheses $\beta_i=0$ holds, since I have randomized away all the dependence) but I wondered how 'canonical' it is. I can think of a couple of other ways to do it, for example:

  1. Refit the hyperparameters as well as the coefficients on each inner loop.
  2. Use sampling with replacement to generate the new data set, instead of shuffling.

Are either of those better?

And, more importantly, can the numbers I'm generating be interpreted as p-values with any degree of reliability?

  • 1
    $\begingroup$ Since p-value is a probability, I think that it is impossible to calculate it exactly from a finite number of iterations. For example, let us observe beta-j>beta-* 2 times from 100 iterations. You are suggesting (mistakenly) that p=0.02. Becouse even if you have probability of some event 0.08, you is still able to observe only two events during 100 iterations. You still can make inference for example if you do not observe any event on 100 iterations. Or 2 events on 200 iterations - then exact p is uknown, but clearly p<0.05 $\endgroup$
    – O_Devinyak
    Aug 2, 2012 at 16:24
  • $\begingroup$ And if you'll refit lambda then beta-* cannot be compared with beta-j since then beta-j is constrained less or more $\endgroup$
    – O_Devinyak
    Aug 2, 2012 at 16:32
  • 1
    $\begingroup$ I agree with your second comment. As to the first comment - I'm not trying to compute p exactly. If I was trying to compute it exactly I wouldn't be using monte carlo. However, with $n$ large enough I can (hopefully) make the error in my approximation to p small, which is what I'm really asking about - will this method compute an approximation to the true p value or not? $\endgroup$ Aug 2, 2012 at 16:35
  • $\begingroup$ I think "yes", but let us wait for an answer of a statistician (cause I'm not) $\endgroup$
    – O_Devinyak
    Aug 2, 2012 at 16:43

1 Answer 1


Great set of alternatives for generating a discussion!

1) Your primary suggestion is actually a form of permutation test rather than a bootstrap, and does indeed generate p-values. The p-values are exact under the assumption of exchangeability, which I can't describe in terms specific to your model as I don't know your model, and if you can do all possible permutations. In a linear regression, assuming the residuals are i.i.d. is sufficient for exchangeability to hold. However, it would be proper to refit the hyperparameters inside the loop.

Note that in unbalanced ANOVA problems with missing cells it may not be possible to devise a permutation test because of the exchangeability requirement (Good, p. 138, for an example.)

2) Using sampling with replacement turns the permutation test into a nonparametric bootstrap-based test. Once again,it would be proper to refit the hyperparameters inside the loop, but exchangeability is not required.

Generally a nonparametric bootstrap test isn't exact or conservative, and it's less powerful than a permutation test, so I'd prefer the permutation test (if I have a choice.) There are also parametric bootstraps which can be used when the distribution of the test statistic under the null hypothesis is known. These are more powerful than their nonparametric cousins. In your case, if you have a linear regression and you are willing to assume i.i.d. Normal residuals, the distribution of $\hat{\beta}$ under the null hypothesis is known and you can use that fact to construct a parametric bootstrap. Still, you're not going to beat the permutation test.

A good reference for this is Permutation, Parametric, and Bootstrap Tests of Hypotheses (Good).

  • 3
    $\begingroup$ For googling, it might help that option 2 is often called the "model-based bootstrap" or "model-based resampling" with modifications like the "wild bootstrap". It is typically discussed in the context of bootstrapping linear models. $\endgroup$
    – caracal
    Aug 2, 2012 at 18:18
  • 1
    $\begingroup$ Correcting for multiple testing also sprang to mind. What is your thought about looking at the coefficient itself as the OP mentioned to calculate the p-value versus looking at the t-test (for example) from the full data fitting versus the distribution of those from each permutation? $\endgroup$
    – B_Miner
    Aug 2, 2012 at 20:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.