A new permutation test statistic for complete block designs
aa r X i v : . [ m a t h . S T ] N ov The Annals of Statistics (cid:13)
Institute of Mathematical Statistics, 2015
A NEW PERMUTATION TEST STATISTIC FOR COMPLETEBLOCK DESIGNS
By Inga Samonenko and John Robinson University of Sydney
We introduce a nonparametric test statistic for the permutationtest in complete block designs. We find the region in which the statis-tic exists and consider particularly its properties on the boundary ofthe region. Further, we prove that saddlepoint approximations for tailprobabilities can be obtained inside the interior of this region. Finally,numerical examples are given showing that both accuracy and powerof the new statistic improves on these properties of the classical F -statistic under some non-Gaussian models and equals them for theGaussian case.
1. Introduction.
Randomized designs and permutation tests originatedin the work of Fisher (1935). Kolassa and Robinson (2011) obtained theo-rems on the distribution of a general likelihood ratio like statistic under weakconditions and applied these to the one-way or k -sample permutation tests,obtaining saddlepoint approximations generalizing the Lugananni–Rice andBarndorff–Nielsen approximations for one-dimensional means. Here, we usetheir general result and apply their approach to permutation tests for com-plete block designs, paying particular attention to the region in which thestatistic exists and in the interior of which saddlepoint approximations canbe obtained. This interior is the admissible domain, following Borovkov andRogozin (1965). We examine the properties of the test statistic in this regionand on its boundary, and obtain results on the relative errors of saddlepointapproximations inside the admissible domain. We also give numerical re-sults for comparisons of the new statistic with the commonly used F -statisticwhich demonstrate the accuracy of the saddlepoint approximation and show,for long tailed error distributions, an improvement in power relative to the F -statistic with no loss of power for near normal errors. Received December 2013; revised April 2014. Supported by UPA Scholarship. Supported by ARC DP0773345.
AMS 2000 subject classifications.
Primary 62G09, 62G10, 62G20; secondary 60F10.
Key words and phrases.
Saddlepoint approximations, admissible domain.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Statistics ,2015, Vol. 43, No. 1, 90–101. This reprint differs from the original in paginationand typographic detail. 1
I. SAMONENKO AND J. ROBINSON
A randomized complete block design is used to compare the effect of k different treatments in b blocks, usually selected to reduce the variationwithin subunits of the block. The analysis of variance is used to test the nullhypothesis that the treatments have the same effect, with the test statistic F , the ratio of the treatment and error mean squares. Under the assumptionthat the errors are normally distributed, the null distribution of F is the F distribution with k − k − b − F testis equivalent to an unconditional likelihood ratio test.The random assignment of k treatments to each block allows us to usea permutation test based on means which is distribution-free and does notrely either on the assumption of normality or on asymptotics. This test canbe performed using the F -statistic and will have correct size, conditionallyon the order statistics in each block, and so unconditionally, for any dis-tribution of errors under the null hypothesis of no treatment effects in astandard two-way model or for a model based on randomization prior to theexperiment. Under the null hypothesis, the permutation distribution of thisstatistic can be calculated exactly by evaluating all possible values of the teststatistic under permutations in each block and taking these as equi-probable.When this is numerically infeasible, Monte Carlo methods are widely usedto approximate the exact distribution by using a large random sample ofthe possible permutations. A chi-squared distribution with ( k − F distribution with ( k −
1) and ( k − b − F , based on exponen-tial tilting. We show that this statistic can be calculated on the admissibledomain, an open convex set, the closure of which contains the support of thetreatment means. We consider the boundary of the admissible domain andshow that the statistic can be obtained on the boundary as a limit whichcan be calculated using lower dimensional versions of the statistic on lowerdimensional versions of the admissible domain. We then obtain saddlepointapproximations for the tail probability of this statistic with relative errorsof order 1 /n in the admissible domain, based on Theorems of Kolassa andRobinson (2011). The results generalize the saddlepoint approximations ofRobinson (1982) in the case of permutation tests of paired units, which canbe regarded as a block design with blocks of size 2, where the admissibledomain is the interval between the mean of the absolute values of differencesof the pairs and its negative. ERMUTATION TESTS FOR COMPLETE BLOCK DESIGN In the next section, we introduce the notation for a complete block design,obtain the likelihood ratio like statistic and define its admissible domain. InSection 3, we describe the admissible domain and give three theorems givingexplicit results for the test statistic on the boundary of the domain, withproofs given in Section 6. In Section 4, we use the theorems of Kolassa andRobinson (2011) to show that tail probabilities for the statistic under per-mutations can be approximated in the admissible domain by an integral ofa formal saddlepoint density given in forms like those of Lugananni–Riceand Barndorff–Nielsen in the one-dimensional case. In Section 5, we presentnumerical calculations illustrating the accuracy of the approximations com-pared to those obtained using the standard test statistics and give powercomparisons showing an improvement in power over the standard F -test forobservations from long tailed distributions. The code used is available from .
2. The test statistic Λ and its admissible domain. Let ( X ij ) be a matrixof observed experimental values, normalized to have row means zero, where i = 1 , . . . , b is the block number and j = 1 , . . . , k is the treatment number. Letmatrix A = ( a ij ) have rows of the matrix ( X ij ) each set in ascending orderand let A i be its i th row. Define the means ¯ X j = P bi =1 X ij /b for j = 1 , . . . , k and let ¯ X = ( ¯ X , . . . , ¯ X k − ) T . Then, given A , under the null hypothesis ofequal treatment effects, the conditional cumulative generating function fortreatment means is bκ ( τ ) = log E ( e P kj =1 τ j ¯ X j | A ) = b X i =1 log E ( e P kj =1 τ j X ij /b | A i ) . Set t i = ( τ i − τ k ) /b for i = 1 , . . . , k −
1. Then we can reduce the problem ofdefining the average cumulative generating function to a ( k − κ ( t ) = 1 b b X i =1 log 1 k ! X π ∈ Π e t T a iπ , (1)where t = ( t , . . . , t k − ) T , Π is the set of possible vectors ( π (1) , . . . , π ( k − k − { , . . . , k } and a iπ = ( a iπ (1) , . . . , a iπ ( k − ) T .Consider the test statistic Λ( ¯ X ), whereΛ( x ) = sup t { t T x − κ ( t ) } , (2)for t, x ∈ R k − . Let us define an admissible domain Ω ⊂ R k − as the set ofall x for which t T x − κ ( t ) attains its maximum. Then there exists a uniquevalue t x such thatΛ( x ) = t Tx x − κ ( t x ) if and only if κ ′ ( t x ) = x, (3) I. SAMONENKO AND J. ROBINSON since t T x − κ ( t ) is strictly convex by noting that − κ ′′ ( t ) is negative definite.In the case k = 2, the admissible domain is ( − P bi =1 | a i − a i | , P bi =1 | a i − a i | ) and the properties of Λ and the saddlepoint approximation are dis-cussed for the two special cases of the binomial and the Wilcoxon signed-rank test in Jin and Robinson (1999). The situation is more complex for k > Remark.
Exact randomization tests have restricted application to de-signed experiments. The only two designs for which we know how to obtaina statistic of our form are the complete block design considered here andthe one-way or k -sample design considered in Kolassa and Robinson (2011).An extension to some other cases such as balanced incomplete block designsor in testing for main effects using restricted randomization as suggested byBrown and Maritz (1982) may be possible but do not seem to be straight-forward.
3. The properties of Λ . First, we will describe the admissible domainand give some results which make it possible to calculate Λ( x ) on the bound-ary of the domain where the solution of the saddlepoint equations (3) doesnot exist. Let ¯ A = ( ¯ A , . . . , ¯ A k ) T be a vector of column means of A and write¯ A π = ( ¯ A π (1) , . . . , ¯ A π ( k − ) T , for any π ∈ Π. Then the support of ¯ X contains¯ A π and the set of ¯ A π , for all π ∈ Π, is the set of extreme points of the convexhull of the support of ¯ X , which is a ( k − P . Theorem 1.
The set Ω is the interior of the ( k − -polytope P . Theorem 2.
The function Λ( x ) is finite on the boundary of Ω and takesits maximum value log k ! at its extreme points. The boundary of P consists of all x ∈ P for which there exists an integer l and distinct integers s , . . . , s k − from the set { , . . . , k − } satisfying oneof the equalities l X j =1 x s j = l X j =1 ¯ A j or l X j =1 x s j = l X j =1 ¯ A k − j +1 . (4) Theorem 3.
On the boundary of Ω corresponding to the value l we have Λ( x ) = Λ ( x ) + Λ ( x ) + log (cid:18) kl (cid:19) , for Λ ( x ) = sup u ,...,u l − l − X j =1 x s j u j − b b X i =1 log 1 l ! X ˆ π ∈ ˆΠ e P l − j =1 a i ˆ π j ) u j ! ERMUTATION TESTS FOR COMPLETE BLOCK DESIGN Fig. 1.
Examples of the admissible domain in the cases k = 3 and k = 4 . and Λ ( x ) = sup u l +1 ,...,u k − k − X j = l +1 x s j u j − b b X i =1 log 1( k − l )! X ˆ π ∈ ˆΠ e P k − j = l +1 a i ˆ π j ) u j ! , where ˆΠ and ˆΠ are sets of all permutations ˆ π and ˆ π of integers { , . . . , l } and { l + 1 , . . . , k } , respectively. Remark.
The result of Theorem 3 demonstrates that the boundary ofΩ ⊂ R k − consists of lower dimensional polytopes, each made up of a crossproduct of two sets of dimension l − k − l −
1, for l = 1 , . . . , k − l − and Λ are defined on these subsets as is Λ in (2). Toillustrate this, in Figure 1 we have given two diagrams showing the polytope P for the cases k = 3 and k = 4. In the first picture, we have 6 vertexesand 6 sides with boundaries made up of lines representing the dimensionreduction to one dimension. In this case, one of Λ and Λ is identicallyzero. In the second picture, the two-dimensional boundaries are either six-sided, corresponding to one of Λ and Λ being identically zero, and theother a two-dimensional function, or are rectangles corresponding to bothΛ and Λ being one-dimensional functions.
4. Saddlepoint approximations for Λ . Consider P(Λ( ¯ X ) ≥ u / A , and define r ( x ) = e − b Λ( x ) (2 π/b ) − k/ | κ ′′ ( t x ) | − / I. SAMONENKO AND J. ROBINSON for x ∈ Ω. In the case of identically distributed random vectors X i = ( X i ,. . . , X i,k − ) ∈ Ω, i = 1 , . . . , N with known densities, r ( ¯ X ) is a saddlepointdensity approximation for ¯ X , obtained by Borovkov and Rogozin (1965).In our case, the lack of a density requires the application of Theorem 1 ofKolassa and Robinson (2011). We considerP(Λ( ¯ X ) ≥ u /
2) = P( ¯ X ∈ F ) , where F = { x : Λ( x ) ≥ u / } and u ∈ Ω − ε = (Ω cε ) c = { x ∈ R k − , y ∈ Ω c : | x − y | < ε } c . Whenever F c ⊂ Ω − ε , our case must only meet the necessary con-ditions (A1)–(A4) stated in Kolassa and Robinson (2011). The cumulativegenerating function (1) exists throughout R k − , therefore, the first conditionis met. The average variance κ ′′ ( t ) is a positive definite matrix which equalsthe identity matrix at the origin. Thus, the second condition is met. Thethird condition only requires the existence of some moments and the fourthis a smoothness condition, which we assume holds. It will hold, for example,if the observations are from a distribution with a continuous component.Thus we can apply Theorems 1 and 2 of Kolassa and Robinson (2011) as inthat paper to get the following result. Theorem 4.
For ε > and u / < log k − ε , under conditions (A1)–(A4) of Kolassa and Robinson (2011), P (Λ( ¯ X ) ≥ u /
2) = Q k − ( bu )(1 + O (1 /b )) + c b b u k − e − bu / G ( u ) − u , (5) P (Λ( ¯ X ) ≥ u /
2) = Q k − ( bu ∗ )(1 + O (1 /b )) , (6) where Q k − ( x ) = P( X k − ≥ x ) = 12 ( k − / Γ(( k − / Z ∞ x z ( k − / − e − z/ dz,u ∗ = u − log( G ( u )) /bu , c b = b ( k − / / ( k − / Γ( k − ) , δ ( u, s ) = Γ(( k − / | κ ′′ ( t x ) | − / | κ ′′ (0) | / r k − π ( k − / u k − | s T κ ′′ (0) / t x | and G ( u ) = Z S k − δ ( u, s ) ds, for S k − the k − -dimensional unit sphere centered at zero, and where, foreach s ∈ S k − , r is chosen so Λ( rκ ′′ (0) / s ) = u / . Here, t x is a solutionto (3) at the point x = rκ ′′ (0) / s . ERMUTATION TESTS FOR COMPLETE BLOCK DESIGN Table 1
Accuracy for exponentially squared errors, b = 10 and k = 4 u MC F F We note that the constraint u / < log k − ε , ensures that the level setof Λ( x ) corresponding to u lies entirely in Ω, since the minimum value ofΛ( x ) for an x on the boundary occurs for l = 1 and Λ ( x ) = Λ ( x ) = 0 inTheorem 3. The remainder of the proof then follows in the same way asin Theorem 2 of Kolassa and Robinson (2011), so we omit it. The theoremgives approximations of the tail probabilities of the test statistic Λ underpermutations, in forms like those of Lugananni–Rice and Barndorff–Nielsenin the one-dimensional case.
5. Numerical results.
Accuracy.
For each of the simulation experiments, we obtained asingle matrix A by sampling from a distribution, that of squared exponentialrandom variables for our Tables 1, 2 and 3. Then we used 100,000 replicatesof random permutations of each block to obtain Monte Carlo approxima-tions to the tail probabilities of the permutation tests for the statistics F and Λ, shown as MC F and MC Λ in the tables. We compared these to thetail probabilities from the F distribution for the F -statistic and to the sad-dlepoint approximations for the Λ statistic obtained using formulas (5) (SPLR) and (6) (SP BN), respectively, with 100 Monte Carlo samples used toapproximate integrals on the sphere S k − , as in the Remark in Section 2 of Table 2
Accuracy for exponential squared errors, b = 5 and k = 3 u MC F F I. SAMONENKO AND J. ROBINSON
Table 3
Average of 1000 permutation test results for exponential squared errors with b = 10 and k = 4 compared to the F distribution u E MC F F Kolassa and Robinson (2011). We also used the method from Genz (2003),for approximation of the integral on the sphere, obtaining effectively thesame accuracy as with Monte Carlo sampling.From Tables 1 and 2, we note that the accuracy is high for the Λ test,even for only 5 blocks of size 3. We note that for Table 2 the theorem holdsfor u less than √ .
48, so we are restricted to this region. Resultsfrom other simulations show even greater accuracy under normal errors orerrors that are not from long tailed distributions.The F -statistic has less accuracy in the tails, partly because the F -statistic approximates the average of tail probabilities conditioned on thematrix A , using 100,000 permutations for each A , so that even in the caseof normal errors, it may not agree with the conditional tail probabilitiesapproximated by MC F -values from the tables of this section, which giveproportions in the tails obtained from 100,000 Monte Carlo simulations froma particular sample and is an approximation of the conditional distribution.To consider the accuracy of the unconditional test, we obtained 1000 repli-cates from each of a normal and exponential squared distribution, obtainedtail probabilities for these from the permutation test for the F -statistic, av-eraged these over the 1000 replicates and compared these approximations tothe F distribution. In the normal case, the results were very accurate, es-sentially replicating the theoretical results, as expected, and for the squaredexponential case the results are given in Table 3 indicating that errors remainunsatisfactory in the tails.5.2. Power results for the saddlepoint approximations.
We compare the F -statistic and the saddle point approximations using (5) and (6) using 100Monte Carlo uniform samples from S k . There were 2000 samples with er-rors drawn from the exponential and the exponential squared distributions,and for each of these p -values were calculated using the saddlepoint approx-imations for the Λ statistic obtained using formulas (5) (PowerLR) and (6)(PowerBN), respectively, and using 10,000 permutations to approximate the p -values for the F -statistic, for a design with 10 blocks of size 4. We selectedtreatment effects µ , as (0 , , , , ( − / , , , / , . . . , ( − / , , , / F -statistic for small P µ and under ERMUTATION TESTS FOR COMPLETE BLOCK DESIGN Table 4
Power calculation under the exponential distribution P µ PowerF 0.056 0.081 0.209 0.424 0.668 0.861 0.982 1PowerLR 0.049 0.091 0.235 0.451 0.679 0.862 0.981 1PowerBN 0.053 0.093 0.238 0.455 0.681 0.866 0.982 1 the exponentially squared distribution, in Table 5, the Λ-statistic gives asubstantial increase in power compared to F -statistic for moderate values of P µ . In both cases there is no difference for higher powers. We note thatthe tests have essentially equal power up to computational accuracy underthe Normal, Uniform and Gamma (shape parameter 5) distributions. Theincrease in power becomes noticeable in long tail distributions like Expo-nential, Exponential Squared, Gamma (shape parameter 0.5) distributions.
6. Proofs of theorems of Section 3.
Proof of Theorem 1.
From (3), Ω = { x : κ ′ ( t ) = x, for some t ∈ R k − } ,the image of κ ′ ( · ). Using equation (1) we get the j th component of κ ′ ( t ), κ ′ j ( t ) = 1 b b X i =1 P π ∈ Π a iπ ( j ) exp( t T a iπ ) P π ∈ Π exp( t T a iπ ) . Here, a iπ ( j ) is the j th component of a iπ and¯ A = 1 b b X i =1 a i < b b X i =1 P π ∈ Π a iπ ( j ) exp( t T a iπ ) P π ∈ Π exp( t T a iπ ) < b b X i =1 a ik = ¯ A k . Using the same approach, we can conclude that for all distinct integers j , j , . . . , j k − taking values 1 , , . . . , k −
1, and for l = 1 , . . . , k − l X j =1 ¯ A j < b b X i =1 P π ∈ Π P lm =1 a iπ ( j m ) exp( t T a iπ ) P π ∈ Π exp( t T a iπ ) < k X j = k − l +1 ¯ A j . (7) Table 5
Power calculation under the exponentially squared distribution P µ PowerF 0.051 0.101 0.230 0.490 0.711 0.840 0.987 1PowerLR 0.057 0.169 0.319 0.545 0.727 0.832 0.976 1PowerBN 0.063 0.178 0.328 0.550 0.731 0.832 0.977 1 I. SAMONENKO AND J. ROBINSON
So Ω ⊂ P .Let us prove that Ω is a convex set. Let x, y ∈ Ω and c + d = 1, c , d > t ∈ R k − t T ( cx + dy ) − κ ( t ) = c ( t T x − κ ( t )) + d ( t T y − κ ( t )) ≤ c Λ( x ) + d Λ( y ) < ∞ . Since the expression t T ( cx + dy ) − κ ( t ) is bounded and convex, it has amaximum, so that cx + dy ∈ Ω andΛ( cx + dy ) ≤ c Λ( x ) + d Λ( y ) , so Ω is convex.To see that each vertex of the polytope is a limiting point of Ω, con-sider any vertex ¯ A ˆ π = ( ¯ A ˆ π (1) , ¯ A ˆ π (2) , . . . , ¯ A ˆ π ( k − ). Suppose ˆ π ( k ) = j and de-fine l , . . . , l k such that ˆ π ( l i ) = i , so l j = k . We can show thatlim t lj +1 →∞ · · · lim t lk →∞ lim t lj − →−∞ · · · lim t l →−∞ κ ′ ( t ) = ( ¯ A ˆ π (1) , . . . , ¯ A ˆ π ( k − ) . (8)To see this, write κ ′ ( t ) = 1 b b X i =1 P π ∈ Π a iπ exp( t T ( a iπ − a i ˆ π )) P π ∈ Π exp( t T ( a iπ − a i ˆ π ))(9) = 1 b b X i =1 a i ˆ π + P π ∈ Π ,π =ˆ π a iπ exp( t T ( a iπ − a i ˆ π ))1 + P π ∈ Π ,π =ˆ π exp( t T ( a iπ − a i ˆ π )) . Note that a i ˆ π ( l ) = a i is the smallest entry in the i th row so the coefficientof t l , a iπ ( l ) − a i ˆ π ( l ) = a iπ ( l ) − a i , is either positive or zero for any π ∈ Π.As t l → −∞ , only the permutations with π ( l ) = 1 give nonzero terms, solim t l →−∞ κ ′ ( t ) = 1 b b X i =1 a i ˆ π + P π ∈{ Π : π ( l )=1 } ,π =ˆ π a iπ exp( t T ( a iπ − a i ˆ π ))1 + P π ∈{ Π : π ( l )=1 } ,π =ˆ π exp( t T ( a iπ − a i ˆ π )) . (10)Continuing to take limits, in the order given in (8), removes all but the firstterm in the numerator and denominator of (10), to prove (8). Thus, thevertex ¯ A ˆ π is a limiting point of Ω. Since ˆ π is arbitrary, this holds for eachvertex of the polytope. Since Ω is a convex set enclosed by the edges of thepolytope, Ω is the interior of the polytope. (cid:3) Proof of Theorem 2.
Consider the expression t T x − κ ( t ). Using pre-vious notation set x = ¯ A ˆ π . Using the definition (1), we can write t T ¯ A ˆ π − κ ( t ) = − b b X i =1 log 1 k ! X π ∈ Π e t T ( a iπ − a i ˆ π ) , (11) ERMUTATION TESTS FOR COMPLETE BLOCK DESIGN since ¯ A ˆ π = b P bi =1 a i ˆ π . Then by the same argument used in Theorem 1, weget lim t lj +1 →∞ · · · lim t lk →∞ lim t lj − →−∞ · · · lim t l →−∞ [ t T ¯ A ˆ π − κ ( t )] = log k ! . From the definition of supremum and equation (2),Λ( x ) = sup t { t T ¯ A ˆ π − κ ( t ) } = log k ! . Since ˆ π is chosen arbitrarily Λ( x ) is equal to log k ! on all vertexes. Theseare the extreme points of Ω and Λ( x ) is convex, so Λ( x ) takes its maximumon the vertexes and is finite on all points of P and so on the boundary ofΩ. (cid:3) Proof of Theorem 3.
Using (4), choose x so that l X j =1 x s j = l X j =1 ¯ A j is true for some { s j } and l . The alternative choice will follow in the sameway. So x T t = l − X j =1 x s j ( t s j − t s l ) + l X j =1 ¯ A j t s l + k − X j = l +1 x s j t s j and, from (1), κ ( t ) can be written1 b b X i =1 log 1 k ! X π ∈ Π exp l − X j =1 a iπ ( j ) ( t s j − t s l ) + l X j =1 a iπ ( j ) t s l + k − X j = l +1 a iπ ( j ) t s j ! . Make the substitution u j = (cid:26) t s j − t s l , for 1 ≤ j < l,t s j , for l ≤ j ≤ k − , and use the first equality in (4), to write x T t − κ ( t ) as k − X j =1 ,j = l x s j u j − b b X i =1 log 1 k ! X π ∈ Π exp k − X j =1 ,j = l a iπ ( j ) u j + u l l X j =1 ( a iπ ( j ) − a ij ) ! , where for each 1 ≤ j ≤ l , P lj =1 ( a iπ ( j ) − a ij ) ≥
0. Let u l → −∞ and we havelim u l →−∞ ( x T t − κ ( t )) = k − X j =1 ,j = l x s j u j − b b X i =1 log 1 k ! X π ∈ ˆΠ exp k − X j =1 ,j = l a iπ ( j ) u j ! , I. SAMONENKO AND J. ROBINSON where ˆΠ = { π ∈ Π : P lj =1 ( a iπ ( j ) − a ij ) = 0 } . Let ˆΠ and ˆΠ be sets of allpermutations ˆ π and ˆ π of integers { , . . . , l } and { l + 1 , . . . , k } , respectively.Then the above expression can be rewrittenlim u l →−∞ ( x T t − κ ( t )) = l − X j =1 x s j u j − b b X i =1 log 1 l ! X ˆ π ∈ ˆΠ e P l − j =1 a i ˆ π j ) u j + k − X j = l +1 x s j u j − b b X i =1 log 1( k − l )! X ˆ π ∈ ˆΠ e P k − j = l +1 a i ˆ π j ) u j (12) + log k ! l !( k − l )! . Now taking suprema over u , . . . , u l − and u l +1 , . . . , u k − in the first twoterms on the right in (12), the statement of the theorem follows. (cid:3) REFERENCES
Borovkov, A. A. and
Rogozin, B. A. (1965). On the multi-dimensional central limittheorem.
Theory Probab. Appl. Brown, B. M. and
Maritz, J. S. (1982). Distribution-free methods in regression.
Austral.J. Statist. Fisher, R. A. (1935).
The Design of Experiments . Oliver and Boyd, Edinburgh.
Genz, A. (2003). Fully symmetric interpolatory rules for multiple integrals over hyper-spherical surfaces.
J. Comput. Appl. Math.
Jin, R. and
Robinson, J. (1999). Saddlepoint approximation near the endpoints of thesupport.
Statist. Probab. Lett. Kolassa, J. and
Robinson, J. (2011). Saddlepoint approximations for likelihood ra-tio like statistics with applications to permutation tests.
Ann. Statist. Robinson, J. (1982). Saddlepoint approximations for permutation tests and confidenceintervals.