Hotelling’s T² test
Section of Estimation and inference for the multivariate normal, part of Statistical learning
Hotelling’s T² test
| Type | Section |
| Taught in | Estimation and inference for the multivariate normal |
| Part | Statistical learning |
A sample mean, however accurately computed, is not a test. A test needs a null value, an alternative, a statistic, and the distribution of that statistic when the null is true. The statistic derived here is the likelihood ratio for a normal mean vector. Its null distribution is read from the normal-and-Wishart pair.
The univariate t test as a template
Main article: Univariate t test as a template for T²
The shape of the result is already visible when p=1. For independent \operatorname{N}(\mu,\sigma^2) data, the test of H_0:\mu=\mu_0 against H_a:\mu\neq\mu_0 uses
t=\frac{\bar x-\mu_0}{s/\sqrt{n}}\sim t_{n-1}
under H_0, where t_{n-1} is the Student t distribution. Under H_0 the distribution of t does not depend on the unknown (\mu,\sigma^2). A quantity with that property is a pivot, so the critical value can be taken from the t_{n-1} table alone. At level \alpha the two-sided rule rejects H_0 when |t|>t_{\alpha/2;\,n-1}. Throughout this page the first subscript of a critical value is the upper-tail probability, and the remaining subscripts are degrees of freedom. The same point is the quantile of order 1-\alpha in the notation of Statistical inference. The algebra below is arranged so that the multivariate statistic reduces to the square of this t when p=1.
Hypotheses and the likelihood ratio
Main article: Hotelling’s T² as a likelihood-ratio test
Take n independent observations from \operatorname{N}_p(\boldsymbol{\mu},\Sigma), with both parameters unknown, and test
H_0:\boldsymbol{\mu}=\boldsymbol{\mu}_0 \qquad\text{against}\qquad H_1:\boldsymbol{\mu}\neq\boldsymbol{\mu}_0.
The unrestricted fit uses \bar{\mathbf{x}} and A/n. Under the null, the mean is fixed and only the covariance can adapt. The likelihood-ratio statistic compares those two best fits. The construction does not by itself establish uniform optimality against every multivariate alternative. Compare the maximised likelihood with \boldsymbol{\mu} fixed at \boldsymbol{\mu}_0 to the maximised likelihood with \boldsymbol{\mu} free:
\lambda=\frac{\max_\Sigma L(\boldsymbol{\mu}_0,\Sigma)}{\max_{\boldsymbol{\mu},\Sigma}L(\boldsymbol{\mu},\Sigma)}\in[0,1].
The numerator maximises over a subset of the denominator’s parameter space, so the ratio cannot exceed one. If it is small, imposing the proposed mean costs a great deal of fit.
Maximisation under H_0. Under H_0 only \Sigma is free. The derivative that produced \widehat\Sigma=A/n measures residuals from the centring point in use. With residuals measured from \boldsymbol{\mu}_0 rather than from a free mean, the same derivative gives
\widehat\Sigma_\omega=\frac1n\sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu}_0)(\mathbf{x}_i-\boldsymbol{\mu}_0)^\top.
Split those residuals at \bar{\mathbf{x}}:
\begin{aligned} \widehat\Sigma_\omega &=\frac1n\sum_{i=1}^n\bigl[(\mathbf{x}_i-\bar{\mathbf{x}})+(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)\bigr] \bigl[(\mathbf{x}_i-\bar{\mathbf{x}})+(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)\bigr]^\top\\ &=\frac1n\Bigl[A+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top\Bigr]. \end{aligned}
The cross terms vanish because \sum_i(\mathbf{x}_i-\bar{\mathbf{x}})=\mathbf{0}. The constrained covariance estimate is A/n plus the outer product of \bar{\mathbf{x}}-\boldsymbol{\mu}_0. The added term is positive semi-definite, and it is zero only when the sample mean equals the null mean.
Maximisation without restriction. With both parameters free, \widehat{\boldsymbol{\mu}}_\Omega=\bar{\mathbf{x}} and \widehat\Sigma_\Omega=A/n. The denominator of \lambda uses those two estimators.
The ratio. Substitute the two maximisers. The powers of 2\pi agree. The exponents also agree, for a reason that separates from the determinant algebra: at each maximiser the quadratic form in the exponent equals np. A scalar equals its own trace, and the trace is cyclic, so \mathbf{v}^\top M^{-1}\mathbf{v}=\operatorname{tr}(M^{-1}\mathbf{v}\mathbf{v}^\top). Apply that identity to the constrained exponent:
\begin{aligned} \sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu}_0)^\top\widehat\Sigma_\omega^{-1}(\mathbf{x}_i-\boldsymbol{\mu}_0) &=\operatorname{tr}\Bigl(\widehat\Sigma_\omega^{-1}\sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu}_0)(\mathbf{x}_i-\boldsymbol{\mu}_0)^\top\Bigr)\\ &=\operatorname{tr}\bigl(\widehat\Sigma_\omega^{-1}\cdot n\widehat\Sigma_\omega\bigr) =n\operatorname{tr}(I_p)=np. \end{aligned}
At the constrained maximiser, n\widehat\Sigma_\omega equals \sum_i(\mathbf{x}_i-\boldsymbol{\mu}_0)(\mathbf{x}_i-\boldsymbol{\mu}_0)^\top. The same steps with \widehat\Sigma_\Omega and residuals from \bar{\mathbf{x}} also give np, because n\widehat\Sigma_\Omega=A. Both exponential factors equal e^{-np/2} and cancel. Normality was used to write the likelihood. The cancellation itself is this trace identity at the two critical points. The likelihood ratio is a ratio of determinants,
\lambda=\frac{|\widehat\Sigma_\Omega|^{n/2}}{|\widehat\Sigma_\omega|^{n/2}} =\frac{|A|^{n/2}}{\bigl|A+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top\bigr|^{n/2}}.
A determinant identity. For conformable blocks with B invertible,
\begin{vmatrix}B&C\\ D&E\end{vmatrix} =\begin{vmatrix}B&C\\ 0&E-DB^{-1}C\end{vmatrix} =|B|\,|E-DB^{-1}C|.
The first equality subtracts DB^{-1} times the first block row from the second, which does not change the determinant. The same determinant, eliminated on the lower-right block instead, has the companion form |B|\,|E-DB^{-1}C|=|E|\,|B-CE^{-1}D|. Apply the identity with B=A, C=-\sqrt{n}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0), D=\sqrt{n}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top and E=1:
\begin{vmatrix} A & -\sqrt{n}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)\\ \sqrt{n}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top & 1 \end{vmatrix} =|A|\Bigl(1+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top A^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)\Bigr).
With E=1 the companion evaluation is |A+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top|, because -CE^{-1}D=+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top. The two evaluations therefore give
\bigl|A+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top\bigr| =|A|\Bigl(1+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top A^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)\Bigr).
This comparison is algebraic. Raise \lambda to the power 2/n:
\lambda^{2/n} =\frac{1}{1+n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top A^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)} =\frac{1}{1+t^2/(n-1)},
where A=(n-1)S has been used to name the quadratic form
t^2=n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top S^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0) =(n-1)\,n(\bar{\mathbf{x}}-\boldsymbol{\mu}_0)^\top A^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}_0).
The two expressions agree because S^{-1}=(n-1)A^{-1}. Thus \lambda is a strictly decreasing function of t^2. Rejecting for small \lambda is the same test as rejecting for large t^2. The block determinant converts the ratio of covariance determinants into that quadratic form, and the conversion does not use normality. Normality enters when the threshold is given a probability. The statistic t^2 is n times the squared Mahalanobis distance from \bar{\mathbf{x}} to \boldsymbol{\mu}_0, computed with S. When p=1, S reduces to s^2 and t^2 is the square of the univariate t statistic.
The null distribution
Main article: Null distribution of Hotelling’s
Write T^2 for the random version of t^2, so that \overline{\mathbf{X}}, S and A are random. The event \lambda<\lambda_0 is the event T^2>t_0^2, where t_0^2=(n-1)(\lambda_0^{-2/n}-1). The threshold is chosen so that the rejection probability equals \alpha when H_0 is true. Under H_0,
\sqrt{n}\,(\overline{\mathbf{X}}-\boldsymbol{\mu}_0)\sim\operatorname{N}_p(\mathbf{0},\Sigma)
is independent of (n-1)S=A\sim W_p(\Sigma,n-1). For that pair,
F=\frac{n-p}{p(n-1)}\,T^2\sim F_{p,\,n-p}.
The F distribution is the same relationship that, when p=1, turns a normal numerator and an independent chi-squared denominator into a t statistic, because the square of a t_{n-1} variable is F_{1,n-1}. In that case
F=\frac{n-1}{1\cdot(n-1)}T^2=T^2=t^2,
so the F test is the two-sided t test. Reject H_0 at level \alpha when F>F_{\alpha;\,p,\,n-p}. The denominator degrees of freedom are n-p, so the classical statistic needs n>p. If n\le p then A is singular and T^2, written with S^{-1}, is undefined.
Theorem 1 (Null distribution of Hotelling’s T^2) For n>p independent observations from \operatorname{N}_p(\boldsymbol{\mu},\Sigma), the statistic
T^2=n(\overline{\mathbf{X}}-\boldsymbol{\mu}_0)^\top S^{-1}(\overline{\mathbf{X}}-\boldsymbol{\mu}_0), \qquad S=\frac{1}{n-1}A,
satisfies
F=\frac{n-p}{p(n-1)}\,T^2\sim F_{p,\,n-p}
when \boldsymbol{\mu}=\boldsymbol{\mu}_0. Reject H_0:\boldsymbol{\mu}=\boldsymbol{\mu}_0 at level \alpha when F>F_{\alpha;\,p,\,n-p}.
T^2 is the squared Mahalanobis distance of the sample mean from the hypothesised mean, in units of the sample covariance, multiplied by n. It is a strictly decreasing function of the likelihood ratio. When p=1 it is the square of the ordinary t statistic.
A confidence region for the mean
Main article: Confidence region for a normal mean vector
A confidence region can now be built without inventing another statistic. Run the level-\alpha test for every proposed mean and retain the values it does not reject. Under the Gaussian model, the resulting random region contains the fixed true mean in a proportion 1-\alpha of repeated samples. Once the data are observed, this statement is not a probability that the parameter lies in this particular region.
n(\bar{\mathbf{x}}-\boldsymbol{\mu})^\top S^{-1}(\bar{\mathbf{x}}-\boldsymbol{\mu}) \le\frac{(n-1)p}{n-p}\,F_{\alpha;\,p,\,n-p}.
A level set of (\mathbf{x}-\boldsymbol{\mu})^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu}) is an ellipsoid whose axes are the eigenvectors of \Sigma. The same algebra, with S in place of \Sigma and with centre \bar{\mathbf{x}}, shows that the inequality is an ellipsoid centred at the sample mean, with axes along the eigenvectors of S. The quadratic form carries a factor n, so the ellipsoid describes uncertainty about the mean and shrinks as n grows. The scatter of the individual observations is what S itself measures. The numbers in the next example are a case in which this ellipsoid excludes the origin while each separate t interval still contains it. The two measurements move together, and the sample mean leaves the origin along a direction in which the sample hardly varies.
Example 1 (Two paired measurements on twelve subjects) Twelve subjects each provide the change in two correlated measurements, (x_1,x_2):
| Subject | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| x_1 | -1.8 | 6.2 | -1.5 | -0.9 | 0.3 | 1.4 | 2.6 | 1.6 | -2.3 | 0.5 | 3.4 | 0.9 |
| x_2 | -2.1 | 3.8 | -2.2 | -1.4 | -1.3 | -1.9 | 0.8 | 0.8 | -2.5 | -0.2 | 0.3 | -1.0 |
Question. Test H_0:\boldsymbol{\mu}=\mathbf{0} at the 5\% level, and compare the result with two separate t tests.
Solution. The sums are \sum x_1=10.4 and \sum x_2=-6.9, so
\bar{\mathbf{x}}=(10.4/12,\ -6.9/12)^\top=(0.8667,\ -0.5750)^\top.
The corrected cross-product matrix and the sample covariance S=A/11 are
A=\begin{pmatrix}65.0067&43.9600\\43.9600&35.6425\end{pmatrix}, \qquad S=\begin{pmatrix}5.9097&3.9964\\3.9964&3.2402\end{pmatrix}.
The sample correlation is
r=\frac{3.9964}{\sqrt{5.9097\times 3.2402}}=0.913.
Then
|S|=5.9097\times 3.2402-3.9964^2=3.1778
and
S^{-1}=\frac{1}{3.1778}\begin{pmatrix}3.2402&-3.9964\\-3.9964&5.9097\end{pmatrix} =\begin{pmatrix}1.0196&-1.2576\\-1.2576&1.8597\end{pmatrix}.
Hence
\begin{aligned} T^2 &=12\,\bar{\mathbf{x}}^\top S^{-1}\bar{\mathbf{x}}\\ &=12\Bigl[1.0196(0.8667)^2-2(1.2576)(0.8667)(-0.5750)+1.8597(0.5750)^2\Bigr]\\ &=12\times 2.6342=31.61. \end{aligned}
The three products round to 0.7659, 1.2535 and 0.6149. Their unrounded sum is 2.6342, and 12\times 2.6342=31.61.
The scaled statistic is
F=\frac{12-2}{2\times 11}\times 31.61=\frac{10}{22}\times 31.61=14.37
on (2,10) degrees of freedom. The critical value is F_{0.05;\,2,\,10}=4.103 and the p-value is 0.0011, so H_0 is rejected.
The separate t tests give
t_1=\frac{0.8667}{\sqrt{5.9097/12}}=1.235 \quad(p=0.243), \qquad t_2=\frac{-0.5750}{\sqrt{3.2402/12}}=-1.107 \quad(p=0.292).
Neither coordinate is individually significant. The two coordinates are strongly positively correlated (r=0.913) while the sample mean departs from the origin in opposite directions, (+,-). Such a departure lies along the shorter principal axis of the data cloud, where there is little variability. The T^2 test takes the correlation into account. The marginal tests cannot. The 95\% Hotelling ellipsoid for the mean excludes the origin. The rectangle of Bonferroni-adjusted t intervals includes it.
Used in
See also
One-sample t test · Univariate t test as a template for T² · Iris data · Invariance of the maximum likelihood estimator · Wishart distribution · Sampling law of the normal mean and scatter · Multivariate central limit theorem · The limit of Hotelling’s statistic