SLIDE 1
Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk
Axel Gandy
Department of Mathematics Imperial College London a.gandy@imperial.ac.uk
useR! 2009, Rennes July 8-10, 2009
SLIDE 2 Introduction
◮ Test statistic T, reject for large values. ◮ Observation: t. ◮ p-value:
p = P(T ≥ t) Often not available in closed form.
◮ Monte Carlo Test:
ˆ pnaive = 1 n
n
I(Ti ≥ t), where T, T1, . . . Tn i.i.d.
◮ Examples:
◮ Bootstrap, ◮ Permutation tests.
◮ Goal: Estimate p using few Xi
Mainly interested in deciding if p ≤ α for some α.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 2
SLIDE 3 Introduction
◮ Test statistic T, reject for large values. ◮ Observation: t. ◮ p-value:
p = P(T ≥ t) Often not available in closed form.
◮ Monte Carlo Test:
ˆ pnaive = 1 n
n
I(Ti ≥ t)
, where T, T1, . . . Tn i.i.d.
◮ Examples:
◮ Bootstrap, ◮ Permutation tests.
◮ Goal: Estimate p using few Xi
Mainly interested in deciding if p ≤ α for some α.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 2
SLIDE 4 Introduction
◮ Test statistic T, reject for large values. ◮ Observation: t. ◮ p-value:
p = P(T ≥ t) Often not available in closed form.
◮ Monte Carlo Test:
ˆ pnaive = 1 n
n
I(Ti ≥ t)
, where T, T1, . . . Tn i.i.d.
◮ Examples:
◮ Bootstrap, ◮ Permutation tests.
◮ Goal: Estimate p using few Xi
Mainly interested in deciding if p ≤ α for some α.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 2
SLIDE 5 Introduction
◮ Test statistic T, reject for large values. ◮ Observation: t. ◮ p-value:
p = P(T ≥ t) Often not available in closed form.
◮ Monte Carlo Test:
ˆ pnaive = 1 n
n
I(Ti ≥ t)
, where T, T1, . . . Tn i.i.d.
◮ Examples:
◮ Bootstrap, ◮ Permutation tests.
◮ Goal: Estimate p using few Xi
Mainly interested in deciding if p ≤ α for some α.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 2
SLIDE 6 Sequential approaches based on Sn = n
i=1 Xi
0 1 2 3 4 5 6 7 8 9 10 n 1 2 3 4 5 6 7 8 9 10 Sn
◮ Stop once Sn ≥ Un or
Sn ≤ Ln
◮ τ: hitting time ◮ Compute ˆ
p based on Sτ and τ.
◮ Hit BU: decide p > α, ◮ Hit BL: decide p ≤ α,
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 3
SLIDE 7 Sequential approaches based on Sn = n
i=1 Xi
0 1 2 3 4 5 6 7 8 9 10 n 1 2 3 4 5 6 7 8 9 10 Sn
◮ Stop once Sn ≥ Un or
Sn ≤ Ln
◮ τ: hitting time ◮ Compute ˆ
p based on Sτ and τ.
◮ Hit BU: decide p > α, ◮ Hit BL: decide p ≤ α,
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 3
SLIDE 8 Sequential approaches based on Sn = n
i=1 Xi
0 1 2 3 4 5 6 7 8 9 10 n 1 2 3 4 5 6 7 8 9 10 Sn
◮ Stop once Sn ≥ Un or
Sn ≤ Ln
◮ τ: hitting time ◮ Compute ˆ
p based on Sτ and τ.
◮ Hit BU: decide p > α, ◮ Hit BL: decide p ≤ α,
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 3
SLIDE 9 Previous Approaches
◮ Besag & Clifford (1991):
m n h Sn
◮ (Truncated) Sequential Probability Ratio Test, Fay et al. (2007)
m n h Sn
◮ R-package MChtest.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 4
SLIDE 10 What do we really want?
Is p ≤ α? Two individuals using the same statistical method on the same data should arrive at the same conclusion. First law of applied statistics, Gleser (1996) Consider the resampling risk RRp(ˆ p) ≡
p(ˆ
p > α) if p ≤ α, P
p(ˆ
p ≤ α) if p > α. Want: sup
p∈[0,1]
RRp(ˆ p) ≤ ǫ for some (small) ǫ > 0. For Besag & Clifford (1991), SPRT: supp RRP ≥ 0.5
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 5
SLIDE 11 What do we really want?
Is p ≤ α? Two individuals using the same statistical method on the same data should arrive at the same conclusion. First law of applied statistics, Gleser (1996) Consider the resampling risk RRp(ˆ p) ≡
p(ˆ
p > α) if p ≤ α, P
p(ˆ
p ≤ α) if p > α. Want: sup
p∈[0,1]
RRp(ˆ p) ≤ ǫ for some (small) ǫ > 0. For Besag & Clifford (1991), SPRT: supp RRP ≥ 0.5
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 5
SLIDE 12 What do we really want?
Is p ≤ α? Two individuals using the same statistical method on the same data should arrive at the same conclusion. First law of applied statistics, Gleser (1996) Consider the resampling risk RRp(ˆ p) ≡
p(ˆ
p > α) if p ≤ α, P
p(ˆ
p ≤ α) if p > α. Want: sup
p∈[0,1]
RRp(ˆ p) ≤ ǫ for some (small) ǫ > 0. For Besag & Clifford (1991), SPRT: supp RRP ≥ 0.5
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 5
SLIDE 13 What do we really want?
Is p ≤ α? Two individuals using the same statistical method on the same data should arrive at the same conclusion. First law of applied statistics, Gleser (1996) Consider the resampling risk RRp(ˆ p) ≡
p(ˆ
p > α) if p ≤ α, P
p(ˆ
p ≤ α) if p > α. Want: sup
p∈[0,1]
RRp(ˆ p) ≤ ǫ for some (small) ǫ > 0. For Besag & Clifford (1991), SPRT: supp RRP ≥ 0.5
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 5
SLIDE 14 What do we really want?
Is p ≤ α? Two individuals using the same statistical method on the same data should arrive at the same conclusion. First law of applied statistics, Gleser (1996) Consider the resampling risk RRp(ˆ p) ≡
p(ˆ
p > α) if p ≤ α, P
p(ˆ
p ≤ α) if p > α. Want: sup
p∈[0,1]
RRp(ˆ p) ≤ ǫ for some (small) ǫ > 0. For Besag & Clifford (1991), SPRT: supp RRP ≥ 0.5
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 5
SLIDE 15 What do we really want?
Is p ≤ α? Two individuals using the same statistical method on the same data should arrive at the same conclusion. First law of applied statistics, Gleser (1996) Consider the resampling risk RRp(ˆ p) ≡
p(ˆ
p > α) if p ≤ α, P
p(ˆ
p ≤ α) if p > α. Want: sup
p∈[0,1]
RRp(ˆ p) ≤ ǫ for some (small) ǫ > 0. For Besag & Clifford (1991), SPRT: supp RRP ≥ 0.5
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 5
SLIDE 16 Recursive Definition of the Boundaries
Want: sup
p RRp(ˆ
p) ≤ ǫ Suffices to ensure P
α(hit BU) ≤ ǫ
P
α(hit BL) ≤ ǫ
Recursive definition: Given U1, . . . , Un−1 and L1, . . . , Ln−1, define
◮ Un as the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ and Ln as the maximal value such that
P
α(hit BL until n) ≤ ǫn
where ǫn ≥ 0 with ǫn ր ǫ (spending sequence).
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 6
SLIDE 17 Recursive Definition of the Boundaries
Want: sup
p RRp(ˆ
p) ≤ ǫ Suffices to ensure P
α(hit BU) ≤ ǫ
P
α(hit BL) ≤ ǫ
Recursive definition: Given U1, . . . , Un−1 and L1, . . . , Ln−1, define
◮ Un as the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ and Ln as the maximal value such that
P
α(hit BL until n) ≤ ǫn
where ǫn ≥ 0 with ǫn ր ǫ (spending sequence).
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 6
SLIDE 18 Recursive Definition of the Boundaries
Want: sup
p RRp(ˆ
p) ≤ ǫ Suffices to ensure P
α(hit BU) ≤ ǫ
P
α(hit BL) ≤ ǫ
Recursive definition: Given U1, . . . , Un−1 and L1, . . . , Ln−1, define
◮ Un as the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ and Ln as the maximal value such that
P
α(hit BL until n) ≤ ǫn
where ǫn ≥ 0 with ǫn ր ǫ (spending sequence).
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 6
SLIDE 19 Recursive Definition of the Boundaries
Want: sup
p RRp(ˆ
p) ≤ ǫ Suffices to ensure P
α(hit BU) ≤ ǫ
P
α(hit BL) ≤ ǫ
Recursive definition: Given U1, . . . , Un−1 and L1, . . . , Ln−1, define
◮ Un as the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ and Ln as the maximal value such that
P
α(hit BL until n) ≤ ǫn
where ǫn ≥ 0 with ǫn ր ǫ (spending sequence).
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 6
SLIDE 20 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
k= 3 k= 2 k= 1 k= 0 1 ǫn Un 1 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 21 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 k= 3 k= 2 k= 1 .2 k= 0 1 .8 ǫn .07 Un 1 2 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 22 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 2 k= 3 k= 2 .04 k= 1 .2 .32 k= 0 1 .8 .64 ǫn .07 .11 Un 1 2 2 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 23 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 2 3 k= 3 k= 2 .04 .06 k= 1 .2 .32 .38 k= 0 1 .8 .64 .51 ǫn .07 .11 .15 Un 1 2 2 2 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 24 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 2 3 4 k= 3 k= 2 .04 .06 .08 k= 1 .2 .32 .38 .41 k= 0 1 .8 .64 .51 .41 ǫn .07 .11 .15 .18 Un 1 2 2 2 3 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 25 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 2 3 4 5 k= 3 .02 k= 2 .04 .06 .08 .14 k= 1 .2 .32 .38 .41 .41 k= 0 1 .8 .64 .51 .41 .33 ǫn .07 .11 .15 .18 .20 Un 1 2 2 2 3 3 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 26 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 2 3 4 5 6 k= 3 .02 .03 k= 2 .04 .06 .08 .14 .20 k= 1 .2 .32 .38 .41 .41 .39 k= 0 1 .8 .64 .51 .41 .33 .26 ǫn .07 .11 .15 .18 .20 .22 Un 1 2 2 2 3 3 3 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 27 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 2 3 4 5 6 7 k= 3 .02 .03 .04 k= 2 .04 .06 .08 .14 .20 .24 k= 1 .2 .32 .38 .41 .41 .39 .37 k= 0 1 .8 .64 .51 .41 .33 .26 .21 ǫn .07 .11 .15 .18 .20 .22 .23 Un 1 2 2 2 3 3 3 3 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 28 Recursive Definition - Example
◮ α = 0.2, ǫn = 0.4 n 5+n. ◮ Un=the minimal value such that
P
α(hit BU until n) ≤ ǫn ◮ Ln = maximal value such that
P
α(hit BL until n) ≤ ǫn
n =
P
α(Sn =k, τ≥n)
1 2 3 4 5 6 7 8 k= 3 .02 .03 .04 .05 k= 2 .04 .06 .08 .14 .20 .24 .26 k= 1 .2 .32 .38 .41 .41 .39 .37 .29 k= 0 1 .8 .64 .51 .41 .33 .26 .21 ǫn .07 .11 .15 .18 .20 .22 .23 .25 Un 1 2 2 2 3 3 3 3 3 Ln
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 7
SLIDE 29 Sequential Decision Procedure - Example
α = 0.2, ǫn = 0.4
n 5+n.
10 20 30 40 50 60 70 80 90 100 n 10 20
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 8
SLIDE 30 Influence of ǫ on the stopping rule
ǫ = 0.1, 0.001, 10−5, 10−7; ǫn = ǫ
n 1000+n
1000 2000 3000 4000 5000 50 100 150 200 250 300 350 n
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 9
SLIDE 31 Sequential Estimation based on the MLE
ˆ p = Sτ τ , τ < ∞ α, τ = ∞,
◮ One can show:
◮ hitting the upper boundary implies ˆ
p > α,
◮ hitting the lower boundary implies ˆ
p < α.
Hence, sup
p RRp(ˆ
p) ≤ ǫ
◮ Furthermore, ∃ random interval In s.t.
◮ In only depends on X1, . . . , Xn, ◮ ˆ
p ∈ In.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 10
SLIDE 32 Example - Two-way sparse contingency table
1 2 2 1 1 1 2 2 3 1 1 1 2 7 3 1 1 2 1 1 1 1 1
◮ H0: variables are independent. ◮ Reject for large values of the likelihood ratio test statistic T ◮ T d
→ χ2
(7−1)(5−1) under H0. Based on this: p = 0.031. ◮ Matrix sparse - approximation poor? ◮ Use parametric bootstrap based on row and column sums. ◮ Naive test statistic ˆ
pnaive with n = 1,000 replicates: p = 0.041 < 0.05. Probability of reporting p > 0.05: roughly 0.08.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 11
SLIDE 33 Example - Two-way sparse contingency table
1 2 2 1 1 1 2 2 3 1 1 1 2 7 3 1 1 2 1 1 1 1 1
◮ H0: variables are independent. ◮ Reject for large values of the likelihood ratio test statistic T ◮ T d
→ χ2
(7−1)(5−1) under H0. Based on this: p = 0.031. ◮ Matrix sparse - approximation poor? ◮ Use parametric bootstrap based on row and column sums. ◮ Naive test statistic ˆ
pnaive with n = 1,000 replicates: p = 0.041 < 0.05. Probability of reporting p > 0.05: roughly 0.08.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 11
SLIDE 34 Example - Two-way sparse contingency table
1 2 2 1 1 1 2 2 3 1 1 1 2 7 3 1 1 2 1 1 1 1 1
◮ H0: variables are independent. ◮ Reject for large values of the likelihood ratio test statistic T ◮ T d
→ χ2
(7−1)(5−1) under H0. Based on this: p = 0.031. ◮ Matrix sparse - approximation poor? ◮ Use parametric bootstrap based on row and column sums. ◮ Naive test statistic ˆ
pnaive with n = 1,000 replicates: p = 0.041 < 0.05. Probability of reporting p > 0.05: roughly 0.08.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 11
SLIDE 35
Example - Bootstrap and Sequential Algorithm
> dat <- matrix(c(1,2,2,1,1,0,1, 2,0,0,2,3,0,0, 0,1,1,1,2,7,3, 1,1,2,0,0,0,1, + 0,1,1,1,1,0,0), nrow=5,ncol=7,byrow=TRUE) > loglikrat <- function(data){ + cs <- colSums(data);rs <- rowSums(data); mu <- outer(rs,cs)/sum(rs) + 2*sum(ifelse(data<=0.5, 0,data*log(data/mu))) + } > resample <- function(data){ + cs <- colSums(data);rs <- rowSums(data); n <- sum(rs) + mu <- outer(rs,cs)/n/n + matrix(rmultinom(1,n,c(mu)),nrow=dim(data)[1],ncol=dim(data)[2]) + } > t <- loglikrat(dat); > library(simctest) > res <- simctest(function(){loglikrat(resample(dat))>=t},maxsteps=1000) > res No decision reached. Final estimate will be in [ 0.02859135 , 0.07965451 ] Current estimate of the p.value: 0.041 Number of samples: 1000 > cont(res, steps=10000) p.value: 0.04035456 Number of samples: 8574
SLIDE 36 Further Uses of the Algorithm
◮ Simulation study to evaluate whether a test is
liberal/conservative.
◮ Determining the sample size to achieve a certain power. ◮ Iterated Use:
◮ Determining the power of a bootstrap test. ◮ Simulation study to evaluate whether a bootstrap test is
liberal/conservative.
◮ Double bootstrap test. Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 13
SLIDE 37
Expected Hitting Time
Result: Ep(τ) < ∞ ∀p = α Example with α = 0.05, ǫn = ǫ
n 1000+n: 400 800 Ep(τ) ε = 0.001 ε = 1e−05 ε = 1e−07 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.2 1.4 p Ep(τ) µp
µp = theoretical lower bound on Ep(τ).
◮ Note:
1
0 µpdp = ∞; ◮ for iterated use: Need to limit the number of steps.
SLIDE 38
Expected Hitting Time
Result: Ep(τ) < ∞ ∀p = α Example with α = 0.05, ǫn = ǫ
n 1000+n: 400 800 Ep(τ) ε = 0.001 ε = 1e−05 ε = 1e−07 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.2 1.4 p Ep(τ) µp
µp = theoretical lower bound on Ep(τ).
◮ Note:
1
0 µpdp = ∞; ◮ for iterated use: Need to limit the number of steps.
SLIDE 39 Summary
◮ Sequential implementation of Monte Carlo Tests and
computation of p-values.
◮ Useful when implementing tests in packages. ◮ After a finite number of steps:
◮ ˆ
p or
◮ interval [ˆ
pL
n, ˆ
pU
n ] in which ˆ
p will lie.
◮ Guarantee (up to a very small error probability):
ˆ p is on the “correct side” of α.
◮ R-package simctest available on CRAN.
(efficient implementation with C-code)
◮ For details see Gandy (2009).
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 15
SLIDE 40 References
Besag, J. & Clifford, P. (1991). Sequential Monte Carlo p-values. Biometrika 78, 301–304. Davison, A. & Hinkley, D. (1997). Bootstrap methods and their application. Cambridge University Press. Fay, M. P., Kim, H.-J. & Hachey, M. (2007). On using truncated sequential probability ratio test boundaries for Monte Carlo implementation of hypothesis
- tests. Journal of Computational & Graphical Statistics 16, 946 – 967.
Gandy, A. (2009). Sequential implementation of Monte Carlo tests with uniformly bounded resampling risk. Accepted for publication in JASA. Gleser, L. J. (1996). Comment on Bootstrap Confidence Intervals by
- T. J. DiCiccio and B. Efron. Statistical Science 11, 219–221.
Axel Gandy Sequential Implementation of Monte Carlo Tests with Uniformly Bounded Resampling Risk 16
