Introduction

We consider a general linear model with fixed effects:

$$Y = X\mbox{\boldmath $\beta$}+ \mbox{\boldmath $\epsilon$}.$$

We assume that we are given an $N\times 1$ data vector $Y$, $N\times p$ design matrix $X$, with unknown $p\times 1$ parameter vector $\mbox{\boldmath$\beta$}$, and an $N\times 1$ error vector $\mbox{\boldmath$\epsilon$}$, with i.i.d. normally distributed components with unknown variance $\sigma^2$. The setup for multiple comparison problems (see Hochberg and Tamhane, 1987 and Hsu, 1992, 1996, Somerville, 1997, 1998, Stoline, 1981) provides an $m \times p$ comparison (or contrast) matrix $C$. The covariance matrix for the multiple comparison problem is then $\mbox{\boldmath$\Sigma$}= C(X^t X)^{-1}C^t$, an $m \times m$ positive semi-definite matrix. The resulting basic numerical problem that is the focus for this paper is the determination of confidence intervals (CI's) for

$$x_i = \sum_{j=1}^p c_{i,j}\beta_j, \ \ i=1,\ldots,m .$$

The distribution for ${\bf x}$ is an $m$-variate Student's t (MVT), with covariance matrix $\mbox{\boldmath$\Sigma$}$ and degrees of freedom $\nu = N - rank(X)$. We use the Dunnett (1954) definition of the MVT distribution given by

$${\bf T}({\bf a}, {\bf b}, \mbox{\boldmath $\Sigma$}, \nu) = ... ...u}},\frac{s{\bf b}}{\sqrt{\nu}},\mbox{\boldmath $\Sigma$}) ds,$$

where the multivariate normal distribution function

$$\mbox{\boldmath $\Phi$}({\bf a},{\bf b}, \mbox{\boldmath $\S... ...2} {\bf x}^t \mbox{\boldmath $\Sigma$}^{-1} {\bf x}} d{\bf x},$$

${\bf x}= (x_1, x_2, ..., x_m)^t$, $-\infty \leq a_i < b_i \leq \infty$ for all $i$, and $\mbox{\boldmath$\Sigma$}$ is a positive semi-definite symmetric $m \times m$ matrix. Efficient and robust numerical software is now available for MVT distribution computations for $1 \leq m \leq 20$ (see Genz and Bretz, 1999, 2000, also Genz and Kwong, 1999, for how to handle singular $\mbox{\boldmath$\Sigma$}$'s). The final input for the CI determination problem is a confidence level $\alpha$. The actual numerical problem consists of finding the critical value $t_\alpha$ where $P(t_\alpha) = 1- \alpha$, with

• $P(t) = {\bf T}(-{\bf\infty},{\bf t}, \nu, \mbox{\boldmath$\Sigma$})$ for one-sided CI's, or
• $P(t) = {\bf T}(\ \ -{\bf t},{\bf t}, \nu, \mbox{\boldmath$\Sigma$})$ for two-sided CI's.

Here ${\bf t}= (t,...,t)^t$ and ${\mbox{\boldmath$\infty$}} = (\infty, ..., \infty)^t$. The numerical problem is therefore a combined problem of using an appropriate numerical optimization method to determine $t_\alpha$ with an efficient numerical integration method for evaluating $P(t)$.