The Exact Null Law of Roy's Largest Root

Why the Joint Root Law Is Needed

Roy’s largest root records the strongest hypothesis direction, so its null probability depends on whether every generalized root stays below a threshold. Knowing the distribution of one root at a time is not enough: under the central MANOVA null, the roots are dependent. The exact null CDF is obtained by integrating their joint density over the region where all of them are at most the threshold.

The integration is especially useful for understanding what an exact calibration entails. In general, it produces a multivariate integral rather than a familiar one-dimensional reference distribution. Some dimensions simplify to closed forms; in others, the displayed integral itself is the basis for numerical evaluation.

Central-Null Roots and Their Joint Density

Let $p$ be the response dimension, $h$ the hypothesis degrees of freedom, and $\nu$ the error degrees of freedom. Under the classical central MANOVA null, take independent matrices $H\sim W_p(h,\Sigma)$ and $E\sim W_p(\nu,\Sigma)$, where $\Sigma$ is positive definite, and assume $\nu>p-1$. The generalized roots are the nonzero eigenvalues of $E^{-1}H$. There are $m=\min(p,h)$ positive roots with probability one. Write $R$ for their maximum, Roy’s largest root.

Definition: Set $m=\min(p,h)$, $n=\max(p,h)$, $a=n/2$, and $b=(h+\nu-n)/2$. The central-null joint density of the $m$ roots $\lambda_1,\ldots,\lambda_m$, treated as unordered coordinates on $(0,\infty)^m$, is $$ f(\lambda_1,\ldots,\lambda_m) = C_{m,a,b} \prod_{i=1}^{m}\lambda_i^{\,a-(m+1)/2}(1+\lambda_i)^{-(a+b)} \prod_{i<j}|\lambda_i-\lambda_j|, $$ where $$ C_{m,a,b} = \frac{\pi^{m^2/2}\Gamma_m(a+b)} {m!\,\Gamma_m(a)\Gamma_m(b)\Gamma_m(m/2)}, \qquad \Gamma_m(c)=\pi^{m(m-1)/4}\prod_{j=1}^{m}\Gamma\!\left(c-\frac{j-1}{2}\right). $$ Here $a>(m-1)/2$ and $b>(m-1)/2$, as ensured by the stated degrees-of-freedom conditions.

This is the central MANOVA root density, with the normalization for unordered roots. The absolute differences in the last product are the eigenvalue interaction term: the joint density does not factor into $m$ separate densities. Also, $a+b=(h+\nu)/2$, so the denominator exponent is determined by the hypothesis and error degrees of freedom together. When $h<p$, only $h$ roots are positive; the remaining eigenvalues of $E^{-1}H$ are zero and are not included in this density.

The Exact CDF as a Bounded-Region Integral

Theorem: Under the central-null assumptions above, for $t\leq 0$, $P(R\leq t)=0$. For $t>0$, $$ P(R\leq t) = C_{m,a,b} \int_{[0,t]^m} \prod_{i=1}^{m}x_i^{\,a-(m+1)/2}(1+x_i)^{-(a+b)} \prod_{i<j}|x_i-x_j| \,dx_1\cdots dx_m. $$

Proof. Since $R=\max_i\lambda_i$, the event $R\leq t$ is exactly the event $\lambda_i\leq t$ for every $i$. All roots are positive with probability one, so when $t>0$ this event is the region $[0,t]^m$ in root coordinates. Integrating the joint density over that region gives the stated expression. If $t\leq0$, the event has probability zero because all roots are positive. $\square$

This formula is an exact CDF, not an approximation. Its upper limit is the same in each coordinate because the event concerns the maximum. The integral over $[0,\infty)^m$ is one by the joint-density normalization, so the CDF tends to one as $t$ tends to infinity.

A Beta-Type Coordinate Transformation

For $t>0$, set $z_i=\lambda_i/(1+\lambda_i)$. This maps each positive root to the interval $(0,1)$, and the event $\lambda_i\leq t$ becomes $z_i\leq s$, where $s=t/(1+t)$. In these coordinates the CDF has a particularly transparent integrand.

Result: With $s=t/(1+t)$, the exact central-null CDF is $$ P(R\leq t) = C_{m,a,b} \int_{[0,s]^m} \prod_{i=1}^{m}z_i^{\,a-(m+1)/2}(1-z_i)^{\,b-(m+1)/2} \prod_{i<j}|z_i-z_j| \,dz_1\cdots dz_m. $$

Derivation. The inverse change of variables is $\lambda_i=z_i/(1-z_i)$, with $d\lambda_i=dz_i/(1-z_i)^2$. For each pair, $$ |\lambda_i-\lambda_j| = \frac{|z_i-z_j|}{(1-z_i)(1-z_j)}. $$ Thus the product of pairwise differences contributes a factor $(1-z_i)^{-(m-1)}$ for each $i$. The power of $1-z_i$ from the one-coordinate terms, the pairwise differences, and the Jacobian is $$ -\left(a-\frac{m+1}{2}\right) -(m-1) +(a+b)-2 = b-\frac{m+1}{2}. $$ The power of $z_i$ remains $a-(m+1)/2$, which gives the displayed density and CDF. $\square$

The transformed integral is over a bounded part of the unit cube, with beta-type powers at its boundaries. It also makes threshold calculations convenient: once $s$ is computed, the upper limit is simply $s$ in every coordinate.

Worked Example: An Exact CDF for Two Roots

Take $p=2$, $h=3$, and $\nu=3$. Then $m=2$, $n=3$, $a=b=3/2$, and $a+b=3$. The exponents in the transformed density are both zero. The normalizing constant is $$ C_{2,3/2,3/2} = \frac{\pi^2\Gamma_2(3)} {2\,\Gamma_2(3/2)^2\Gamma_2(1)} =3, $$ because $\Gamma_2(3/2)=\pi/2$, $\Gamma_2(3)=3\pi/2$, and $\Gamma_2(1)=\pi$.

Consequently, the transformed joint density is $3|z_1-z_2|$ on the unit square. For $s=t/(1+t)$, integrate over the square $[0,s]^2$. Scaling both coordinates by $s$ shows that $$ P(R\leq t) = 3\int_{[0,s]^2}|z_1-z_2|\,dz_1\,dz_2 = 3s^3\int_{[0,1]^2}|u-v|\,du\,dv. $$ The last integral is $2\int_0^1\int_0^u(u-v)\,dv\,du=2\int_0^1u^2/2\,du=1/3$. Hence $P(R\leq t)=s^3$. At $t=1$, $s=1/2$, so $P(R\leq1)=1/8$. At $t\to\infty$, $s\to1$, and the CDF tends to one, as required.

Worked Example: A Different Root Interaction

Let $p=3$, $h=2$, and $\nu=3$. Now $m=2$, $n=3$, $a=3/2$, and $b=1$. The constant is $$ C_{2,3/2,1} = \frac{\pi^2\Gamma_2(5/2)} {2\,\Gamma_2(3/2)\Gamma_2(1)\Gamma_2(1)} = \frac34, $$ using $\Gamma_2(5/2)=3\pi/4$, $\Gamma_2(3/2)=\pi/2$, and $\Gamma_2(1)=\pi$. The transformed density is therefore $$ \frac34\,|z_1-z_2|(1-z_1)^{-1/2}(1-z_2)^{-1/2}. $$

For a threshold $t=1$, the transformed upper limit is $s=1/2$. To evaluate the integral, put $y=1-\sqrt{1-z}$, so $z=2y-y^2$ and $(1-z)^{-1/2}dz=2\,dy$. The threshold $z\leq1/2$ corresponds to $0\leq y\leq c$, where $c=1-1/\sqrt{2}$. Also, $$ |z_1-z_2|=|y_1-y_2|(2-y_1-y_2). $$ The CDF becomes $$ P(R\leq1) = 3\int_{[0,c]^2}|y_1-y_2|(2-y_1-y_2)\,dy_1\,dy_2 = 6\int_0^c\int_0^y(y-v)(2-y-v)\,dv\,dy. $$ The inner integral is $y^2-2y^3/3$. Therefore the result is $$ 6\int_0^c\left(y^2-\frac{2y^3}{3}\right)dy = 2c^3-c^4 \approx 0.0429. $$ As a normalization check, replacing $c$ by 1 gives $2-1=1$, the full probability.

When Exact Evaluation Still Requires Integration

Worked Example: Setting Up a Higher-Dimensional Calculation

Take $p=4$, $h=2$, and $\nu=6$. Here $m=2$, $n=4$, $a=b=2$, and $a+b=4$. From $\Gamma_2(2)=\pi/2$, $\Gamma_2(4)=45\pi/4$, and $\Gamma_2(1)=\pi$, the constant is $$ C_{2,2,2} = \frac{\pi^2\Gamma_2(4)} {2\,\Gamma_2(2)^2\Gamma_2(1)} = \frac{45}{2}. $$ The powers in the transformed density are $1/2$ for both $z_i$ and $1-z_i$. Thus at $t=1$, for which $s=1/2$, $$ P(R\leq1) = \frac{45}{2} \int_{[0,1/2]^2}|z_1-z_2| \sqrt{z_1z_2(1-z_1)(1-z_2)} \,dz_1\,dz_2. $$

Using symmetry across the diagonal, this can be written as the one-dimensional outer integral $$ 45\int_0^{1/2}\sqrt{x(1-x)} \left[\int_0^x(x-y)\sqrt{y(1-y)}\,dy\right]dx. $$ The inner integral is finite and the expression is an exact probability; it can be evaluated by quadrature when a numerical CDF value is needed. This example also shows why an exact null law need not reduce to a simple named univariate distribution.

Common Pitfall and Practical Meaning

Common mistake: Treat the roots as independent, find the marginal CDF of one root, and raise it to the $m$th power to obtain the CDF of their maximum. That would require the event probabilities for the roots to factor, which they do not: the joint density contains $\prod_{i<j}|\lambda_i-\lambda_j|$. The correction is to integrate the joint density over $[0,t]^m$, or over $[0,t/(1+t)]^m$ after transformation.

The exact integral explains both the strength and the limitation of Roy’s criterion. Its null probability targets the event that no root exceeds the threshold, directly matching the largest-root statistic. But outside tractable cases, evaluation requires multidimensional integration. The integral is the exact central-null reference probability; it should not be confused with a later approximation or with a product of marginal probabilities.

Check Your Understanding

Use the definitions and integration arguments above to answer these questions.

Why does the event $R\leq t$ correspond to the full cube $[0,t]^m$, rather than a region in which only one root is constrained?
For $t>0$, what is the upper limit in each transformed coordinate $z_i$, and how is it obtained?
In the case $p=2$, $h=3$, and $\nu=3$, calculate $P(R\leq2)$ using the exact CDF derived in the example.
Which part of the joint density shows that the roots cannot be treated as independent?
For $p=4$, $h=2$, and $\nu=6$, state the exponents of $z_i$ and $1-z_i$ in the transformed joint density.