R help - Apr 2009 - 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]]

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

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"

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 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