a Department of Mathematics, Washington State University, Pullman, WA 99164-3113, USA;

A problem, which frequently arises in statistical analysis, is to compute the standardised q-variate central t-distribution function

Here , R is a q ´ q positive definite correlation matrix, and . There are many applications. One example is the area of multiple comparisons, where selection and ranking procedures or many-to-one comparisons require the evaluation of (see Hochberg and Tamhane (1987) for a comprehensive overview).

Until recently a direct evaluation of the multiple integral for arbitrary correlation matrices was considered computationally infeasible. Wang and Kennedy (1997) used interval analysis by applying the multivariate Taylor expansion to the density. But even low accuracy results required calculation times that were too large for practical purposes. Two different methodologies were developed by Somerville (1997, 1998). The first approach applies Monte Carlo techniques on a transformed multivariate t-integral. The second proposal uses a binning procedure combined with an univariate quadrature method to estimate the integral. Dunnett (1955) reduced to an inner q-variate normal integral combined with an outer chi integral. However, this method requires the evaluation of an additional integral. The advantage of this approach lies in the larger number of articles published about the multinormal integral (see, for example, Schervish, 1984, Solow, 1990, Genz, 1992, Joe, 1995, or Hajivassiliou et al., 1996). If further on R has a special structure (for example the product correlation structure r _i,j = l _{il j}, -1 < l _i, l _j < 1), the multivariate normal distribution reduces to a single integral over q multiplied univariate normal integrals. But since this approach is limited to problems with correlation matrices that have a special structure, we cannot use this approach for the more general problems considered in this paper. In this context Hsu (1992) proposed to approximate a given arbitrary correlation matrix R by the ‘closest’ , which still satisfies above structural condition. This can be done by use of factor analytic methods (Hsu, 1992) or by linear programming techniques (Hsu and Nelson, 1998). Finally, Edwards and Berry (1987) introduced crude Monte Carlo for equation (1) in order to calculate critical constants for different multiple comparison procedures.

Note that we could have defined by introducing the noncentral multivariate t-distribution. But we leave this generalisation for Section 5 and focus ourselves in the next sections at first on the alternate definition in central form without the noncentrality vector d . In the following section a new method for evaluating will be introduced. With this method, we transform the q-variate integrand into a product of univariate integrands. In Section 3, we describe algorithms for evaluating the transformed integral. An illustration of the use of this transformation by computing simultaneous confidence intervals is given in Section 4. In Section 5, we apply these methods to integrals required for power calculation of multiple contrast tests. We provide an expression in closed form and illustrate the numerical computation by an example.

We assume that in equation (1) for each i at least a_i or b_i is finite, because otherwise an appropriate transformation reduces the dimension of the problem by integrating the i^th variable explicitly. In a series of substitutions based on a talk held by Genz (1993) we transform to an integral over the (q - 1)-dimensional unit hypercube. The method is similar to the one used by Genz (1992) for the multivariate normal integral. In this form several numerical procedures can be used for the final evaluation, and some of these procedures are presented as algorithms in the next section.

Let CC^t = R be the Cholesky decomposition for R, where C is a lower triangular matrix. If we let x = Cy, then x^tR^-1x = y^tC^tC^-tC^-1Cy = y^ty and dx = |C|dy = |R| ^½dy and so

,

with and . We now split the integrand into a product of q factors by use of

.

When applied to T(a, b) this transformation yields

.

At this stage each of the y_i’s can be substituted using . Beginning with i = q the q-fold substitution results in

where follows from induction, and and are the new integration limits. For the last steps we set where denotes the univariate t-distribution with n degrees of freedom. Using equation (1) becomes

where and . With the sequence of transformations described here becomes now an integral of over the (q - 1)-dimensional unit hypercube, which can be evaluated with different multidimensional numerical computation methods.

In this section we present three numerical algorithms that use equation (2) to estimate for a given error requirement e . For the last algorithm, a lattice rule algorithm, the theoretical research is still ongoing. The other two algorithms are an acceptance-rejection sampling and a crude Monte Carlo algorithm. These methods are simpler and known to be reliable. Other Monte Carlo techniques could be applied but will not be analysed here. Deak (1990) gives a good overview with application to the multinormal case.

This procedure generates q-dimensional uniform random vectors w₁, ..., w_N and estimates T(a, b) using

,

with and , i = 1, ..., q, l = 1, ..., N. To control the simulated error we make use of the usual error estimate of the means

.

Further we denote by g the Monte Carlo confidence factor for the standard error. If, for example, g = 3, we then expect the actual error of to be less than the error bound e in 99.7% of the cases.

1. INPUT q, g , n , R, a, b, e .
2. Compute lower triangular Cholesky factor C for R.
3. Initialise N = 0, Intval = 0, Varsum = 0.
4. REPEAT

1. Generate uniform random w₁, ..., w_q Î [0, 1].
2. Set f = 1.
3. FOR i = 1, 2, ..., q

1. Set f = 0.
1. Set N = N + 1,

5. OUTPUT Intval, ErrEst, N.

The acceptance-rejection method is widely used and may be the most intuitive way to deal with equation (2). Deak (1990), however, showed that among the various Monte Carlo methods it is the one with the worst efficiency and therefore other approaches to evaluate T(a, b) are needed.

Monte Carlo Method

For estimating (2) by the crude Monte Carlo method we let w₁, ..., w_N be uniformly and independently distributed on [0, 1]^q-1. Then the random variables , l = 1, ..., N, are independent and their expected value is . Consequently the arithmetic average

is an unbiased estimator of the integral .

The Monte Carlo and Lattice rule algorithms that we describe in this section both require the evaluation of for particular values of w, so we give here an algorithm for that will be used by both of the numerical integration algorithms. In the following algorithm note that the initializations of and are required only for the first evaluation of (assuming a, b, , and C are fixed for a particular integral), and we would avoid wasteful computation of t-values be setting d_i = 0 or e_i = 1 if a_i = -¥ or b_i = ¥ , respectively.

1. Input a, b, w, , C.
2. Initialise f₁ = e₁ - d₁.
3. For i = 1, 2, ..., q - 1

In the following we give the algorithm when applying crude Monte Carlo on equation (2).

1. INPUT q, g , n , R, a, b, e .
2. Compute lower triangular Cholesky factor C for R.
3. Initialise N = 0, Intval = 0, Varsum = 0.
4. REPEAT

1. Generate uniform random w₁, ..., w_q-1 Î [0, 1].
2. Evaluate f_q =

1. Set N = N + 1,

5. OUTPUT Intval, ErrEst, N.

As seen from the central limit theorem, the Monte Carlo integration yields a probabilistic error bound in O(N^-½) (note that the order of magnitude does not depend on q). This means halving the error requires quadrupling the number of sample points. However, Quasi-Monte Carlo integration methods use sequences of nodes that are designed to be more uniform than random, while still using an integration formula with equal weights similar to (3). Niederreiter (1992) and others showed that under suitable conditions a deterministic error bound is given by O(N^-1log^q-1N). In the sequel we will make use of one approach out of this broad class, the so-called randomised lattice rule (Joe, 1990; Sloan and Joe, 1994, p.170)

Here N is the simulation size, usually being very small (e.g. 10 - 20), p corresponds to the fineness of the lattice and z Î Ñ ^q-1 denotes the strategically chosen lattice vector. Braces around vectors indicate that each component has to be replaced by its fractional part. Finally, w₁, ..., w_N denote again [0, 1]^q-1-uniform random vectors. The error estimate for is

.

For more details concerning lattice rules in general and the approach used here the reader is referred to the book of Sloan and Joe (1994).

The algorithm listed below consists basically of two loops. The l-loop stands for the outer sum in equation (4), and j is the variable of the inner sum. This is repeated several times until a pre-assigned error level e is reached. At each step the counting variable n is incremented by 1 and the number of lattice points p = p_n is increased in dependence on n.

1. INPUT q, g , n , R, a, b, e .
2. Compute lower triangular Cholesky factor C for R.
3. Initialise N = 10, n = 0.
4. REPEAT

1. Set n = n + 1, Intval = 0, Varsum = 0.
2. FOR l = 1, 2, ..., N

1. Set Latsum = 0.
2. Generate uniform random w₁, ..., w_q-1 Î [0, 1].
3. FOR j = 1, 2, ..., p_n

• Set .
• Evaluate f_q = .
• Set Latsum = Latsum + (f_q - Latsum)/j.

1. Set Varsum = Varsum + (l - 1)(Latsum - Intval)²/l,

1. Set ErrEst =

5. OUTPUT Intval, ErrEst, Np_n.

The best choice of z is still an open research question. Several proposals have been published in the literature, the most common of which is to choose z of the form , , proposed by Korobov (1960). This leaves us with the problem of how to choose h. The method we used minimises

with (x) = F(x)/4 + ¾ and F(x) = 1 + 2p ²(x² - x + 1/6), x Î [0, 1] (see Sloan and Joe, 1994, p. 173). A SAS/IML-program, which implements the lattice rule algorithm presented above together with an error estimation is available from the website with URL http://www.bioinf.uni-hannover.de/mcp_home/software.html. Fortran source codes for all three algorithms and an adaptive integration algorithm and an algorithm based on transforming (1) into a spherical coordinate system similar to the system used by Deak (1990) for multivariate normal integrals can be found the website with URL http://www.sci.wsu.edu/ math/faculty/genz.

4. Computation of Simultaneous confidence Intervals

In order to illustrate the use of the transformation method presented in Section 2 for evaluating multivariate t-probabilities, we consider one important application, the computation of simultaneous confidence intervals (see, for example, Edwards and Berry, 1987, or Hsu and Nelson, 1998). An example of the use of t-probability confidence intervals will be described in detail in the next section, where one-sided intervals are required. With these problems, we are given a desired confidence level and we need to determine an upper limit vector t = (t, ..., t) Î Ñ ^q, so that . We will illustrate the use of our algorithms with one of the problems to be considered in the next section (the Dunnett contrast case), which has q = 3, and a correlation matrix R given by

.

If we want to apply our algorithms on this problem, we need to combine our algorithms with a iterated-nonlinear-equation-solving algorithm. We let , so that we need to find t such that . We have successfully used various modified secant algorithms for this problem. The results given below were produced by the Pegasus method (see Ralston and Rabinowitz, 1978, p. 341). Using this algorithm we need a starting interval which contains our desired solution t, where and . At each stage an estimated t is computed, along with , and a new interval is produced. In order to simplify the start of the algorithm we used a = - 4 and b = 4, assuming and . The following Table I shows the behaviour of our Lattice rule algorithm () combined with the Pegasus method with termination if in each of the q components, for .

Þ insert TABLE I about here

The final t value was t = 2.1664, apparently correct to at least four digits, and this required a total of 22144 values.

A commonly used algorithm for the confidence interval problem is a type of rejection algorithm (see Edwards and Berry, 1987), which we will call %-rejection. The basic idea is to generate a large number, say m, of random vectors for given correlation matrix R and degrees of freedom n . For each let . Then sort the and let , where . The result is that after sorting we reject of the (smallest) and pick the smallest one left as an estimator for t. Several refinements to this basic algorithm have been suggested (see Hsu and Nelson, 1990, 1998, for discussion and further references), but we have found that rejection type algorithms are generally not very efficient for computing multivariate t-probabilities, and it is beyond the scope of this paper to provide a detailed comparison between our algorithms and %-rejection algorithms for this problem. However, we will describe how our algorithms could be used for %-rejection, too. The vectors needed for %-rejection are typically generated using a combination of multivariate normal and chi variables, but our transformations can also be used to generate vectors with correct distribution. The acceptance-rejection algorithm given in the last section can be easily modified to provide a simple %-rejection algorithm

1. INPUT m, a , n , R.
2. Compute lower triangular Cholesky factor C for R.
3. FOR l = 1, 2, ..., m

1. Generate uniform random w₁, ..., w_q Î [0, 1].
2. FOR i = 1, 2, ..., q

1. Set

4. Sort the .
1. OUTPUT

Using this algorithm with m = 22139, so that , we found . A standard 99% confidence interval obtained using a normal approximation to the binomial distribution is given by .

In the general many-to-one context of comparing several treatments with a control the following one-way layout model is widely used, consisting of one control group, for example a placebo or a negative control, and k treatments, labelled 0, 1, ..., k respectively. Let be the sample values, identically and independently N(µ_i, s ²) distributed. Denote by the sample mean and by the pooled variance estimator with degrees of freedom. For testing the following two alternatives are frequently considered. On the one hand the (one-sided) simple tree ordering is given by

H_A¹: $ i: µ₀ < µ_i, i = 1, ..., k.

But in practice it often occurs that several dose groups have to be compared with a negative control group, whereas a monotonous dose-response trend can be assumed a-priori. In this case, the total order restricted alternative

H_A²: µ₀ £ µ₁ £ ... £ µ_k, µ₀ < µ_k,

may be more appropriate. However, caution can be advisable, if one is not sure, whether this kind of underlying dose-response relationship really holds. Bauer (1997) showed that already small departures from this assumption may lead to inappropriate results if trend tests designed for such situations are used. These tests, such as the ones presented in the sequel, then are not useful in the sense that they do not control the probability of incorrectly declaring a dose to be effective when it is not.

We now consider the statistic

of multiple contrast tests according to Mukerjee et al. (1986, 1987), Westfall (1997) and Hothorn et al. (1997). This is just the maximum over q single contrasts tests

where denotes the central univariate t-distribution and c_i the contrast coefficients under the condition . Besides this, the choice of the c_i’s is free and numerous proposals concerning their optimal choice have been published (see, for example, Tamhane et al., 1996).

To test H₀ versus H_A¹ one can use, for example, the one-sided many-to-one test of Dunnett (1955), which can be formulated as a multiple contrast test with the weights

For the total restricted case H_A² Bretz and Hothorn (1998) developed a multiple contrast approach based on Williams’ test (1971), using in the balanced set-up the k contrasts

As each Tⁱ in formula (5) is univariate t-distributed, their joint distribution will follow per definition (Cornish, 1954; Dunn and Massey, 1965) a central q-variate t-distribution with n degrees of freedom and correlation matrix with entries of the form

The power for the multiple contrast test can now be derived as follows:

The probability term of equation (6) is q-variate noncentral t-distributed with the noncentrality parameter and n degrees of freedom. We denote the corresponding multivariate distribution function by , where a_l = -¥ " l. In equation (6) stands for the (1 - a )-equipercentage point of the central q-variate t-distribution, so that where a is the pre-specified type-I-error. We now make use of the relationship

.

Here U and are independent random variables distributed as standardised q-variate normal and chi with n degrees of freedom, . This leads to the explicit representation of the noncentral q-variate t-distribution function by

where F _q(.) is the q-variate normal integral with expectation 0, correlation matrix R and upper integration bound . Note that above formulae may be applied to both alternatives H_A¹ or H_A², whereas the first case for Dunnett’s many-to-one test was already considered by several authors, e.g. Bristol (1989) or Hayter and Liu (1992).

A correct evaluation of the probability term in equation (6) therefore consists of the following numerical problems. First we need the ability of computing the noncentral q-variate t-distribution function . This can be done by using the transformation indicated in equation (7). For the evaluation of the arising q-variate standard normal distribution F _q(.) several methods from the literature can be applied. Some of them are quoted in the introductory Section 1. The remaining single integral can be computed by using standard univariate integration techniques. Finally, we need to compute the critical value t in equation (6). The vector is the implicit solution of , where . Therefore, the computation reduces to the problem discussed in Section 4 and the Pegasus method (or the %-rejection algorithm) stated there can be applied. The evaluation of , on the other hand, requires the computation of , what can be done by using the algorithms provided in Section 3.

To illustrate the ideas discussed above, suppose k = 3 groups are to be compared with a negative control, for example increasing dose levels of a certain drug, where the control corresponds to the zero dose. Frequently in this situation it is possible to assume monotonicity a-priori from the experimenter point of view, i.e., if any effect exists, then it will increase with respect to the treatment levels. Therefore, the goal is to assess, whether a global trend exists. The interesting alternative is then chosen to be H_A². For this alternative, however, we do have the problem of defining appropriate contrast matrices and corresponding test statistics. While no general rule seems to be available yet, we compare some of the widely used single and multiple contrast tests throughout the literature for different dose-response relationships.

For experimental reasons we let the practical relevance shift be d = 1 and the size be set to a = 0.05. The variances of all 4 groups are assumed to be homogeneous and estimated from preceding pilot studies as s ² = 1. For simplicity we allocate the total sample size of 38 experimental units roughly according to the widely used -rule, resulting in the sample size vector (14, 8, 8, 8). The power, or, equivalently, the type-II-error, was calculated using a SAS/IML program based on the algorithms presented here. The results with an estimated error of 10^-4 are summarised in Table II.

Þ insert TABLE II about here

>From the computations done we can see that there is a strong dependency of the power of tests on the underlying dose-response shape. With increasing number of contrasts covering more possible dose-response shapes, the problem is reduced but not alleviated. Hence, it would be desirable to have some measure of goodness for a pre-chosen (multiple) contrast test. But even if Hayter and Liu (1992) succeeded in determining the least favourable configuration for Dunnett’s two-sided test in the partial balanced set-up, yielding the unweighted minimum power over all shapes, in the general case no solution exists. Several other aspects, such as the number of contrasts, sample size allocation or the ratio variance/relevance shift, influence the power function in complicated fashion as well. So, even if we managed to calculate the power of a multiple contrast in a specific situation, little can be assessed about their optimal choice.

TABLE I Iteration steps for the Pegasus method using the lattice rule algorithm .

TABLE II Power of several contrast tests for different dose-response profiles.