Johan Jackson
2009-Apr-12 21:47 UTC
[R] p-values from bootstrap - what am I not understanding?
Dear stats experts: Me and my little brain must be missing something regarding bootstrapping. I understand how to get a 95%CI and how to hypothesis test using bootstrapping (e.g., reject or not the null). However, I'd also like to get a p-value from it, and to me this seems simple, but it seems no-one does what I would like to do to get a p-value, which suggests I'm not understanding something. Rather, it seems that when people want a p-value using resampling methods, they immediately jump to permutation testing (e.g., destroying dependencies so as to create a null distribution). SO - here's my thought on getting a p-value by bootstrapping. Could someone tell me what is wrong with my approach? Thanks: STEPS TO GETTING P-VALUES FROM BOOTSTRAPPING - PROBABLY WRONG: 1) sample B times with replacement, figure out theta* (your statistic of interest). B is large (> 1000) 2) get the distribution of theta* 3) the mean of theta* is generally near your observed theta. In the same way that we use non-centrality parameters in other situations, move the distribution of theta* such that the distribution is centered around the value corresponding to your null hypothesis (e.g., make the distribution have a mean theta = 0) 4) Two methods for finding 2-tailed p-values (assuming here that your observed theta is above the null value): Method 1: find the percent of recentered theta*'s that are above your observed theta. p-value = 2 * this percent Method 2: find the percent of recentered theta*'s that are above the absolute value of your observed value. This is your p-value. So this seems simple. But I can't find people discussing this. So I'm thinking I'm wrong. Could someone explain where I've gone wrong? J Jackson [[alternative HTML version deleted]]
Peter Dalgaard
2009-Apr-12 22:38 UTC
[R] p-values from bootstrap - what am I not understanding?
Johan Jackson wrote:> Dear stats experts: > Me and my little brain must be missing something regarding bootstrapping. I > understand how to get a 95%CI and how to hypothesis test using bootstrapping > (e.g., reject or not the null). However, I'd also like to get a p-value from > it, and to me this seems simple, but it seems no-one does what I would like > to do to get a p-value, which suggests I'm not understanding something. > Rather, it seems that when people want a p-value using resampling methods, > they immediately jump to permutation testing (e.g., destroying dependencies > so as to create a null distribution). SO - here's my thought on getting a > p-value by bootstrapping. Could someone tell me what is wrong with my > approach? Thanks: > > STEPS TO GETTING P-VALUES FROM BOOTSTRAPPING - PROBABLY WRONG: > > 1) sample B times with replacement, figure out theta* (your statistic of > interest). B is large (> 1000) > > 2) get the distribution of theta* > > 3) the mean of theta* is generally near your observed theta. In the same way > that we use non-centrality parameters in other situations, move the > distribution of theta* such that the distribution is centered around the > value corresponding to your null hypothesis (e.g., make the distribution > have a mean theta = 0) > > 4) Two methods for finding 2-tailed p-values (assuming here that your > observed theta is above the null value): > Method 1: find the percent of recentered theta*'s that are above your > observed theta. p-value = 2 * this percent > Method 2: find the percent of recentered theta*'s that are above the > absolute value of your observed value. This is your p-value. > > So this seems simple. But I can't find people discussing this. So I'm > thinking I'm wrong. Could someone explain where I've gone wrong?There's nothing particularly wrong about this line of reasoning, or at least not (much) worse than the calculation of CI. After all, one definition of a CI at level 1-alpha is that it contains values of theta0 for which the hypothesis theta=theta0 is accepted at level alpha. (Not the only possible definition, though.) The crucial bit in both cases is the assumption of approximate translation invariance, which holds asymptotically, but maybe not well enough in small samples. There are some braintwisters connected with the bootstrap; e.g., if the bootstrap distribution is skewed to the right, should the CI be skewed to the right or to the left? The answer is that it cannot be decided based on the distribution of theta* alone since that depends only on the true theta, and we need to know what the distribution would have been had a different theta been the true one. The point is that these things get tricky, so most people head for the safe haven of permutation testing, where it is rather more easy to feel that you know what you are doing. For a rather different approach, you might want to look into the theory of empirical likelihood (book by Art Owen, or just Google it). -- O__ ---- Peter Dalgaard ?ster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
Robert A LaBudde
2009-Apr-13 00:22 UTC
[R] p-values from bootstrap - what am I not understanding?
There is really nothing wrong with this approach, which differs primarily from the permutation test in that sampling is with replacement instead of without replacement (multinomial vs. multiple hypergeometric). One of the issues that permutation tests don't have is bias in the statistic. In order for bootstrap p-values to be reasonably accurate, you need a reasonable dataset size, so that sampling with replacement isn't a big effect, and so that enough patterns arise in resampling. It also helps if the data is continuous instead of categorical or binary. The same issues affect permutation tests, but untroubled by bias. The usual methods for p-values (e.g., see Fisher's test in Agresti's Categorical Analysis) work here. Typically there is some ambiguity on how to treat the values equal to the observed statistic. If you include it, the p-value is conservative for rejection. If you don't, it's liberal for rejection. If you include 1/2 weight, it averages correctly in the long run. Ditto for 2-tailed p-values vs. single tails. Several different methods (some of which you listed) are used. As a general rule, if you have data from which you wish a p-value, a permutation (i.e., without replacement) test is used, but for confidence intervals, bootstrapping (i.e., with replacement) is used. For reasonably large datasets, both methods will agree closely. But permutation tests are typically used for smaller size datasets. (Think binomial vs. hypgeometric distributions for p-values, and when they agree.) At 05:47 PM 4/12/2009, Johan Jackson wrote:>Dear stats experts: >Me and my little brain must be missing something regarding bootstrapping. I >understand how to get a 95%CI and how to hypothesis test using bootstrapping >(e.g., reject or not the null). However, I'd also like to get a p-value from >it, and to me this seems simple, but it seems no-one does what I would like >to do to get a p-value, which suggests I'm not understanding something. >Rather, it seems that when people want a p-value using resampling methods, >they immediately jump to permutation testing (e.g., destroying dependencies >so as to create a null distribution). SO - here's my thought on getting a >p-value by bootstrapping. Could someone tell me what is wrong with my >approach? Thanks: > >STEPS TO GETTING P-VALUES FROM BOOTSTRAPPING - PROBABLY WRONG: > >1) sample B times with replacement, figure out theta* (your statistic of >interest). B is large (> 1000) > >2) get the distribution of theta* > >3) the mean of theta* is generally near your observed theta. In the same way >that we use non-centrality parameters in other situations, move the >distribution of theta* such that the distribution is centered around the >value corresponding to your null hypothesis (e.g., make the distribution >have a mean theta = 0) > >4) Two methods for finding 2-tailed p-values (assuming here that your >observed theta is above the null value): >Method 1: find the percent of recentered theta*'s that are above your >observed theta. p-value = 2 * this percent >Method 2: find the percent of recentered theta*'s that are above the >absolute value of your observed value. This is your p-value. > >So this seems simple. But I can't find people discussing this. So I'm >thinking I'm wrong. Could someone explain where I've gone wrong? > > >J Jackson > > [[alternative HTML version deleted]] > >______________________________________________ >R-help at r-project.org mailing list >https://stat.ethz.ch/mailman/listinfo/r-help >PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >and provide commented, minimal, self-contained, reproducible code.===============================================================Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: ral at lcfltd.com Least Cost Formulations, Ltd. URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239 Fax: 757-467-2947 "Vere scire est per causas scire"
Viechtbauer Wolfgang (STAT)
2009-Apr-15 21:58 UTC
[R] p-values from bootstrap - what am I not understanding?
Some useful comments have already been made. I would like to comment on the two definitions of the p-value under (4) -- since I thought exactly about this issue a while back. Maybe this will be useful ... Suppose the distribution of a test statistic Z under H0 is given by f(Z) and that the distribution is not symmetric. A bootstrap distribution (under H0) may in fact look like that. We can define the critical values for alpha = .05 in two different ways: (1) We find a single value of z so that the sum of the lower and upper tail is equal to .05 (red shaded region, corresponding to z = 9.51 and z = -9.51). Since the distribution is very skewed, there is no area in the lower tail, so that the red-shaded region in the upper tail is actually equal to .05. (2) We put .025 in the lower and .025 in the upper tail (blue shaded regions, corresponding to z = -3.82 and z = 11.53 in the figure). The first definition treats deviations on both ends of the sampling distribution equally. The second definition acknowledges the fact that deviations of a certain magnitude are more likely to occur in the upper than in the lower tail. Hence, the critical value on the upper tail is more extreme than on the lower tail. This is, in fact, how critical values are typically defined for non-symmetric distributions. Now given some data, we observe the value z of the test statistic. How should we now calculate the two-sided p-value? option 1: p = prob(|Z| > |z|) ----------------------------- If z = 11.53, then p = .025 (reject H0) If z = -11.53, then p = .025 (reject H0) If z = 9.51, then p = .050 (reject H0) If z = -9.51, then p = .050 (reject H0) If z = -3.82, then p = .303 (do not reject H0) If z = 1, then p = .779 Note that one could never actually get a p-value below .05 if z is negative. option 2: p = 2*prob(Z > z) if z is positive or 2*prob(Z < z) if z is negative --------------------------------------------------------------------- If z = 11.53, then p = .050 (reject H0) If z = -11.53, then p = .000 (reject H0) If z = 9.51, then p = .100 (do not reject H0) If z = -9.51, then p = .000 (reject H0) If z = -3.82, then p = .050 (reject H0) If z = 1, then p = 1.07 !!! For z = 9.51, we should not reject H0 according to the second definition of the critical values. On the other hand, for z = -3.82, we should reject H0 according to the second definition. So option 2 is more in line how we typically define critical values in non-symmetric distributions. However, note that the p value can be above 1 according to the second definition! ###################################################################### alpha <- .05 shape <- 4 scale <- 2 xs <- seq(0, 30, .1) #ys <- dchisq(xs, df=4) ys <- dgamma(xs, shape=shape, scale=scale) offset <- xs[which( ys == max(ys) )] xs <- xs - offset par(mar=c(5,4,2,2)) plot(xs, ys, type="l", lwd=2, xlim=c(-15,25), xlab="Z", ylab="f(Z)") abline(v = 0) ######################################################################## ##### ### some p-value calculations round( pgamma(3.82+offset, shape=shape, scale=scale, lower.tail=F) + pgamma(-3.82+offset, shape=shape, scale=scale, lower.tail=T) , 3 ) round( pgamma(0+offset, shape=shape, scale=scale, lower.tail=F) + pgamma(-0+offset, shape=shape, scale=scale, lower.tail=T) , 3 ) round( pgamma(1+offset, shape=shape, scale=scale, lower.tail=F) + pgamma(-1+offset, shape=shape, scale=scale, lower.tail=T) , 3 ) 2*round( pgamma(1+offset, shape=shape, scale=scale, lower.tail=F) , 3 ) ######################################################################## ##### upp.crit <- qgamma(alpha, shape=shape, scale=scale, lower.tail=F) x.r <- seq(upp.crit, max(xs)+offset, length = 100) y.r <- c( dgamma(x.r, shape=shape, scale=scale), 0, 0) x.r <- c(x.r, max(xs), upp.crit) - offset polygon(x.r, y.r, density = 50, col="red", angle=135) text(upp.crit-offset, .03, paste("z = ", round(upp.crit-offset,2)), pos=4, offset=0, col="red") text(-1*(upp.crit-offset), .03, paste("z = ", -1*round(upp.crit-offset,2)), pos=2, offset=0, col="red") low.crit <- qgamma(alpha/2, shape=shape, scale=scale, lower.tail=T) x.l <- seq(min(xs), low.crit, length = 100) y.l <- c( dgamma(x.l, shape=shape, scale=scale) , 0, 0) x.l <- c(x.l, low.crit, min(xs)) - offset polygon(x.l, y.l, density = 50, col="blue") text(low.crit-offset, .05, paste("z = ", round(low.crit-offset,2)), pos=2, offset=0, col="blue") upp.crit <- qgamma(alpha/2, shape=shape, scale=scale, lower.tail=F) x.r <- seq(upp.crit, max(xs)+offset, length = 100) y.r <- c( dgamma(x.r, shape=shape, scale=scale), 0, 0) x.r <- c(x.r, max(xs), upp.crit) - offset polygon(x.r, y.r, density = 50, col="blue") text(upp.crit-offset, .05, paste("z = ", round(upp.crit-offset,2)), pos=4, offset=0, col="blue") abline(h = 0) ###################################################################### I hope this helps! Best, -- Wolfgang Viechtbauer Department of Methodology and Statistics University of Maastricht, The Netherlands http://www.wvbauer.com/ ----Original Message---- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Johan Jackson Sent: Sunday, April 12, 2009 23:48 To: r-help at r-project.org Subject: [R] p-values from bootstrap - what am I not understanding?> Dear stats experts: > Me and my little brain must be missing something regardingbootstrapping.> I understand how to get a 95%CI and how to hypothesis test using > bootstrapping (e.g., reject or not the null). However, I'd also liketo> get a p-value from it, and to me this seems simple, but it seemsno-one> does what I would like to do to get a p-value, which suggests I'm not > understanding something. Rather, it seems that when people want ap-value> using resampling methods, they immediately jump to permutation testing > (e.g., destroying dependencies so as to create a null distribution).SO -> here's my thought on getting a p-value by bootstrapping. Could someone > tell me what is wrong with my approach? Thanks: > > STEPS TO GETTING P-VALUES FROM BOOTSTRAPPING - PROBABLY WRONG: > > 1) sample B times with replacement, figure out theta* (your statisticof> interest). B is large (> 1000) > > 2) get the distribution of theta* > > 3) the mean of theta* is generally near your observed theta. In thesame> way that we use non-centrality parameters in other situations, movethe> distribution of theta* such that the distribution is centered aroundthe> value corresponding to your null hypothesis (e.g., make thedistribution> have a mean theta = 0) > > 4) Two methods for finding 2-tailed p-values (assuming here that your > observed theta is above the null value): Method 1: find the percent of > recentered theta*'s that are above your observed theta. p-value = 2 * > this percent Method 2: find the percent of recentered theta*'s thatare> above the absolute value of your observed value. This is your p-value.> > So this seems simple. But I can't find people discussing this. So I'm > thinking I'm wrong. Could someone explain where I've gone wrong? > > > J Jackson