Computing Exact Max-t Critical Values

From a Max-t Statistic to a Probability Calculation

A max-t procedure protects a family of contrasts by comparing the largest absolute standardized contrast with one common critical value. For correlated contrasts, the critical value depends on their joint distribution: treating the statistics as independent can give the wrong familywise coverage. Computing the value therefore means evaluating a joint probability for a specified correlation matrix and degrees of freedom, then solving for the value that gives the desired coverage.

Suppose there are \(m\) prespecified contrasts and their standardized estimation errors have a central multivariate t distribution with degrees of freedom \(\nu\) and correlation matrix \(R\). Here \(R\) is positive definite and has ones on its diagonal. Write \(T=(T_1,\ldots,T_m)^\mathsf{T}\), and define the two-sided maximum \(M=\max_{1\leq j\leq m}|T_j|\).

Definition: For a desired familywise error rate \(\alpha\), the two-sided max-t critical value \(c_\alpha\) is the value satisfying \(P(M\leq c_\alpha)=1-\alpha\). Equivalently, \(c_\alpha\) is the \(1-\alpha\) quantile of \(M\). The simultaneous intervals for the contrasts use this same multiplier on their standard errors.

The phrase exact max-t critical value refers to calibration under the specified finite-sample multivariate t distribution, rather than to a Bonferroni bound or a large-sample normal approximation. The probability usually has no simple closed form. It can, however, be evaluated by numerical integration and then inverted to find \(c_\alpha\).

The Multivariate-t Mixture Representation

A useful route to the probability is the scale-mixture representation of the central multivariate t distribution. Let \(Z\sim N_m(0,R)\) and \(U\sim\chi^2_\nu\) be independent. Then \(T=Z/\sqrt{U/\nu}\) has the central multivariate t distribution with correlation matrix \(R\) and degrees of freedom \(\nu\). The same random scale \(U\) divides every component; this shared denominator matters even when \(R\) is the identity matrix.

Theorem: If \(R\) is positive definite and \(U,Z\) have the distributions and independence just stated, then for \(c\geq0\), \(P(M\leq c)=\int_0^\infty G_R(c\sqrt{u/\nu})f_{\chi^2_\nu}(u)\,du\), where \(G_R(a)=P(-a\leq Z_j\leq a\text{ for every }j)\) is a centered multivariate-normal rectangle probability.

Proof. Conditional on \(U=u>0\), the event \(M\leq c\) is equivalent to \(|Z_j|\leq c\sqrt{u/\nu}\) for every component \(j\). Its conditional probability is therefore \(G_R(c\sqrt{u/\nu})\). Averaging that conditional probability over the chi-square density of \(U\), by the law of total probability, gives the displayed integral. \(\square\)

The integral has a clear computational interpretation: for each possible scale \(u\), evaluate a normal probability over a box whose half-width is \(c\sqrt{u/\nu}\), then average over the distribution of the scale. One can integrate directly over \(u\), or transform using the chi-square quantile function \(q_\nu(v)\):

P(M\leq c)=\int_0^1 G_R\!\left(c\sqrt{\frac{q_\nu(v)}{\nu}}\right)\,dv.

The second form places the integration on the fixed interval \((0,1)\). In either form, the multivariate-normal rectangle probability is itself a numerical task. For small \(m\), adaptive quadrature may be practical; for larger \(m\), algorithms for multivariate-normal probabilities commonly use transformations and randomized quasi-Monte Carlo integration.

Worked Examples

Worked Example: Constructing the Correlation Matrix for Three Contrasts

Let a four-dimensional observation have covariance matrix \(\Sigma=\operatorname{diag}(1,4,9,16)\). Consider the three contrast vectors \(a_1=(1,-1,0,0)^\mathsf{T}\), \(a_2=(0,0,1,-1)^\mathsf{T}\), and \(a_3=(1,0,-1,0)^\mathsf{T}\). Find the correlation matrix needed for max-t calibration.

For contrasts \(a_i^\mathsf{T}\mu\), their covariance entries are \(a_i^\mathsf{T}\Sigma a_j\). The variances are \(a_1^\mathsf{T}\Sigma a_1=1+4=5\), \(a_2^\mathsf{T}\Sigma a_2=9+16=25\), and \(a_3^\mathsf{T}\Sigma a_3=1+9=10\).

The covariances are \(a_1^\mathsf{T}\Sigma a_2=0\), \(a_1^\mathsf{T}\Sigma a_3=1\), and \(a_2^\mathsf{T}\Sigma a_3=-9\). Dividing each covariance by the product of the corresponding standard deviations gives

R= \begin{pmatrix} 1&0&1/\sqrt{50}\\ 0&1&-9/\sqrt{250}\\ 1/\sqrt{50}&-9/\sqrt{250}&1 \end{pmatrix} \approx \begin{pmatrix} 1&0&0.1414\\ 0&1&-0.5692\\ 0.1414&-0.5692&1 \end{pmatrix}.

For example, the \((2,3)\) entry is \(-9/(\sqrt{25}\sqrt{10})=-9/\sqrt{250}\). This matrix, together with the relevant degrees of freedom and familywise error rate, specifies the probability calculation. It is important to use the correlation matrix of the standardized contrasts, not the unscaled covariance matrix of their estimates.

Worked Example: Reducing a Finite-df Calculation to One Integral

Suppose two standardized contrasts have degrees of freedom \(\nu=8\) and correlation matrix \(R=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\). To find a two-sided 95% max-t critical value, set the rectangle probability equal to \(0.95\). For a bivariate normal vector with correlation \(\rho=0.5\), condition on its first coordinate to obtain

G_R(a)=\int_{-a}^{a}\phi(z) \left[ \Phi\!\left(\frac{a-0.5z}{\sqrt{0.75}}\right) -\Phi\!\left(\frac{-a-0.5z}{\sqrt{0.75}}\right) \right]\,dz.

Here \(\phi\) and \(\Phi\) are the standard normal density and distribution function; the conditional variance is \(1-0.5^2=0.75\). The required critical value is the root of

\int_0^\infty G_R\!\left(c\sqrt{u/8}\right)f_{\chi^2_8}(u)\,du=0.95.

A numerical calculation proceeds in two nested stages: evaluate the one-dimensional integral for \(G_R(a)\), then integrate over \(u\) and find the root in \(c\). For example, bisection can start with \(c=0\), where the probability is zero, and increase an upper endpoint until the computed probability exceeds \(0.95\). Each midpoint divides the bracket according to whether its probability is below or above \(0.95\). The final reported value should be accompanied by the numerical tolerances used for the inner and outer integrations and for the root search.

Worked Example: An Independent-Normal Benchmark

As a limiting check, consider \(m=2\) independent standard normal statistics. This is the \(\nu\to\infty\) benchmark with \(R=I_2\), not a finite-degrees-of-freedom t calculation. Independence gives \(P(M\leq c)=[2\Phi(c)-1]^2\). Setting this equal to \(0.95\) yields

2\Phi(c)-1=\sqrt{0.95},\qquad \Phi(c)=\frac{1+\sqrt{0.95}}{2}=0.9873397, \qquad c\approx 2.237.

The calculation checks both the two-sided probability and the familywise target: the probability that both absolute values are at most \(c\) is \(0.95\). It is a useful limiting benchmark for software, but it must not be substituted for the finite-df calculation when the statistics have a shared t scale.

Root Finding and Numerical Error Checks

For a positive-definite \(R\), the function \(H(c)=P(M\leq c)\) is continuous and strictly increasing from zero to one as \(c\) ranges from zero to infinity. Continuity follows because \(M\) has no point masses: conditional on \(U=u>0\), the normal vector has a density, and the boundary of the rectangle has probability zero. Strict increase follows because for any \(0\leq c_1<c_2\), the region where \(c_1<M<c_2\) contains an open set of positive volume; the conditional normal density is positive there. Averaging over \(U\) preserves positive probability. Thus the equation \(H(c)=1-\alpha\) has a unique solution for \(0<\alpha<1\).

A practical workflow is to evaluate \(H(c)\) at trial values, bracket the target, and use bisection or a safeguarded root solver. Bisection is especially transparent: once \(H(c_L)<1-\alpha<H(c_U)\), retain the half-interval whose endpoints still straddle the target. Continue until the bracket is narrower than the desired critical-value tolerance. A solver should not be allowed to hide an inaccurate probability calculation behind a small change in its proposed root.

Check the result in several ways. Confirm that \(R\) is symmetric with unit diagonal and positive definite; verify that probabilities increase with \(c\); and evaluate the final probability at the reported value, as well as at nearby values on either side. Tighten integration tolerances and confirm that the resulting root is stable. If possible, compare two numerical integration settings or methods. Small integration error matters most when the probability curve is shallow near the target, because a small probability error can then produce a noticeable critical-value error.

Worked Example: Interpreting a Simulation-Based Check

For the two-contrast, \(\nu=8\), \(\rho=0.5\) model above, suppose an independent simulation at a proposed value \(c_0\) generates \(100{,}000\) draws from the multivariate-t representation. Assume \(95{,}012\) draws satisfy \(\max_j|T_j|\leq c_0\). The estimated coverage is

\widehat H(c_0)=\frac{95012}{100000}=0.95012.

Treating the indicators as independent Bernoulli observations, the estimated standard error is \(\sqrt{0.95012(1-0.95012)/100000}\approx0.000688\). A rough 95% Monte Carlo interval is therefore \(0.95012\pm1.96(0.000688)\), or approximately \([0.94877,0.95147]\). This check is consistent with a target of \(0.95\), but it is not precise enough by itself to establish a highly accurate critical value. More simulation draws narrow the Monte Carlo uncertainty; deterministic or quasi-Monte Carlo integration can provide a more efficient primary calculation.

Simulation checks should use the joint model, including the common chi-square scale. At each replicate, draw one \(U\), draw one \(Z\) with correlation matrix \(R\), and form every component of \(T=Z/\sqrt{U/\nu}\) using that same \(U\). Drawing separate denominators for different contrasts would check a different distribution.

Why the Computation Matters

The max-t multiplier translates a joint coverage requirement into a common interval width. A critical value that is too small loses familywise coverage; one that is unnecessarily large makes all intervals wider. Numerical accuracy is therefore part of the inferential procedure, not merely a software detail. Reporting the degrees of freedom, the contrast-correlation matrix, the target probability, and the integration and root-finding tolerances makes the calibration reproducible.

Common pitfall: “The correlations are zero, so the t statistics are independent; I can multiply their marginal t probabilities.” Zero correlation does not imply independence for a multivariate t vector with finite degrees of freedom. In its mixture representation, the components share \(U\), so their magnitudes are linked even when the components of \(Z\) are independent. The correction is to evaluate the joint multivariate-t probability—or its chi-square mixture of joint normal rectangle probabilities—using the specified \(R\) and \(\nu\). Only in a setting that genuinely supplies independence, such as independent normal coordinates in the infinite-df limit, is the product calculation justified.

Key takeaway: Compute an exact max-t critical value by evaluating the multivariate-t rectangle probability for the specified contrast-correlation matrix and degrees of freedom, then solving \(H(c)=1-\alpha\). Validate both the probability calculation and the numerical root, and use simulation as an error check rather than as an unexamined substitute for the joint model.

Check Your Understanding

Use these questions to check the logic of the computation.

Why does the multivariate-t rectangle probability average a normal rectangle probability over a chi-square scale?
What additional information, beyond the number of contrasts and the target error rate, determines a finite-df max-t critical value?
Why can zero off-diagonal entries in \(R\) fail to make finite-df t components independent?
How can bisection verify that a numerical root brackets the desired familywise coverage?
What should be held common across contrast components in a simulation from the multivariate-t representation?