Let \(\mathbf{Z}_1, \ldots, \mathbf{Z}_m\) be independent \(N_p(\mathbf{0}, \boldsymbol\Sigma)\) vectors. The random matrix
\[ \mathbf{A} = \sum_{i=1}^{m} \mathbf{Z}_i\mathbf{Z}_i' \]is called a Wishart matrix, and its distribution is the Wishart distribution \(W_p(m, \boldsymbol\Sigma)\) with \(m\) degrees of freedom and scale matrix \(\boldsymbol\Sigma\).
Start from \(p = 1\). There \(\mathbf{Z}_i\) is a scalar \(Z_i \sim N(0, \sigma^{2})\), and \(A = \sum_i Z_i^{2}\). Writing \(Z_i = \sigma U_i\) with \(U_i \sim N(0,1)\) gives \(A = \sigma^{2}\sum_i U_i^{2} \sim \sigma^{2}\chi^{2}_{m}\). So
\[ W_1(m, \sigma^{2}) = \sigma^{2}\chi^{2}_{m}. \]The Wishart is the chi-square, generalised to matrices. Every property below is the matrix version of something already known about \(\chi^{2}\).
(i) The mean. \(E(\mathbf{A}) = m\boldsymbol\Sigma\). Proof: \(E(\mathbf{Z}_i\mathbf{Z}_i') = \boldsymbol\Sigma\) by the definition of a dispersion matrix about a zero mean, and expectation is additive over the \(m\) independent terms. (\(p = 1\): \(E(\chi^{2}_m) = m\).)
(ii) Additivity. If \(\mathbf{A}_1 \sim W_p(m_1, \boldsymbol\Sigma)\) and \(\mathbf{A}_2 \sim W_p(m_2, \boldsymbol\Sigma)\) are independent with the same \(\boldsymbol\Sigma\), then \(\mathbf{A}_1 + \mathbf{A}_2 \sim W_p(m_1+m_2, \boldsymbol\Sigma)\). Proof: concatenate the two lists of \(\mathbf{Z}\) vectors; the sum of the outer products over the combined list is the sum of the two matrices. (\(p = 1\): chi-squares add.) This is exactly what licenses pooling two sample covariance matrices in Unit 3.
(iii) Linear transformation. If \(\mathbf{A} \sim W_p(m, \boldsymbol\Sigma)\) and \(\mathbf{C}\) is a constant \(q \times p\) matrix, then \(\mathbf{CAC}' \sim W_q(m, \mathbf{C}\boldsymbol\Sigma\mathbf{C}')\). Proof: \(\mathbf{CAC}' = \sum_i (\mathbf{CZ}_i)(\mathbf{CZ}_i)'\), and each \(\mathbf{CZ}_i\) is \(N_q(\mathbf{0}, \mathbf{C}\boldsymbol\Sigma\mathbf{C}')\) by Result (a) of Unit 1.
(iv) Any fixed direction gives a chi-square. Taking \(\mathbf{C} = \mathbf{a}'\) in (iii),
\[ \frac{\mathbf{a}'\mathbf{A}\mathbf{a}}{\mathbf{a}'\boldsymbol\Sigma\mathbf{a}} \sim \chi^{2}_{m} \qquad \text{for every fixed } \mathbf{a} \ne \mathbf{0}. \]In particular each diagonal entry \(a_{jj}/\sigma_{jj} \sim \chi^{2}_m\): the Wishart's diagonal carries the familiar univariate result, one variable at a time. The words for every fixed \(\mathbf{a}\) are essential — if \(\mathbf{a}\) is chosen after seeing the data, as it is in discriminant analysis, the result fails, and that failure is precisely why Hotelling's \(T^{2}\) exists.
For \(m \ge p\) and \(\boldsymbol\Sigma\) positive definite the Wishart density is
\[ f(\mathbf{A}) = \frac{|\mathbf{A}|^{(m-p-1)/2} \exp\left\{-\tfrac12\operatorname{tr}(\boldsymbol\Sigma^{-1}\mathbf{A})\right\}} {2^{mp/2}|\boldsymbol\Sigma|^{m/2}\,\Gamma_p(m/2)}, \qquad \mathbf{A} \text{ positive definite}, \]where \(\Gamma_p\) is the multivariate gamma function \(\Gamma_p(u) = \pi^{p(p-1)/4}\prod_{j=1}^{p}\Gamma\!\left(u - \tfrac{j-1}{2}\right)\). At \(p = 1\) this collapses to the \(\sigma^{2}\chi^{2}_{m}\) density, term for term.
If \(m < p\) the matrix \(\mathbf{A}\) is singular with probability 1 — it is a sum of \(m\) rank-one matrices, so its rank cannot exceed \(m\) — and no density on the space of positive definite matrices exists. This is the reason a multivariate analysis needs more observations than variables, and it is a hard constraint rather than a rule of thumb.
Given. The four bivariate observations of Example 1.3, from which
\[ \mathbf{A} = \begin{pmatrix}20 & 16\\ 16 & 20\end{pmatrix}. \]Step 1 — name its distribution. By Result (ii) of Unit 1, section 5, \(\mathbf{A} \sim W_2(n-1, \boldsymbol\Sigma) = W_2(3, \boldsymbol\Sigma)\). One degree of freedom has gone into \(\bar{\mathbf{x}}\), exactly as \(n-1\) appears in the univariate \(S^{2}\).
Step 2 — check the degrees of freedom are enough. Here \(m = 3\) and \(p = 2\), so \(m \ge p\) and \(\mathbf{A}\) is non-singular with probability 1. Indeed \(|\mathbf{A}| = 400 - 256 = 144 \ne 0\). With only two observations we would have had \(m = 1 < 2 = p\) and a singular \(\mathbf{A}\), whatever the data.
Step 3 — use property (i). \(E(\mathbf{A}) = 3\boldsymbol\Sigma\), so \(\mathbf{A}/3\) is unbiased for \(\boldsymbol\Sigma\) — which is the matrix \(\mathbf{S}\) of Example 1.3, and is the reason for the divisor \(n-1\).
Step 4 — use property (iv) on a direction chosen in advance. Take \(\mathbf{a} = (1, -1)'\), the difference of the two variables. Then
\[ \mathbf{a}'\mathbf{A}\mathbf{a} = 20 - 16 - 16 + 20 = 8, \]and \(8/\mathbf{a}'\boldsymbol\Sigma\mathbf{a} \sim \chi^{2}_{3}\). If a null hypothesis specified \(\mathbf{a}'\boldsymbol\Sigma\mathbf{a} = 4\), the observed value would be \(8/4 = 2\), against \(E(\chi^{2}_3) = 3\) — unremarkable.
Interpretation. The Wishart is not an extra distribution to memorise. It is the bookkeeping device that lets the single matrix \(\mathbf{A}\) answer chi-square questions about any pre-specified linear combination of the variables, instead of forcing a separate calculation for each.
The generalized variance of a sample is \(|\mathbf{S}|\), and the corresponding quantity for the sum-of-squares matrix is \(|\mathbf{A}|\). If \(\mathbf{A} \sim W_p(m, \boldsymbol\Sigma)\) with \(m \ge p\), then
\[ \frac{|\mathbf{A}|}{|\boldsymbol\Sigma|} \;\stackrel{d}{=}\; \chi^{2}_{m}\cdot\chi^{2}_{m-1}\cdots\chi^{2}_{m-p+1}, \]a product of \(p\) independent chi-square variables with degrees of freedom falling by one at each step. The proof is by the Bartlett decomposition: write \(\mathbf{A} = \mathbf{TT}'\) with \(\mathbf{T}\) lower triangular; the diagonal entries \(t_{jj}^{2}\) are independent \(\chi^{2}_{m-j+1}\) variables, the off-diagonal entries are independent standard normals, and \(|\mathbf{A}| = \prod_j t_{jj}^{2}\) because the determinant of a triangular matrix is the product of its diagonal.
Two immediate consequences, since independent expectations multiply:
\[ E\!\left(\frac{|\mathbf{A}|}{|\boldsymbol\Sigma|}\right) = m(m-1)\cdots(m-p+1), \qquad E|\mathbf{S}| = \frac{m(m-1)\cdots(m-p+1)}{m^{p}}\,|\boldsymbol\Sigma| \quad \text{with } \mathbf{S} = \mathbf{A}/m. \]At \(p = 1\) the product has one factor and this is \(E(S^{2}) = \sigma^{2}\). At \(p > 1\) the factor is less than 1, so \(|\mathbf{S}|\) underestimates \(|\boldsymbol\Sigma|\) — and increasingly so as \(p\) grows. A determinant is a product of \(p\) quantities each estimated with error, and the errors compound.
Given. \(p = 2\), \(n = 4\), so \(m = n - 1 = 3\), as in Example 2.1.
Step 1 — write the product. \(|\mathbf{A}|/|\boldsymbol\Sigma| \stackrel{d}{=} \chi^{2}_{3}\cdot\chi^{2}_{2}\), the two factors independent.
Step 2 — its expectation. \(E = 3 \times 2 = 6\).
Step 3 — the bias of \(|\mathbf{S}|\). With \(\mathbf{S} = \mathbf{A}/3\) and \(p = 2\), \(|\mathbf{S}| = |\mathbf{A}|/3^{2}\), so
\[ E|\mathbf{S}| = \frac{6}{9}\,|\boldsymbol\Sigma| = 0.666667\,|\boldsymbol\Sigma|. \]Step 4 — read how bad that is. The unbiased \(\mathbf{S}\) produces a generalized variance that is on average a third too small. Each diagonal entry of \(\mathbf{S}\) is unbiased; the determinant of an unbiased matrix is not an unbiased determinant, because the determinant is a non-linear function.
Step 5 — watch it worsen with \(p\). With \(m = 10\):
| \(p\) | \(m(m-1)\cdots(m-p+1)\) | \(m^{p}\) | \(E|\mathbf{S}|/|\boldsymbol\Sigma|\) |
|---|---|---|---|
| 1 | 10 | 10 | 1.000000 |
| 2 | 90 | 100 | 0.900000 |
| 3 | 720 | 1000 | 0.720000 |
| 5 | 30,240 | 100,000 | 0.302400 |
Interpretation. At \(p = 5\) with ten degrees of freedom the generalized variance is on average less than a third of the truth. This is the first appearance of a theme that dominates high-dimensional statistics: quantities that behave perfectly well one variable at a time can degrade sharply when \(p\) is not small compared with \(n\).
Let \(r\) be the sample correlation from \(n\) pairs drawn from a bivariate normal with population correlation \(\rho\). When \(\rho = 0\), the density of \(r\) is
\[ f(r) = \frac{\Gamma\!\left(\frac{n-1}{2}\right)} {\sqrt{\pi}\;\Gamma\!\left(\frac{n-2}{2}\right)} \left(1 - r^{2}\right)^{(n-4)/2}, \qquad -1 < r < 1, \]which is symmetric about 0 and free of every unknown parameter — which is exactly what makes a table possible.
The usable form. Substituting \(t = r\sqrt{n-2}/\sqrt{1-r^{2}}\) turns that density into Student's \(t\) with \(n-2\) degrees of freedom, so
\[ t = \frac{r\sqrt{n-2}}{\sqrt{1-r^{2}}} \sim t_{n-2} \qquad \text{under } H_0: \rho = 0. \]This is the same \(t\) as the test of a regression slope, and it is the same test: the slope is zero precisely when the correlation is.
When \(\rho \ne 0\) the density of \(r\) is skewed and depends on \(\rho\), so the \(t\) transformation fails. Fisher's transformation repairs it:
\[ z = \tfrac12\ln\frac{1+r}{1-r} = \tanh^{-1} r \;\;\dot\sim\;\; N\!\left(\zeta = \tanh^{-1}\rho,\ \frac{1}{n-3}\right). \]Two features earn its keep: the variance no longer involves \(\rho\), and the skewness is removed. Intervals are built on the \(z\) scale and transformed back with \(r = \tanh z\), which keeps them inside \((-1, 1)\) automatically.
Given. \(r = 0.8\) from \(n = 20\) pairs.
(a) Test \(H_0: \rho = 0\).
\[ \sqrt{n-2} = \sqrt{18} = 4.242641, \qquad \sqrt{1 - 0.64} = \sqrt{0.36} = 0.6, \] \[ t = \frac{0.8 \times 4.242641}{0.6} = \frac{3.394113}{0.6} = 5.6569, \]on 18 degrees of freedom, with two-sided \(p\)-value \(0.0000229\). Reject decisively.
(b) Test \(H_0: \rho = 0.5\). The \(t\) test cannot be used, because its derivation assumed \(\rho = 0\). Use Fisher's \(z\):
\[ z = \tfrac12\ln\frac{1.8}{0.2} = \tfrac12\ln 9 = 1.098612, \qquad \zeta_0 = \tfrac12\ln\frac{1.5}{0.5} = \tfrac12\ln 3 = 0.549306, \] \[ SE = \frac{1}{\sqrt{n-3}} = \frac{1}{\sqrt{17}} = 0.2425356, \] \[ Z = \frac{1.098612 - 0.549306}{0.2425356} = \frac{0.549306}{0.2425356} = 2.264847, \]with two-sided \(p\)-value \(0.023522\). Reject at \(5\%\) but not at \(1\%\): the data are consistent with a correlation well above \(0.5\), but not overwhelmingly so.
(c) A \(95\%\) confidence interval for \(\rho\). On the \(z\) scale,
\[ 1.098612 \pm 1.96 \times 0.2425356 = 1.098612 \pm 0.475370 = (0.623242,\ 1.573982). \]Transforming back with \(r = \tanh z\),
\[ \left(\tanh 0.623242,\ \tanh 1.573982\right) = (0.5534,\ 0.9177). \]Interpretation. The interval is strikingly asymmetric about \(r = 0.8\): it reaches \(0.247\) below and only \(0.118\) above, because the transformation compresses the region near \(\pm 1\) where a correlation cannot go. An interval computed as \(r \pm 1.96\,SE\) on the \(r\) scale would have been symmetric, and wrong — and with \(r\) closer to 1 it would have run past 1 altogether. Note also that \(0.5\) lies outside the interval, agreeing with the rejection in (b), as it must.
Partition as in Unit 1, with \(X_1\) singled out and \(\mathbf{X}_2 = (X_2, \ldots, X_p)'\) the rest.
The multiple correlation \(R_{1\cdot 23\ldots p}\) is the largest correlation between \(X_1\) and any linear combination of the others, and
\[ R^{2}_{1\cdot 23\ldots p} = \frac{\boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\boldsymbol\Sigma_{21}} {\sigma_{11}} = 1 - \frac{\sigma_{11\cdot 2}}{\sigma_{11}} = 1 - \frac{|\boldsymbol\Sigma|}{\sigma_{11}|\boldsymbol\Sigma_{22}|}. \]The last equality is worth keeping: it computes a multiple correlation from two determinants, with no matrix inverse at all.
The partial correlation \(\rho_{ij\cdot\text{rest}}\) is the ordinary correlation between \(X_i\) and \(X_j\) computed from the conditional dispersion matrix \(\boldsymbol\Sigma_{11\cdot 2}\) — that is, the correlation that remains once the other variables have been held fixed. In the three-variable case this reduces to the formula of Statistical Methods, Unit 3:
\[ \rho_{12\cdot 3} = \frac{\rho_{12} - \rho_{13}\rho_{23}} {\sqrt{1-\rho_{13}^{2}}\sqrt{1-\rho_{23}^{2}}}. \]Their null distributions. Both are pleasantly simple:
The second is the overall \(F\) of a regression table, reached from a different direction. Note that \(R\) is never negative and its null distribution is not symmetric, so the test is one-sided in \(F\) by construction.
Given. From \(n = 28\) observations,
\[ \mathbf{R} = \begin{pmatrix}1 & 0.6 & 0.5\\ 0.6 & 1 & 0.4\\ 0.5 & 0.4 & 1\end{pmatrix}. \]Asked. \(R_{1\cdot 23}\) and \(r_{12\cdot 3}\), each with a test.
Step 1 — the determinant of \(\mathbf{R}\). Expanding along the first row,
\[ |\mathbf{R}| = 1(1 - 0.16) - 0.6(0.6 - 0.2) + 0.5(0.24 - 0.5) \] \[ = 0.84 - 0.6(0.4) + 0.5(-0.26) = 0.84 - 0.24 - 0.13 = 0.47. \]Step 2 — the determinant of the block for \(X_2, X_3\).
\[ |\mathbf{R}_{22}| = \begin{vmatrix}1 & 0.4\\ 0.4 & 1\end{vmatrix} = 1 - 0.16 = 0.84. \]Step 3 — the multiple correlation, by determinants. With \(\sigma_{11} = 1\) on the correlation scale,
\[ R^{2}_{1\cdot 23} = 1 - \frac{0.47}{0.84} = 1 - 0.559524 = 0.440476, \qquad R_{1\cdot 23} = \sqrt{0.440476} = 0.663684. \]Step 4 — test it. Here \(k = 2\) and \(n - k - 1 = 25\):
\[ F = \frac{0.440476/2}{0.559524/25} = \frac{0.220238}{0.022381} = 9.8404, \]on \((2, 25)\) degrees of freedom, with \(p\)-value \(0.000704\). The two variables together explain \(44\%\) of the variance of \(X_1\), and that is far more than chance.
Step 5 — the partial correlation.
\[ r_{12\cdot 3} = \frac{0.6 - (0.5)(0.4)}{\sqrt{1 - 0.25}\sqrt{1 - 0.16}} = \frac{0.6 - 0.2}{\sqrt{0.75 \times 0.84}} = \frac{0.4}{\sqrt{0.63}} = \frac{0.4}{0.793725} = 0.503953. \]Step 6 — test it. One variable is held fixed, so \(k = 1\) and the degrees of freedom are \(n - 2 - 1 = 25\):
\[ t = \frac{0.503953\sqrt{25}}{\sqrt{1 - 0.253969}} = \frac{0.503953 \times 5}{\sqrt{0.746031}} = \frac{2.519765}{0.863731} = 2.9173, \]with two-sided \(p\)-value \(0.007359\). The association between \(X_1\) and \(X_2\) survives controlling for \(X_3\).
Interpretation. The simple correlation \(r_{12} = 0.6\) fell to \(0.504\) once \(X_3\) was held fixed, because part of what \(X_1\) and \(X_2\) share is shared with \(X_3\). It fell, but it did not vanish — had the data been generated with \(X_3\) as the only common cause, the partial correlation would have been near zero. Distinguishing those two pictures from one correlation matrix is the whole business of path analysis in Unit 4.
Regress \(X_1\) on \(\mathbf{X}_2\). The sample coefficient vector is
\[ \hat{\boldsymbol\beta} = \mathbf{A}_{22}^{-1}\mathbf{a}_{21}, \]with \(\mathbf{A}\) the sum-of-squares-and-cross-products matrix partitioned as \(\boldsymbol\Sigma\) was. Conditional on \(\mathbf{X}_2\), this is exactly the least squares estimator of Linear Algebra and Linear Models, Unit 4, so the Gauss–Markov theorem and everything proved there applies unchanged, and
\[ \hat{\boldsymbol\beta} \sim N\!\left(\boldsymbol\beta,\ \sigma_{11\cdot 2}\,\mathbf{A}_{22}^{-1}\right), \qquad \frac{(n - p)\,\hat\sigma_{11\cdot 2}}{\sigma_{11\cdot 2}} \sim \chi^{2}_{n-p}, \]the two being independent. Dividing a normal by the square root of an independent chi-square over its degrees of freedom gives, for the \(j\)th coefficient,
\[ \frac{\hat\beta_j - \beta_j}{\sqrt{\hat\sigma_{11\cdot 2}\,(\mathbf{A}_{22}^{-1})_{jj}}} \sim t_{n-p}, \]which is the whole of inference about regression coefficients: a test at \(\beta_j = 0\) and a confidence interval \(\hat\beta_j \pm t\,SE\).
The conditional argument is the point. Nothing above requires \(\mathbf{X}_2\) to be random. If it is, the results hold conditionally on it and therefore unconditionally as well, since the distribution quoted does not depend on the conditioning values. That is why the regression theory of a fixed-design linear model can be used verbatim when all the variables are jointly normal.
Given. \(n = 20\) observations on \(X_1\) (dependent), \(X_2\) and \(X_3\), with corrected sums of squares and products
\[ \mathbf{A} = \begin{pmatrix}100 & 60 & 40\\ 60 & 90 & 30\\ 40 & 30 & 80\end{pmatrix}. \]Step 1 — partition and invert.
\[ \mathbf{A}_{22} = \begin{pmatrix}90 & 30\\ 30 & 80\end{pmatrix}, \qquad |\mathbf{A}_{22}| = 7200 - 900 = 6300, \qquad \mathbf{A}_{22}^{-1} = \frac{1}{6300}\begin{pmatrix}80 & -30\\ -30 & 90\end{pmatrix}. \]Step 2 — the coefficients.
\[ \hat{\boldsymbol\beta} = \frac{1}{6300}\begin{pmatrix}80 & -30\\ -30 & 90\end{pmatrix} \begin{pmatrix}60\\40\end{pmatrix} = \frac{1}{6300}\begin{pmatrix}4800 - 1200\\ -1800 + 3600\end{pmatrix} = \frac{1}{6300}\begin{pmatrix}3600\\1800\end{pmatrix} = \begin{pmatrix}0.571429\\ 0.285714\end{pmatrix}. \]Step 3 — the sums of squares.
\[ SS_{\text{reg}} = \hat{\boldsymbol\beta}'\mathbf{a}_{21} = \tfrac{4}{7}(60) + \tfrac{2}{7}(40) = \frac{240 + 80}{7} = \frac{320}{7} = 45.714286, \] \[ SS_{\text{res}} = a_{11} - SS_{\text{reg}} = 100 - \frac{320}{7} = \frac{380}{7} = 54.285714. \]Step 4 — the residual mean square. The degrees of freedom are \(n - p = 20 - 3 = 17\), one for the intercept and one for each coefficient:
\[ \hat\sigma_{11\cdot 2} = \frac{380/7}{17} = \frac{380}{119} = 3.193277. \]Step 5 — the standard errors.
\[ \operatorname{Var}(\hat\beta_1) = 3.193277 \times \frac{80}{6300} = 0.040550, \qquad SE(\hat\beta_1) = 0.201369, \] \[ \operatorname{Var}(\hat\beta_2) = 3.193277 \times \frac{90}{6300} = 0.045618, \qquad SE(\hat\beta_2) = 0.213584. \]Step 6 — the two \(t\) statistics, on 17 degrees of freedom.
\[ t_1 = \frac{0.571429}{0.201369} = 2.8377 \;\;(p = 0.0114), \qquad t_2 = \frac{0.285714}{0.213584} = 1.3377 \;\;(p = 0.1986). \]Step 7 — a \(95\%\) interval for \(\beta_1\). With \(t_{17,\,0.025} = 2.109816\),
\[ 0.571429 \pm 2.109816 \times 0.201369 = 0.571429 \pm 0.424852 = (0.146577,\ 0.996281). \]Step 8 — the overall test, for comparison.
\[ R^{2} = \frac{320/7}{100} = \frac{16}{35} = 0.457143, \qquad F = \frac{R^{2}/2}{(1 - R^{2})/17} = \frac{8/35}{19/595} = \frac{8 \times 17}{19} = \frac{136}{19} = 7.158, \]on \((2, 17)\) degrees of freedom, \(p = 0.0056\).
Interpretation. The model as a whole is significant, and \(X_2\) is carrying it; \(X_3\) adds nothing once \(X_2\) is present, despite having a simple correlation with \(X_1\) of \(40/\sqrt{100 \times 80} = 0.447\). That is not a contradiction: \(X_2\) and \(X_3\) are themselves correlated \(\left(30/\sqrt{90 \times 80} = 0.354\right)\), so they compete for the same explanation. A \(t\) statistic in a multiple regression always answers the question “does this variable add anything to the others?” and never “does this variable matter?”