Given \(\mathbf{X}\) with dispersion matrix \(\boldsymbol\Sigma\), find the linear combination \(Y_1 = \mathbf{e}_1'\mathbf{X}\) of maximum variance. Without a restriction the problem has no answer, since doubling \(\mathbf{e}_1\) quadruples the variance, so require \(\mathbf{e}_1'\mathbf{e}_1 = 1\). The problem is
\[ \max_{\mathbf{e}'\mathbf{e}=1} \operatorname{Var}(\mathbf{e}'\mathbf{X}) = \max_{\mathbf{e}'\mathbf{e}=1} \mathbf{e}'\boldsymbol\Sigma\mathbf{e} \quad\Longleftrightarrow\quad \max_{\mathbf{e}\ne\mathbf{0}} \frac{\mathbf{e}'\boldsymbol\Sigma\mathbf{e}}{\mathbf{e}'\mathbf{e}}, \]which is exactly the Rayleigh quotient. Its maximum is the largest eigenvalue \(\lambda_1\), attained at the corresponding eigenvector. So
\[ Y_i = \mathbf{e}_i'\mathbf{X}, \qquad \operatorname{Var}(Y_i) = \lambda_i, \qquad \operatorname{Cov}(Y_i, Y_j) = 0 \;\; (i \ne j), \]the last because eigenvectors of a symmetric matrix belonging to distinct eigenvalues are orthogonal: \(\operatorname{Cov}(Y_i, Y_j) = \mathbf{e}_i'\boldsymbol\Sigma\mathbf{e}_j = \lambda_j\mathbf{e}_i'\mathbf{e}_j = 0\). The components are uncorrelated by construction, not by assumption.
Three consequences.
which is how a component is interpreted: by seeing which variables it correlates with.
\(\boldsymbol\Sigma\) or \(\boldsymbol\rho\)? Principal components are not invariant under rescaling, because a variance is not. If one variable is measured in millimetres and another in kilometres, the first will dominate every component for a reason that has nothing to do with the data. Unless all the variables are in the same units and of comparable magnitude, analyse the correlation matrix — whose trace is \(p\), so each component is compared against the "one variable's worth" benchmark of \(\lambda = 1\).
Given. \(\boldsymbol\Sigma = \begin{pmatrix}5 & 2\\ 2 & 2\end{pmatrix}\).
Step 1 — the characteristic equation. For a \(2 \times 2\) matrix it is \(\lambda^{2} - (\operatorname{tr})\lambda + \det = 0\), with \(\operatorname{tr} = 7\) and \(\det = 10 - 4 = 6\):
\[ \lambda^{2} - 7\lambda + 6 = 0 \;\Longrightarrow\; (\lambda - 6)(\lambda - 1) = 0 \;\Longrightarrow\; \lambda_1 = 6, \;\; \lambda_2 = 1. \]Step 2 — check against the trace. \(6 + 1 = 7 = 5 + 2\). \(\checkmark\)
Step 3 — the first eigenvector. Solve \((\boldsymbol\Sigma - 6\mathbf{I})\mathbf{e} = \mathbf{0}\):
\[ \begin{pmatrix}-1 & 2\\ 2 & -4\end{pmatrix}\begin{pmatrix}e_1\\e_2\end{pmatrix} = \begin{pmatrix}0\\0\end{pmatrix} \;\Longrightarrow\; -e_1 + 2e_2 = 0 \;\Longrightarrow\; e_1 = 2e_2. \]The second row \(2e_1 - 4e_2 = 0\) says the same thing, as it must for a singular matrix. Taking \(\mathbf{e} \propto (2, 1)'\) and normalising by \(\sqrt{4 + 1} = \sqrt 5\),
\[ \mathbf{e}_1 = \frac{1}{\sqrt 5}\begin{pmatrix}2\\1\end{pmatrix} = \begin{pmatrix}0.894427\\ 0.447214\end{pmatrix}. \]Step 4 — the second eigenvector. By orthogonality it must be \(\mathbf{e}_2 \propto (1, -2)'\), that is
\[ \mathbf{e}_2 = \frac{1}{\sqrt 5}\begin{pmatrix}1\\-2\end{pmatrix} = \begin{pmatrix}0.447214\\ -0.894427\end{pmatrix}, \]and indeed \(\mathbf{e}_1'\mathbf{e}_2 = (2 - 2)/5 = 0\).
Step 5 — the components and how much they carry.
\[ Y_1 = 0.894427X_1 + 0.447214X_2, \qquad Y_2 = 0.447214X_1 - 0.894427X_2, \] \[ \frac{\lambda_1}{\lambda_1 + \lambda_2} = \frac{6}{7} = 0.857143, \qquad \frac{\lambda_2}{\lambda_1+\lambda_2} = \frac{1}{7} = 0.142857. \]Step 6 — the loadings. With \(\sqrt{\lambda_1} = \sqrt 6 = 2.449490\), \(\sqrt{\sigma_{11}} = \sqrt 5 = 2.236068\) and \(\sqrt{\sigma_{22}} = \sqrt 2 = 1.414214\),
\[ \operatorname{Corr}(Y_1, X_1) = \frac{0.894427 \times 2.449490}{2.236068} = 0.979796, \qquad \operatorname{Corr}(Y_1, X_2) = \frac{0.447214 \times 2.449490}{1.414214} = 0.774597. \]Interpretation. One component of the two carries \(85.7\%\) of the variance and correlates \(0.98\) with \(X_1\) and \(0.77\) with \(X_2\): both variables load on it positively, so \(Y_1\) is a "size" or "overall level" component, and \(Y_2\), which contrasts them, is a "shape" component. That pattern — a first component on which everything loads positively, later ones that contrast — is what a positively correlated set of variables almost always produces.
Given. The matrix of Unit 2, Example 2.4:
\[ \mathbf{R} = \begin{pmatrix}1 & 0.6 & 0.5\\ 0.6 & 1 & 0.4\\ 0.5 & 0.4 & 1\end{pmatrix}. \]Step 1 — build the characteristic equation from the invariants. For a \(3 \times 3\) matrix,
\[ \lambda^{3} - c_1\lambda^{2} + c_2\lambda - c_3 = 0, \]where \(c_1 = \operatorname{tr}\mathbf{R}\), \(c_2\) is the sum of the three \(2 \times 2\) principal minors, and \(c_3 = |\mathbf{R}|\). Here
\[ c_1 = 3, \qquad c_2 = (1 - 0.36) + (1 - 0.25) + (1 - 0.16) = 0.64 + 0.75 + 0.84 = 2.23, \]and \(c_3 = 0.47\) from Example 2.4, so the equation is
\[ \lambda^{3} - 3\lambda^{2} + 2.23\lambda - 0.47 = 0. \]Step 2 — solve it. The cubic does not factor over the rationals, so the roots are found numerically:
\[ \lambda_1 = 2.004458, \qquad \lambda_2 = 0.613092, \qquad \lambda_3 = 0.382451. \]Step 3 — two checks before going on.
\[ \lambda_1 + \lambda_2 + \lambda_3 = 3.000000 = \operatorname{tr}\mathbf{R}, \qquad \lambda_1\lambda_2\lambda_3 = 0.470000 = |\mathbf{R}|. \checkmark \]Step 4 — the eigenvectors and loadings.
| \(\lambda\) | \(\mathbf{e}\) | loadings \(e_k\sqrt\lambda\) | proportion | cumulative | |
|---|---|---|---|---|---|
| PC1 | 2.004458 | (0.613324, 0.579903, 0.536233) | (0.868338, 0.821020, 0.759193) | 0.668153 | 0.668153 |
| PC2 | 0.613092 | (−0.179172, −0.559070, 0.809530) | (−0.140292, −0.437752, 0.633863) | 0.204364 | 0.872516 |
| PC3 | 0.382451 | (0.769240, −0.592582, −0.238989) | (0.475718, −0.366468, −0.147797) | 0.127484 | 1.000000 |
Because \(\mathbf{R}\) is a correlation matrix, the loadings are directly the correlations between the components and the variables, with no further division.
Step 5 — how many components to keep. Two rules of thumb, and what each says here:
Interpretation. The two rules disagree, which is usual and is why neither is a substitute for looking at the loadings. PC1 loads \(0.87\), \(0.82\), \(0.76\) — positively and almost equally on all three variables, so it is the common factor running through them, and it alone reproduces \(67\%\) of the correlation structure. PC2 contrasts \(X_3\) against \(X_1\) and \(X_2\). Compare this with the multiple correlation \(R_{1\cdot 23} = 0.664\) found for the same matrix in Example 2.4: both are saying that these three variables share a single dominant dimension, from different directions.
Two sets of variables are measured on the same units — say \(p\) aptitude scores and \(q\) job-performance scores. A correlation matrix between the sets has \(pq\) entries, which is too many to read. Canonical correlation analysis asks for the single pair of linear combinations, one from each set, that are most highly correlated:
\[ U = \mathbf{a}'\mathbf{X}^{(1)}, \qquad V = \mathbf{b}'\mathbf{X}^{(2)}, \qquad \max_{\mathbf{a},\mathbf{b}} \operatorname{Corr}(U, V), \]then for the next such pair uncorrelated with the first, and so on, giving \(\min(p, q)\) pairs in all.
The solution. Partition the correlation matrix as
\[ \mathbf{R} = \begin{pmatrix}\mathbf{R}_{11} & \mathbf{R}_{12}\\ \mathbf{R}_{21} & \mathbf{R}_{22}\end{pmatrix}. \]The squared canonical correlations \(\rho_i^{2}\) are the eigenvalues of
\[ \mathbf{M} = \mathbf{R}_{11}^{-1}\mathbf{R}_{12}\mathbf{R}_{22}^{-1}\mathbf{R}_{21}, \]and the vectors \(\mathbf{a}_i\) are its eigenvectors, scaled so that \(\mathbf{a}_i'\mathbf{R}_{11}\mathbf{a}_i = 1\) (which makes \(\operatorname{Var}(U_i) = 1\)). The vectors \(\mathbf{b}_i\) come from the mirror-image matrix \(\mathbf{R}_{22}^{-1}\mathbf{R}_{21}\mathbf{R}_{11}^{-1}\mathbf{R}_{12}\), which has the same non-zero eigenvalues.
Two special cases tie it to what is already known. If \(q = 1\), there is one canonical correlation and it is the multiple correlation of that single variable on the other set. If \(p = q = 1\), it is the ordinary correlation. Canonical correlation is therefore the common generalisation of both.
Given. Standardised variables with
\[ \mathbf{R}_{11} = \begin{pmatrix}1 & 0.5\\ 0.5 & 1\end{pmatrix}, \qquad \mathbf{R}_{22} = \begin{pmatrix}1 & 0.4\\ 0.4 & 1\end{pmatrix}, \qquad \mathbf{R}_{12} = \begin{pmatrix}0.6 & 0.3\\ 0.3 & 0.6\end{pmatrix}. \]Step 1 — the two inverses. \(|\mathbf{R}_{11}| = 0.75\) and \(|\mathbf{R}_{22}| = 0.84\), so
\[ \mathbf{R}_{11}^{-1} = \frac{4}{3}\begin{pmatrix}1 & -0.5\\ -0.5 & 1\end{pmatrix}, \qquad \mathbf{R}_{22}^{-1} = \frac{25}{21}\begin{pmatrix}1 & -0.4\\ -0.4 & 1\end{pmatrix}. \]Step 2 — the middle product.
\[ \mathbf{R}_{12}\mathbf{R}_{22}^{-1} = \frac{25}{21}\begin{pmatrix}0.6 & 0.3\\ 0.3 & 0.6\end{pmatrix} \begin{pmatrix}1 & -0.4\\ -0.4 & 1\end{pmatrix} = \frac{25}{21}\begin{pmatrix}0.48 & 0.06\\ 0.06 & 0.48\end{pmatrix}. \]Step 3 — multiply by \(\mathbf{R}_{21} = \mathbf{R}_{12}'\).
\[ \frac{25}{21}\begin{pmatrix}0.48 & 0.06\\ 0.06 & 0.48\end{pmatrix} \begin{pmatrix}0.6 & 0.3\\ 0.3 & 0.6\end{pmatrix} = \frac{25}{21}\begin{pmatrix}0.306 & 0.180\\ 0.180 & 0.306\end{pmatrix}. \]Step 4 — premultiply by \(\mathbf{R}_{11}^{-1}\).
\[ \mathbf{M} = \frac{4}{3}\cdot\frac{25}{21} \begin{pmatrix}1 & -0.5\\ -0.5 & 1\end{pmatrix} \begin{pmatrix}0.306 & 0.180\\ 0.180 & 0.306\end{pmatrix} = \frac{100}{63}\begin{pmatrix}0.216 & 0.027\\ 0.027 & 0.216\end{pmatrix} = \begin{pmatrix}\tfrac{12}{35} & \tfrac{3}{70}\\ \tfrac{3}{70} & \tfrac{12}{35}\end{pmatrix}. \]Step 5 — its eigenvalues. A matrix of the form \(\begin{pmatrix}\alpha & \beta\\ \beta & \alpha\end{pmatrix}\) has eigenvalues \(\alpha \pm \beta\) with eigenvectors \((1, 1)'\) and \((1, -1)'\), so
\[ \rho_1^{2} = \frac{12}{35} + \frac{3}{70} = \frac{24 + 3}{70} = \frac{27}{70} = 0.385714, \qquad \rho_2^{2} = \frac{12}{35} - \frac{3}{70} = \frac{21}{70} = 0.300000, \] \[ \rho_1 = \sqrt{0.385714} = 0.621059, \qquad \rho_2 = \sqrt{0.3} = 0.547723. \]Step 6 — the first canonical pair, scaled. With \(\mathbf{a} \propto (1,1)'\), write \(\mathbf{a} = c(1,1)'\) and impose \(\mathbf{a}'\mathbf{R}_{11}\mathbf{a} = 1\):
\[ c^{2}(1 + 2(0.5) + 1) = 3c^{2} = 1 \;\Longrightarrow\; c = \frac{1}{\sqrt 3} = 0.577350. \]Similarly \(\mathbf{b} = d(1,1)'\) with \(d^{2}(1 + 2(0.4) + 1) = 2.8d^{2} = 1\), so \(d = 1/\sqrt{2.8} = 0.597614\). Hence
\[ U_1 = 0.577350\left(X^{(1)}_1 + X^{(1)}_2\right), \qquad V_1 = 0.597614\left(X^{(2)}_1 + X^{(2)}_2\right). \]Step 7 — verify the correlation directly.
\[ \operatorname{Corr}(U_1, V_1) = \mathbf{a}'\mathbf{R}_{12}\mathbf{b} = \frac{1}{\sqrt 3}\cdot\frac{1}{\sqrt{2.8}}\,(0.6 + 0.3 + 0.3 + 0.6) = \frac{1.8}{\sqrt{8.4}} = \frac{1.8}{2.898275} = 0.621059. \checkmark \]Step 8 — the second pair. Now \(\mathbf{a} \propto (1, -1)'\), and \(c^{2}(1 - 2(0.5) + 1) = c^{2} = 1\) gives \(\mathbf{a} = (1, -1)'\); for the second set \(d^{2}(1 - 0.8 + 1) = 1.2d^{2} = 1\) gives \(d = 1/\sqrt{1.2} = 0.912871\). Then
\[ \operatorname{Corr}(U_2, V_2) = (1, -1)\begin{pmatrix}0.6 & 0.3\\ 0.3 & 0.6\end{pmatrix} \begin{pmatrix}0.912871\\ -0.912871\end{pmatrix} = 0.912871\,(0.3 + 0.3) = 0.547723. \checkmark \]Interpretation. The strongest link between the two sets is between their sums — the overall level in one set against the overall level in the other, correlating \(0.621\). The second link is between their differences, correlating \(0.548\). Neither is visible in the raw \(\mathbf{R}_{12}\), whose four entries are \(0.6, 0.3, 0.3, 0.6\); the analysis has replaced four numbers with two, and told you what they mean.
Discriminant analysis in Unit 3 was given the groups and asked for a rule. Clustering is given no groups and asked to find some — which means there is no right answer to check against, and no null distribution to test. It is an exploratory method, and the honest way to read a clustering is as a hypothesis about structure rather than as a finding.
Agglomerative hierarchical clustering starts with \(n\) clusters of one object each and repeatedly merges the closest pair, until one cluster remains. What "closest" means between clusters is the choice that defines the method:
| Linkage | Distance between clusters \(A\) and \(B\) | Tends to produce |
|---|---|---|
| Single (nearest neighbour) | \(\min_{i \in A,\; j \in B} d_{ij}\) | Long straggly clusters; prone to chaining, where a line of intermediate points welds two groups together |
| Complete (furthest neighbour) | \(\max_{i \in A,\; j \in B} d_{ij}\) | Compact, roughly equal-sized clusters; sensitive to outliers |
| Average | \(\dfrac{1}{|A||B|}\sum_{i \in A}\sum_{j \in B} d_{ij}\) | A compromise between the two, and the most commonly used |
The result is drawn as a dendrogram, whose vertical scale is the distance at which each merge happened. Cutting it at a chosen height gives a set of clusters; the largest vertical gap between consecutive merges is the usual, and unavoidably informal, guide to where to cut.
Given. Five objects in the plane:
\[ A(1, 1), \quad B(2, 1), \quad C(4, 5), \quad D(7, 7), \quad E(5, 7). \]Step 1 — the Euclidean distance matrix. Each entry is \(\sqrt{(\Delta x)^{2} + (\Delta y)^{2}}\); for instance \(d_{AC} = \sqrt{3^{2} + 4^{2}} = \sqrt{25} = 5\).
| A | B | C | D | E | |
|---|---|---|---|---|---|
| A | 0 | 1.000000 | 5.000000 | 8.485281 | 7.211103 |
| B | 1.000000 | 0 | 4.472136 | 7.810250 | 6.708204 |
| C | 5.000000 | 4.472136 | 0 | 3.605551 | 2.236068 |
| D | 8.485281 | 7.810250 | 3.605551 | 0 | 2.000000 |
| E | 7.211103 | 6.708204 | 2.236068 | 2.000000 | 0 |
Step 2 — first merge. The smallest entry is \(d_{AB} = 1.000000\). Merge \(\{A, B\}\) at height \(1.000000\). This first merge is the same for every linkage, because with singleton clusters all three definitions coincide.
Step 3 — second merge. The smallest remaining distance is \(d_{DE} = 2.000000\). Merge \(\{D, E\}\).
Step 4 — third merge, which is where the linkages separate. The candidates are \(C\) with \(\{D,E\}\), \(\{A,B\}\) with \(C\), and \(\{A,B\}\) with \(\{D,E\}\):
| Pair | Single | Complete | Average |
|---|---|---|---|
| \(C\) – \(\{D,E\}\) | \(\min(3.605551, 2.236068) = 2.2361\) | \(\max(3.605551, 2.236068) = 3.6056\) | \(\tfrac12(3.605551 + 2.236068) = 2.9208\) |
| \(\{A,B\}\) – \(C\) | \(4.4721\) | \(5.0000\) | \(\tfrac12(5.000000 + 4.472136) = 4.7361\) |
| \(\{A,B\}\) – \(\{D,E\}\) | \(6.7082\) | \(8.4853\) | \(\tfrac14(8.485281 + 7.211103 + 7.810250 + 6.708204) = 7.5537\) |
All three agree that \(C\) joins \(\{D, E\}\), at heights \(2.2361\), \(3.6056\) and \(2.9208\) respectively. The merge heights are given to four decimals throughout this comparison: they are used only to order and to place the merges, and the distance matrix above carries the six-decimal figures they come from.
Step 5 — last merge. Only \(\{A,B\}\) and \(\{C,D,E\}\) remain:
\[ \text{single } 4.4721, \qquad \text{complete } 8.4853, \qquad \text{average } \frac{39.686974}{6} = 6.6145. \]The average is over all six pairs \(AC, AD, AE, BC, BD, BE\).
Step 6 — read the dendrogram. Under single linkage the merge heights are \(1.0000,\; 2.0000,\; 2.2361,\; 4.4721\). The largest gap is between \(2.2361\) and \(4.4721\), so cutting there gives two clusters, \(\{A, B\}\) and \(\{C, D, E\}\).
Interpretation. All three linkages give the same tree here, and differ only in the heights — which is a comfortable outcome, and not one to rely on. Single linkage merged the last two groups at \(4.4721\) and complete linkage at \(8.4853\), nearly twice as high, because single linkage reports the closest pair of members and complete linkage the furthest. When a data set contains a chain of intermediate points the two methods produce genuinely different trees, and the choice of linkage then decides the answer.
\(K\)-means is not hierarchical: the number of clusters \(K\) is fixed in advance, and the algorithm alternates two steps until nothing changes.
Each step can only decrease the within-cluster sum of squares \(W = \sum_k \sum_{i \in C_k}\|\mathbf{x}_i - \bar{\mathbf{x}}_k\|^{2}\) — the assignment step because every object moves to a closer centroid, the update step because the mean minimises the sum of squared distances. Since there are finitely many assignments, the algorithm terminates. It converges to a local minimum, however, so the starting centroids matter and the algorithm is run from several starts.
\(K\)-nearest-neighbour clustering is a different idea using the same letter: an object is linked to its \(K\) nearest neighbours, and connected groups of links become clusters. It makes no assumption about cluster shape, which is its advantage over \(K\)-means — the latter, by using centroids and squared distances, effectively looks for roughly spherical clusters of similar size.
Given. The points of Example 4.4, \(K = 2\), starting centroids \(A(1,1)\) and \(D(7,7)\).
Step 1 — first assignment. Compare each point's distance to the two centroids, reading the numbers from the distance matrix above:
| Point | to \(A\) | to \(D\) | Cluster |
|---|---|---|---|
| A | 0 | 8.485281 | 1 |
| B | 1.000000 | 7.810250 | 1 |
| C | 5.000000 | 3.605551 | 2 |
| D | 8.485281 | 0 | 2 |
| E | 7.211103 | 2.000000 | 2 |
Step 2 — update the centroids.
\[ \bar{\mathbf{x}}_1 = \left(\frac{1+2}{2},\ \frac{1+1}{2}\right) = (1.500000,\ 1.000000), \] \[ \bar{\mathbf{x}}_2 = \left(\frac{4+7+5}{3},\ \frac{5+7+7}{3}\right) = \left(\frac{16}{3},\ \frac{19}{3}\right) = (5.333333,\ 6.333333). \]Step 3 — reassign, and check nothing moves.
| Point | to centroid 1 | to centroid 2 | Cluster |
|---|---|---|---|
| A | 0.500000 | 6.871843 | 1 |
| B | 0.500000 | 6.289321 | 1 |
| C | 4.716991 | 1.885618 | 2 |
| D | 8.139410 | 1.795055 | 2 |
| E | 6.946222 | 0.745356 | 2 |
No object changes cluster, so the algorithm has converged in one iteration.
Step 4 — the within-cluster sum of squares. For cluster 1, each point is \(0.5\) from the centroid in \(x\) and \(0\) in \(y\):
\[ W_1 = (0.5)^{2} + (0.5)^{2} = 0.500000. \]For cluster 2, the deviations from \(\left(\tfrac{16}{3}, \tfrac{19}{3}\right)\) are \(C\!\left(-\tfrac43, -\tfrac43\right)\), \(D\!\left(\tfrac53, \tfrac23\right)\), \(E\!\left(-\tfrac13, \tfrac23\right)\), so
\[ W_2 = \frac{16+16}{9} + \frac{25+4}{9} + \frac{1+4}{9} = \frac{32 + 29 + 5}{9} = \frac{66}{9} = 7.333333, \] \[ W = W_1 + W_2 = 0.500000 + 7.333333 = 7.833333. \]Interpretation. The solution matches the two-cluster cut of the dendrogram, which is reassuring but not guaranteed — the two methods optimise different things. Almost all of \(W\) sits in cluster 2, because \(\{C, D, E\}\) is genuinely more spread out than \(\{A, B\}\); \(K\)-means has no objection to unequal spread, but it does implicitly prefer clusters of similar shape, since it uses plain squared distance rather than a Mahalanobis distance.
The input is a matrix of distances \(d_{ij}\) — not the coordinates. The question is whether points in some low-dimensional space could have produced them, and if so, where those points are. The classical solution is four steps of matrix algebra.
Why it works. If the distances really are Euclidean distances between points \(\mathbf{y}_i\) centred at the origin, then \(d_{ij}^{2} = \|\mathbf{y}_i\|^{2} + \|\mathbf{y}_j\|^{2} - 2\mathbf{y}_i'\mathbf{y}_j\), and the double-centring removes the first two terms exactly, leaving \(b_{ij} = \mathbf{y}_i'\mathbf{y}_j\). So \(\mathbf{B} = \mathbf{YY}'\), which is symmetric non-negative definite, and factoring it by its eigen-decomposition recovers \(\mathbf{Y}\) up to a rotation. A rotation is all the ambiguity there can be, since distances determine a configuration only up to rigid motion.
Negative eigenvalues mean the distances are not Euclidean — for instance because they came from a subjective similarity judgement. Small negative values are ignored; large ones say classical scaling is the wrong tool, and non-metric scaling should be used instead. That method keeps only the order of the dissimilarities, and finds a configuration whose distances have the same ranking, by minimising Kruskal's stress
\[ \text{Stress} = \sqrt{\frac{\sum_{i<j}\left(d_{ij} - \hat d_{ij}\right)^{2}} {\sum_{i<j} d_{ij}^{2}}}, \]where \(\hat d_{ij}\) are the fitted monotone values. It is fitted iteratively and has no closed form.
Given. Three objects with \(d_{12} = 3\), \(d_{13} = 5\), \(d_{23} = 4\). No coordinates.
Step 1 — squared distances.
\[ \mathbf{D}^{(2)} = \begin{pmatrix}0 & 9 & 25\\ 9 & 0 & 16\\ 25 & 16 & 0\end{pmatrix}. \]Step 2 — the means. Row means are \(\tfrac{34}{3}\), \(\tfrac{25}{3}\), \(\tfrac{41}{3}\); the matrix is symmetric so the column means are the same; the grand mean is \(\tfrac{100}{9}\).
Step 3 — double-centre. For the \((1,1)\) entry,
\[ b_{11} = -\tfrac12\left(0 - \tfrac{34}{3} - \tfrac{34}{3} + \tfrac{100}{9}\right) = -\tfrac12\left(\frac{-204 + 100}{9}\right) = \frac{52}{9} = 5.777778, \]and for the \((1,3)\) entry,
\[ b_{13} = -\tfrac12\left(25 - \tfrac{34}{3} - \tfrac{41}{3} + \tfrac{100}{9}\right) = -\tfrac12\left(\frac{225 - 102 - 123 + 100}{9}\right) = -\frac{50}{9} = -5.555556. \]Continuing in the same way,
\[ \mathbf{B} = \frac{1}{9}\begin{pmatrix}52 & -2 & -50\\ -2 & 25 & -23\\ -50 & -23 & 73\end{pmatrix}. \]Step 4 — two checks. Every row of \(\mathbf{B}\) sums to zero \(\left(52 - 2 - 50 = 0\right)\), as it must, because the configuration has been centred. And \(\operatorname{tr}\mathbf{B} = (52 + 25 + 73)/9 = 150/9 = 16.666667\), which is the total squared distance from the centroid.
Step 5 — the eigenvalues.
\[ \lambda_1 = 12.964148, \qquad \lambda_2 = 3.702519, \qquad \lambda_3 = 0.000000. \]Their sum is \(16.666667 = \operatorname{tr}\mathbf{B}\). \(\checkmark\) The third is exactly zero because the centring imposed one linear constraint; the first two are positive, so the distances are Euclidean and two dimensions suffice — which is correct, since three points always lie in a plane.
Step 6 — the coordinates. Scaling each eigenvector by \(\sqrt{\lambda_i}\),
| Object | Axis 1 | Axis 2 |
|---|---|---|
| 1 | −2.152311 | −1.070203 |
| 2 | −0.658129 | 1.531223 |
| 3 | 2.810440 | −0.461020 |
Step 7 — check by rebuilding the distances.
\[ d_{12} = \sqrt{(-2.152311 + 0.658129)^{2} + (-1.070203 - 1.531223)^{2}} = 3.000000, \]and likewise \(d_{13} = 5.000000\) and \(d_{23} = 4.000000\). \(\checkmark\)
Interpretation. The recovery is exact, and it had to be: three Euclidean distances can always be realised by three points in a plane. The method earns its keep when \(n\) is large and two axes reproduce \(n(n-1)/2\) distances only approximately — then the proportion \((\lambda_1 + \lambda_2)/\sum_i\lambda_i\) measures how good the picture is, and here it is \(16.666667/16.666667 = 1\).
Factor analysis makes a claim that principal component analysis does not: that the observed correlations are caused by a small number of unobserved variables. The orthogonal \(m\)-factor model for standardised variables is
\[ X_j = l_{j1}F_1 + l_{j2}F_2 + \cdots + l_{jm}F_m + \varepsilon_j, \qquad j = 1, \ldots, p, \]with the factors \(F\) uncorrelated, of mean 0 and variance 1, the specific errors \(\varepsilon\) uncorrelated with each other and with the factors. In matrix form
\[ \boldsymbol\rho = \mathbf{LL}' + \boldsymbol\Psi, \qquad \boldsymbol\Psi = \operatorname{diag}(\psi_1, \ldots, \psi_p). \]Taking the \((j,j)\) entry gives the fundamental decomposition of a variance:
\[ 1 = \underbrace{\sum_{k=1}^{m} l_{jk}^{2}}_{h_j^{2},\ \text{communality}} + \underbrace{\psi_j}_{\text{uniqueness}}, \]and the \((i,j)\) entry with \(i \ne j\) gives the key structural statement:
\[ \rho_{ij} = \sum_{k=1}^{m} l_{ik}l_{jk} \qquad (i \ne j), \]that is, the whole of the correlation between two variables is due to the shared factors. This is a testable claim, and it is the substantive difference from principal components, which are merely a rotation and assume nothing.
Factor loadings are not unique. If \(\mathbf{T}\) is any orthogonal matrix then \((\mathbf{LT})(\mathbf{LT})' = \mathbf{LTT}'\mathbf{L}' = \mathbf{LL}'\), so \(\mathbf{LT}\) fits exactly as well. This indeterminacy is not a defect to be apologised for but a licence: the loadings may be rotated towards an interpretable pattern, and varimax rotation does exactly that by making each variable load heavily on as few factors as possible.
Given. The same correlation matrix again: \(\rho_{12} = 0.6\), \(\rho_{13} = 0.5\), \(\rho_{23} = 0.4\). Asked. Fit a single common factor.
Step 1 — write what the model says. With \(m = 1\) the three off-diagonal equations are
\[ l_1l_2 = 0.6, \qquad l_1l_3 = 0.5, \qquad l_2l_3 = 0.4. \]Three equations in three unknowns — an exactly determined system, which is why \(p = 3\) is the smallest interesting case.
Step 2 — solve it. Multiply the first two and divide by the third:
\[ \frac{(l_1l_2)(l_1l_3)}{l_2l_3} = \frac{l_1^{2}l_2l_3}{l_2l_3} = l_1^{2} = \frac{0.6 \times 0.5}{0.4} = \frac{0.30}{0.40} = 0.750000. \]By the same device,
\[ l_2^{2} = \frac{0.6 \times 0.4}{0.5} = \frac{0.24}{0.50} = 0.480000, \qquad l_3^{2} = \frac{0.5 \times 0.4}{0.6} = \frac{0.20}{0.60} = 0.333333. \]Step 3 — take square roots.
\[ l_1 = 0.866025, \qquad l_2 = 0.692820, \qquad l_3 = 0.577350. \]The common sign is arbitrary: replacing every \(l_j\) by \(-l_j\) leaves every product unchanged, which is the one-dimensional case of the rotational indeterminacy noted above.
Step 4 — reproduce the correlations.
\[ l_1l_2 = 0.866025 \times 0.692820 = 0.6000, \qquad l_1l_3 = 0.866025 \times 0.577350 = 0.5000, \] \[ l_2l_3 = 0.692820 \times 0.577350 = 0.4000. \checkmark \]Four decimals, because the loadings themselves were rounded to six: multiplying two six-decimal figures cannot be trusted in its sixth place.
The fit is exact, because the model had exactly as many parameters as equations.
Step 5 — communalities and uniquenesses.
| Variable | Loading \(l_j\) | Communality \(h_j^{2}\) | Uniqueness \(\psi_j\) |
|---|---|---|---|
| \(X_1\) | 0.866025 | 0.750000 | 0.250000 |
| \(X_2\) | 0.692820 | 0.480000 | 0.520000 |
| \(X_3\) | 0.577350 | 0.333333 | 0.666667 |
Step 6 — the tetrad condition, which is what makes this testable. With four or more variables the one-factor model implies \(\rho_{ij}\rho_{kl} = \rho_{ik}\rho_{jl}\) for all distinct \(i,j,k,l\) — Spearman's tetrad difference vanishes. With only three variables there is no tetrad to check, so a one-factor model always fits three variables exactly and the fit is no evidence at all. That is worth saying plainly: an exact fit here is a consequence of counting, not of the data.
Interpretation, and a comparison worth making. Principal component analysis of this same matrix (Example 4.2) gave a first component with loadings \(0.868,\ 0.821,\ 0.759\) explaining \(66.8\%\) of the variance. The one-factor solution gives \(0.866,\ 0.693,\ 0.577\). They are close for \(X_1\) and diverge for \(X_3\), and the reason is structural: a principal component explains as much total variance as it can, uniquenesses included, whereas a common factor explains only the shared variance and assigns the rest to \(\psi_j\). The sum of communalities here is \(0.750000 + 0.480000 + 0.333333 = 1.563333\), well below \(\lambda_1 = 2.004458\), and the difference is exactly the specific variance that PCA absorbs and factor analysis sets aside.
Path analysis puts a direction on the arrows. A path diagram states which variables are taken as causes of which, and the path coefficients are the standardised regression coefficients of each variable on its stated causes. Its value is the decomposition of a correlation: for a model in which \(X_1\) and \(X_2\) are correlated causes of \(X_3\),
\[ r_{13} = \underbrace{p_{31}}_{\text{direct}} + \underbrace{r_{12}\,p_{32}}_{\text{indirect, through } X_2}, \]and similarly for \(r_{23}\). The coefficients come from the usual two-variable formulae
\[ p_{31} = \frac{r_{13} - r_{12}r_{23}}{1 - r_{12}^{2}}, \qquad p_{32} = \frac{r_{23} - r_{12}r_{13}}{1 - r_{12}^{2}}, \]with residual path \(\sqrt{1 - R^{2}_{3\cdot 12}}\) carrying whatever is unexplained.
What it cannot do. The arrows are assumptions, supplied by the investigator; the data cannot test their direction. Reversing every arrow generally fits the correlation matrix exactly as well. Path analysis quantifies a causal story; it does not establish one.
Given. The same three variables, now read as a model in which \(X_1\) and \(X_2\) are correlated causes of \(X_3\): \(r_{12} = 0.6\), \(r_{13} = 0.5\), \(r_{23} = 0.4\).
Step 1 — the two path coefficients. With \(1 - r_{12}^{2} = 1 - 0.36 = 0.64\),
\[ p_{31} = \frac{0.5 - (0.6)(0.4)}{0.64} = \frac{0.5 - 0.24}{0.64} = \frac{0.26}{0.64} = \frac{13}{32} = 0.406250, \] \[ p_{32} = \frac{0.4 - (0.6)(0.5)}{0.64} = \frac{0.4 - 0.30}{0.64} = \frac{0.10}{0.64} = \frac{5}{32} = 0.156250. \]Step 2 — the explained variance.
\[ R^{2}_{3\cdot 12} = p_{31}r_{13} + p_{32}r_{23} = 0.406250(0.5) + 0.156250(0.4) = 0.203125 + 0.062500 = 0.265625, \]so the residual path is \(\sqrt{1 - 0.265625} = \sqrt{0.734375} = 0.856957\).
Step 3 — decompose \(r_{13}\).
\[ \underbrace{0.406250}_{\text{direct}} + \underbrace{(0.6)(0.156250)}_{\text{through } X_2} = 0.406250 + 0.093750 = 0.500000 = r_{13}. \checkmark \]Step 4 — and \(r_{23}\).
\[ \underbrace{0.156250}_{\text{direct}} + \underbrace{(0.6)(0.406250)}_{\text{through } X_1} = 0.156250 + 0.243750 = 0.400000 = r_{23}. \checkmark \]Interpretation. The two decompositions tell opposite stories from the same matrix. For \(X_1\), \(81\%\) of its correlation with \(X_3\) is direct \((0.40625/0.5)\). For \(X_2\), only \(39\%\) is direct \((0.15625/0.4)\) — most of its apparent association with \(X_3\) is borrowed from \(X_1\), which it happens to be correlated with. That is the same conclusion reached by a quite different route in Example 2.5, where a variable with a respectable simple correlation contributed nothing once its partner was in the model.
Correspondence analysis is principal component analysis for a two-way contingency table of counts. Divide the table by its grand total to get proportions \(p_{ij}\), and form
\[ z_{ij} = \frac{p_{ij} - p_{i\cdot}p_{\cdot j}}{\sqrt{p_{i\cdot}p_{\cdot j}}}, \]the standardised residual from independence — the same quantity whose square, summed and multiplied by \(n\), is Pearson's \(\chi^{2}\). A singular value decomposition of \(\mathbf{Z}\) then places rows and columns on the same plot, so that a row category near a column category is one that occurs together with it more often than independence would predict. The total \(\chi^{2}/n\), called the inertia, is split among the axes exactly as variance is split among principal components.
Conjoint analysis works the other way round: instead of decomposing an observed table, it designs one. Respondents rank or rate a set of hypothetical products, each defined by a combination of attribute levels drawn from a fractional factorial design (the designs of Design and Analysis of Experiments (STS-203)). The ratings are then regressed on attribute-level indicators, and the fitted coefficients — the part-worths — say how much each level contributes. The range of the part-worths within an attribute measures that attribute's importance. It is an ordinary linear model; what makes it multivariate practice is the design of the profiles and the joint reading of the results.