Practical — Multivariate Analysis (STS-205, Section B)

Topics Covered

Section B — The Thirteen Practicals

Practical 1 — MLE of the Mean Vector and Covariance Matrix

1. Problem

Four observations on two variables, \((2, 5), (4, 9), (6, 7), (8, 11)\), are a sample from a bivariate normal population. Find the maximum likelihood estimates \(\hat{\boldsymbol\mu}\) and \(\hat{\boldsymbol\Sigma}\), the unbiased \(\mathbf{S}\), both generalized variances, and the correlation.

2. Aim

To estimate the mean vector and the dispersion matrix of a multivariate normal population by maximum likelihood.

3. Formula

Applying it:

1. Compute the mean vector and the deviations.
2. Form the matrix of sums of squares and products \(\mathbf{A}\).
3. Divide by \(n\) and by \(n-1\); take determinants from the exact values, not the rounded decimals.

4. Calculation

A warning worth taking. Computing \(|\mathbf{S}|\) from the six-decimal entries gives \(44.444449 - 28.444441 = 16.000008\), not 16: a determinant multiplies its entries together, so rounding error multiplies too. Always take a determinant from the exact values.

X <- cbind(c(2, 4, 6, 8), c(5, 9, 7, 11)); n <- nrow(X)
colMeans(X)
A <- crossprod(scale(X, scale = FALSE)); A
A / n                                   # the MLE
c(det_mle = det(A / n), det_S = det(A / (n - 1)), r = cor(X)[1, 2])

[1] 5 8
     [,1] [,2]
[1,]   20   16
[2,]   16   20
     [,1] [,2]
[1,]    5    4
[2,]    4    5
det_mle   det_S       r 
    9.0    16.0     0.8 

5. Result

\(\hat{\boldsymbol\mu} = (5, 8)'\), \(\hat{\boldsymbol\Sigma} = \begin{pmatrix}5 & 4\\ 4 & 5\end{pmatrix}\), \(|\hat{\boldsymbol\Sigma}| = 9\), \(|\mathbf{S}| = 16\) and \(r = 0.8\). The ratio \(9/25 = 0.36 = 1 - r^{2}\) is not a coincidence: for \(p = 2\), \(|\boldsymbol\Sigma| = \sigma_{11}\sigma_{22}(1 - \rho^{2})\), so a generalized variance collapses towards zero as the variables approach an exact linear relation (worked also in Unit 1, Example 1.3).

Practical 2 — Writing the Density, and Reading It Backwards

1. Problem

(a) Write the bivariate normal density with \(\boldsymbol\mu = (5, 8)'\) and \(\boldsymbol\Sigma = \begin{pmatrix}5 & 4\\ 4 & 5\end{pmatrix}\), the estimates of Practical 1. (b) A bivariate normal density is

Find \(\boldsymbol\mu\), \(\boldsymbol\Sigma\), \(\rho\) and the constant \(k\).

2. Aim

To write a multivariate normal density from its parameters, and to read the parameters back out of a given density.

3. Formula

Applying it:

1. (a) Find \(|\boldsymbol\Sigma|\) and \(\boldsymbol\Sigma^{-1}\) (swap the diagonal, negate the off-diagonal, divide by the determinant); expand the quadratic form; assemble.
2. (b) Read \(\boldsymbol\mu\) from the centring; take \(-\tfrac12\) out of the exponent to get \(\boldsymbol\Sigma^{-1}\), halving the printed cross-term; invert; verify.

4. Calculation

(a) \(|\boldsymbol\Sigma| = 25 - 16 = 9\), \(|\boldsymbol\Sigma|^{1/2} = 3\), and \(\boldsymbol\Sigma^{-1} = \frac{1}{9}\begin{pmatrix}5 & -4\\ -4 & 5\end{pmatrix}\). The constant is \(1/(2\pi \times 3) = 1/18.849556 = 0.053052\). With \(u = x_1 - 5\), \(v = x_2 - 8\), \((\mathbf{x}-\boldsymbol\mu)'\boldsymbol\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu) = \frac{1}{9}(5u^{2} - 8uv + 5v^{2})\), the cross-term coefficient being \(2 \times (-4) = -8\) because the matrix contributes \(-4\) twice. So

the 18 being \(2 \times 9\). At \(\mathbf{x} = \boldsymbol\mu\) the exponent is zero, so \(f(5, 8) = 0.053052\), the largest value the density takes. \(\checkmark\)

(b) The exponent is a quadratic form in \(x_1 - 3\) and \(x_2 - 4\), so \(\boldsymbol\mu = (3, 4)'\). Writing \(-\tfrac{1}{12} = -\tfrac12 \times \tfrac16\),

the off-diagonal entry being \(-2\) and not \(-4\), because the printed cross term is the sum of the \((1,2)\) and \((2,1)\) contributions. Its determinant is \((1/6)^{2}(10 - 4) = 1/6\), so \(|\boldsymbol\Sigma| = 6\) and

S <- matrix(c(5, 4, 4, 5), 2)
c(det = det(S), constant = 1 / (2 * pi * sqrt(det(S))))
solve(S) * 9                           # the inverse, times its determinant
Sinv <- matrix(c(2, -2, -2, 5), 2) / 6  # (b): read from the density
Sigma <- solve(Sinv); Sigma
c(rho = Sigma[1, 2] / sqrt(Sigma[1, 1] * Sigma[2, 2]), k = 1 / (2 * pi * sqrt(det(Sigma))))

       det   constant 
9.00000000 0.05305165 
     [,1] [,2]
[1,]    5   -4
[2,]   -4    5
     [,1] [,2]
[1,]    5    2
[2,]    2    2
       rho          k 
0.63245553 0.06497473 

5. Result

(a) \(f(x_1, x_2) = 0.053052\,\exp\{-\frac{1}{18}[5(x_1-5)^{2} - 8(x_1-5)(x_2-8) + 5(x_2-8)^{2}]\}\). (b) \(\boldsymbol\mu = (3, 4)'\), \(\boldsymbol\Sigma = \begin{pmatrix}5 & 2\\ 2 & 2\end{pmatrix}\), \(\rho = 0.6325\), \(k = 0.064975\).

Interpretation. This is the matrix analysed in Practical 7, whose principal components are \(\lambda_1 = 6\) and \(\lambda_2 = 1\): the density's contours have their long axis along \((2, 1)'\), \(\sqrt 6 = 2.449\) times the short one. Every geometric fact about the distribution is in those four numbers.

Practical 3 — Hotelling's \(T^{2}\) for One Sample

1. Problem

\(n = 10\) observations on \(p = 2\) variables have \(\bar{\mathbf{x}} = (7, 12)'\) and \(\mathbf{S} = \begin{pmatrix}4 & 3\\ 3 & 9\end{pmatrix}\). Test \(H_0: \boldsymbol\mu = (5, 9)'\) at 5%.

2. Aim

To test a hypothesis about a mean vector by Hotelling's T-squared.

3. Formula

The test assumes multivariate normality.

Applying it:

1. Form \(\mathbf{d} = \bar{\mathbf{x}} - \boldsymbol\mu_0\) and invert \(\mathbf{S}\).
2. Compute the quadratic form and \(T^{2}\).
3. Convert to \(F\) and compare with \(F_{0.05;\,2,\,8}\).

4. Calculation

Variable by variable, \(t_1 = 2/\sqrt{4/10} = 3.1623\) and \(t_2 = 3/\sqrt{9/10} = 3.1623\), each on 9 d.f., each significant.

xb <- c(7, 12); S <- matrix(c(4, 3, 3, 9), 2); n <- 10; p <- 2
d <- xb - c(5, 9)
T2 <- n * drop(t(d) %*% solve(S) %*% d); F <- (n - p) / ((n - 1) * p) * T2
c(T2 = T2, F = F, F05 = qf(0.95, p, n - p), p = pf(F, p, n - p, lower.tail = FALSE))

         T2           F         F05           p 
13.33333333  5.92592593  4.45897011  0.02637278 

5. Result

\(F = 5.926 > 4.459\) (\(p = 0.026\)): reject \(H_0\); the mean vector is not \((5, 9)'\). Here the separate tests agree. Notice the cancellation \(36 - 36\): because \(x_1\) and \(x_2\) are positively correlated, a joint shift in the same direction is ordinary and counts for little. \(T^{2}\) measures distance in the metric of \(\mathbf{S}^{-1}\) (worked also in Unit 3, Example 3.1).

Practical 4 — Mahalanobis \(D^{2}\) for One Sample

1. Problem

For the data of Practical 3 (\(n = 10\), \(\bar{\mathbf{x}} = (7, 12)'\), \(\mathbf{S} = \begin{pmatrix}4 & 3\\ 3 & 9\end{pmatrix}\), \(\boldsymbol\mu_0 = (5, 9)'\)), find the Mahalanobis distance of the sample mean from \(\boldsymbol\mu_0\), relate it to \(T^{2}\), and compare it with the Euclidean distance.

2. Aim

To measure the distance of a mean vector from a hypothesised value in the metric of the dispersion matrix.

3. Formula

Applying it:

1. Form the deviation vector and \(\mathbf{S}^{-1}\).
2. Compute \(D^{2}\) and \(D\); relate to \(T^{2}\).
3. Compute the Euclidean distance, and \(D^{2}\) for a shift of the same size in the opposite direction.

4. Calculation

\(\mathbf{d} = (2, 3)'\); \(|\mathbf{S}| = 27\), \(\mathbf{S}^{-1} = \tfrac{1}{27}\begin{pmatrix}9 & -3\\ -3 & 4\end{pmatrix}\).

the statistic of Practical 3: \(F = 5.925926\) on \((2, 8)\) d.f., \(p = 0.026373\).

S <- matrix(c(4, 3, 3, 9), 2); d <- c(7, 12) - c(5, 9)
D2 <- drop(t(d) %*% solve(S) %*% d)
c(D2 = D2, D = sqrt(D2), T2 = 10 * D2, euclid = sqrt(sum(d^2)),
  opposite = drop(t(c(2, -3)) %*% solve(S) %*% c(2, -3)))

       D2         D        T2    euclid  opposite 
 1.333333  1.154701 13.333333  3.605551  4.000000 

5. Result

\(D^{2} = 1.3333\), \(D = 1.1547\), and \(T^{2} = nD^{2} = 13.33\), significant at 5% (Practical 3).

Interpretation. The Euclidean distance 3.606 and the Mahalanobis distance 1.155 are not comparable numbers, and the difference is the whole point. The Euclidean figure changes if \(x_1\) is rescaled; the Mahalanobis figure does not. \(D^{2}\) is small because the shift \((2, 3)'\) runs along the direction in which the variables already co-vary; the shift \((2, -3)'\), the same Euclidean size, is three times as far (\(D^{2} = 4\)) in the metric that matters.

Practical 5 — Hotelling's \(T^{2}\) for Equality of Two Mean Vectors

1. Problem

Two samples, \(n_1 = 12\) and \(n_2 = 15\), on \(p = 2\) variables have \(\bar{\mathbf{x}}_1 = (10, 14)'\), \(\bar{\mathbf{x}}_2 = (8, 11)'\) and pooled dispersion matrix \(\mathbf{S}_p = \begin{pmatrix}4 & 2\\ 2 & 5\end{pmatrix}\). Test \(H_0: \boldsymbol\mu_1 = \boldsymbol\mu_2\) at 5%.

2. Aim

To test the equality of two mean vectors by the two-sample Hotelling's T-squared.

3. Formula

The test assumes multivariate normality and equal dispersion matrices.

Applying it:

1. Form \(\mathbf{d} = \bar{\mathbf{x}}_1 - \bar{\mathbf{x}}_2\) and invert \(\mathbf{S}_p\).
2. Compute \(D^{2} = \mathbf{d}'\mathbf{S}_p^{-1}\mathbf{d}\), then \(T^{2}\) with the multiplier \(n_1n_2/(n_1+n_2)\).
3. Convert to \(F\) on \((2, 24)\) d.f. and decide.

4. Calculation

Sp <- matrix(c(4, 2, 2, 5), 2); d <- c(10, 14) - c(8, 11); n1 <- 12; n2 <- 15; p <- 2
D2 <- drop(t(d) %*% solve(Sp) %*% d); T2 <- n1 * n2 / (n1 + n2) * D2
F <- (n1 + n2 - p - 1) / ((n1 + n2 - 2) * p) * T2
c(D2 = D2, T2 = T2, F = F, F05 = qf(0.95, p, n1 + n2 - p - 1), p = pf(F, p, n1 + n2 - p - 1, lower.tail = FALSE))

          D2           T2            F          F05            p 
 2.000000000 13.333333333  6.400000000  3.402826105  0.005920537 

5. Result

\(F = 6.4 > 3.403\) (\(p = 0.0059\)): reject \(H_0\); the two mean vectors differ (worked also in Unit 3, Example 3.2).

Practical 6 — Mahalanobis \(D^{2}\) for Equality of Two Mean Vectors

1. Problem

For the two samples of Practical 5, find the Mahalanobis distance between the two mean vectors, and explain how it differs from \(T^{2}\).

2. Aim

To measure the distance between two population centres in the metric of their common dispersion matrix.

3. Formula

Applying it:

1. Compute \(D^{2}\) from \(\mathbf{d}\) and \(\mathbf{S}_p^{-1}\).
2. Take the square root; relate to \(T^{2}\) through the multiplier.

4. Calculation

Sp <- matrix(c(4, 2, 2, 5), 2); d <- c(10, 14) - c(8, 11)
D2 <- mahalanobis(c(10, 14), c(8, 11), Sp)   # the same quadratic form
c(D2 = D2, D = sqrt(D2), T2 = 12 * 15 / 27 * D2)

       D2         D        T2 
 2.000000  1.414214 13.333333 

5. Result

\(D^{2} = 2\), \(D = 1.414\): the two population centres are about one and a half "standard multivariate units" apart. That is a statement about the populations and does not change with sample size. \(T^{2} = 13.333\) is a statement about the evidence, and it grows with \(n_1\) and \(n_2\) through the multiplier. A large \(T^{2}\) with a small \(D^{2}\) means a reliably detected but practically trivial difference (worked also in Unit 3, Example 3.2, Step 3).

Practical 7 — Computation of Principal Components

1. Problem

(a) Find the principal components of \(\boldsymbol\Sigma = \begin{pmatrix}5 & 2\\ 2 & 2\end{pmatrix}\), the proportion of variance each carries, and the loadings of the first. (b) Find the principal components of the correlation matrix

and decide how many to keep.

2. Aim

To find the principal components of a dispersion or correlation matrix from its eigenvalues and eigenvectors.

3. Formula

\(c_1\) is the trace, \(c_2\) the sum of the three \(2 \times 2\) principal minors and \(c_3\) the determinant. Check every eigenvalue set against the trace and the determinant.

Applying it:

1. Form and solve the characteristic equation.
2. Check: the eigenvalues sum to the trace and multiply to the determinant.
3. Find the normalised eigenvectors; the components, the proportions and the loadings follow.

4. Calculation

(a) \(\operatorname{tr} = 7\), \(\det = 10 - 4 = 6\): \(\lambda^{2} - 7\lambda + 6 = (\lambda - 6)(\lambda - 1) = 0\), so \(\lambda_1 = 6\), \(\lambda_2 = 1\); \(6 + 1 = 7\). \(\checkmark\) From \((\boldsymbol\Sigma - 6\mathbf{I})\mathbf{e} = \mathbf{0}\): \(-e_1 + 2e_2 = 0\), so \(\mathbf{e}_1 = (2, 1)'/\sqrt 5 = (0.894427, 0.447214)'\), and by orthogonality \(\mathbf{e}_2 = (1, -2)'/\sqrt 5 = (0.447214, -0.894427)'\).

(b) \(c_1 = 3\), \(c_2 = 0.64 + 0.75 + 0.84 = 2.23\), \(c_3 = |\mathbf{R}| = 0.47\): \(\lambda^{3} - 3\lambda^{2} + 2.23\lambda - 0.47 = 0\), whose roots, found numerically, are \(2.004458\), \(0.613092\), \(0.382451\) (sum \(3.000000\), product \(0.470000\) \(\checkmark\)).

Kaiser's rule (\(\lambda > 1\)) keeps one component; the cumulative rule (80%) needs two (87.25%).

eigen(matrix(c(5, 2, 2, 2), 2))
R <- matrix(c(1, .6, .5, .6, 1, .4, .5, .4, 1), 3); e <- eigen(R)
round(e$values, 6); c(sum = sum(e$values), prod = prod(e$values))
round(e$vectors, 6)
round(cumsum(e$values) / 3, 6)

eigen() decomposition
$values
[1] 6 1

$vectors
           [,1]       [,2]
[1,] -0.8944272  0.4472136
[2,] -0.4472136 -0.8944272

[1] 2.004458 0.613092 0.382451
 sum prod 
3.00 0.47 
          [,1]      [,2]      [,3]
[1,] -0.613324 -0.179172  0.769240
[2,] -0.579903 -0.559070 -0.592582
[3,] -0.536233  0.809530 -0.238989
[1] 0.668153 0.872516 1.000000

5. Result

(a) \(Y_1\) carries 85.7% of the variance and correlates 0.98 with \(X_1\) and 0.77 with \(X_2\): a "size" component; \(Y_2\), which contrasts them, is a "shape" component. (b) PC1 loads 0.87, 0.82, 0.76 — positively and almost equally on all three variables — and alone reproduces 67% of the correlation structure; PC2 contrasts \(X_3\) with \(X_1\) and \(X_2\). The two rules disagree on how many to keep, which is usual, and why neither is a substitute for looking at the loadings (worked also in Unit 4, Examples 4.1 and 4.2).

Practical 8 — Classification by Discriminant Analysis

1. Problem

For the two populations of Practical 5 (\(\bar{\mathbf{x}}_1 = (10, 14)'\), \(\bar{\mathbf{x}}_2 = (8, 11)'\), \(\mathbf{S}_p = \begin{pmatrix}4 & 2\\ 2 & 5\end{pmatrix}\)), find Fisher's linear discriminant function, the cut-off, the probability of misclassification, and classify \(\mathbf{x}_0 = (9, 12)'\).

2. Aim

To build a linear discriminant function between two normal populations, use it to classify a new observation, and find the probability of misclassification.

3. Formula

Allocate to population 1 if \(y_0 > m\), otherwise to population 2. Check: \(\bar y_1 - \bar y_2 = D^{2}\). The rule assumes normal populations with equal dispersion matrices.

Applying it:

1. Compute the coefficient vector \(\mathbf{a}\).
2. Score the two means; check that their gap is \(D^{2}\).
3. Find the cut-off and the error rate; score and classify \(\mathbf{x}_0\).

4. Calculation

Sp <- matrix(c(4, 2, 2, 5), 2); m1 <- c(10, 14); m2 <- c(8, 11)
a <- solve(Sp, m1 - m2); a
y1 <- sum(a * m1); y2 <- sum(a * m2); y0 <- sum(a * c(9, 12))
c(y1 = y1, y2 = y2, cutoff = (y1 + y2) / 2, y0 = y0, error = pnorm(-sqrt(y1 - y2) / 2))

[1] 0.25 0.50
       y1        y2    cutoff        y0     error 
9.5000000 7.5000000 8.5000000 8.2500000 0.2397501 

5. Result

\(y = 0.25x_1 + 0.5x_2\), cut-off 8.5, misclassification probability 0.240; \(\mathbf{x}_0\) scores 8.25 and is allocated to population 2 — but barely, and honest practice is to report the score rather than only the verdict.

Interpretation. The mean vectors differ beyond reasonable doubt (\(p = 0.0059\), Practical 5), yet nearly a quarter of individuals will be misclassified: a significant difference between group means is not separation between group members. From the table, \(D^{2} = 4\) would be needed to bring the error rate below 16% (worked also in Unit 3, Example 3.5).

Practical 9 — Cluster Analysis by Single, Complete and Average Linkage

1. Problem

Five objects in the plane are \(A(1, 1)\), \(B(2, 1)\), \(C(4, 5)\), \(D(7, 7)\), \(E(5, 7)\). Cluster them hierarchically by single, complete and average linkage, and choose the number of clusters.

2. Aim

To group objects by agglomerative hierarchical clustering, and to see how the linkage changes the merge heights.

3. Formula

Applying it:

1. Compute the distance matrix.
2. Merge the closest pair; recompute the distances to the new cluster by the linkage.
3. Repeat until one cluster remains; cut the dendrogram in its largest gap.

4. Calculation

First merge \(\{A, B\}\) at 1.0000 and then \(\{D, E\}\) at 2.0000, the same for every linkage. The third merge is where the linkages separate:

All three join \(C\) to \(\{D, E\}\). The last merge, \(\{A,B\}\) with \(\{C,D,E\}\), is at 4.4721 (single), 8.4853 (complete) and \(39.686974/6 = 6.6145\) (average, over the six pairs \(AC, AD, AE, BC, BD, BE\)).

P <- rbind(A = c(1, 1), B = c(2, 1), C = c(4, 5), D = c(7, 7), E = c(5, 7))
round(dist(P), 6)
for (m in c("single", "complete", "average")) print(round(hclust(dist(P), m)$height, 4))
cutree(hclust(dist(P), "single"), k = 2)

         A        B        C        D
B 1.000000                           
C 5.000000 4.472136                  
D 8.485281 7.810250 3.605551         
E 7.211103 6.708204 2.236068 2.000000
[1] 1.0000 2.0000 2.2361 4.4721
[1] 1.0000 2.0000 3.6056 8.4853
[1] 1.0000 2.0000 2.9208 6.6145
A B C D E 
1 1 2 2 2 

5. Result

Single-linkage 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\}\). All three linkages give the same tree here and differ only in the heights — single linkage reports the closest pair of members and complete linkage the furthest. With a chain of intermediate points they produce genuinely different trees (worked also, with the dendrogram, in Unit 4, Example 4.4).

Practical 10 — Canonical Variables and Canonical Correlations

1. Problem

Two sets of two standardised variables have

Find the canonical correlations and the canonical variables.

2. Aim

To find the linear combinations of two sets of variables that are most highly correlated with each other.

3. Formula

A matrix \(\begin{pmatrix}\alpha & \beta\\ \beta & \alpha\end{pmatrix}\) has eigenvalues \(\alpha \pm \beta\) with eigenvectors \((1, 1)'\) and \((1, -1)'\).

Applying it:

1. Invert \(\mathbf{R}_{11}\) and \(\mathbf{R}_{22}\).
2. Form \(\mathbf{M}\); its eigenvalues are the squared canonical correlations.
3. Scale each canonical vector to unit variance; verify the correlation directly.

4. Calculation

First pair: \(\mathbf{a} = c(1,1)'\) with \(3c^{2} = 1\), \(c = 0.577350\); \(\mathbf{b} = d(1,1)'\) with \(2.8d^{2} = 1\), \(d = 0.597614\). Check: \(\mathbf{a}'\mathbf{R}_{12}\mathbf{b} = 1.8/\sqrt{8.4} = 0.621059\). \(\checkmark\) Second pair: \(\mathbf{a} = (1, -1)'\), \(\mathbf{b} = 0.912871(1, -1)'\), and \(\operatorname{Corr}(U_2, V_2) = 0.912871(0.3 + 0.3) = 0.547723\). \(\checkmark\)

R11 <- matrix(c(1, .5, .5, 1), 2); R22 <- matrix(c(1, .4, .4, 1), 2)
R12 <- matrix(c(.6, .3, .3, .6), 2)
M <- solve(R11) %*% R12 %*% solve(R22) %*% t(R12); M
sqrt(eigen(M)$values)                   # the canonical correlations

           [,1]       [,2]
[1,] 0.34285714 0.04285714
[2,] 0.04285714 0.34285714
[1] 0.6210590 0.5477226

5. Result

The canonical correlations are \(\rho_1 = 0.621\) and \(\rho_2 = 0.548\), with \(U_1 = 0.5774(X^{(1)}_1 + X^{(1)}_2)\), \(V_1 = 0.5976(X^{(2)}_1 + X^{(2)}_2)\). The strongest link between the two sets is between their sums; the second is between their differences. Neither is visible in the raw \(\mathbf{R}_{12}\): the analysis has replaced four numbers with two, and told you what they mean (worked also in Unit 4, Example 4.3).

Practical 11 — The Orthogonal Factor Model

1. Problem

For three variables with \(\rho_{12} = 0.6\), \(\rho_{13} = 0.5\), \(\rho_{23} = 0.4\), fit the orthogonal factor model with one common factor, and find the communalities and uniquenesses.

2. Aim

To fit a one-factor orthogonal factor model and split each variable's variance into its common and unique parts.

3. Formula

With \(m = 1\) and \(p = 3\) there are three equations in three unknowns, an exactly determined system.

Applying it:

1. Multiply two of the equations and divide by the third, for each loading in turn.
2. Take square roots (the common sign is arbitrary).
3. Reproduce the correlations; compute the communalities and uniquenesses.

4. Calculation

r12 <- .6; r13 <- .5; r23 <- .4
l <- sqrt(c(r12 * r13 / r23, r12 * r23 / r13, r13 * r23 / r12))
round(rbind(loading = l, communality = l^2, uniqueness = 1 - l^2), 6)
round(c(l[1] * l[2], l[1] * l[3], l[2] * l[3]), 4)

                [,1]    [,2]     [,3]
loading     0.866025 0.69282 0.577350
communality 0.750000 0.48000 0.333333
uniqueness  0.250000 0.52000 0.666667
[1] 0.6 0.5 0.4

5. Result

Loadings 0.866, 0.693, 0.577; communalities 0.75, 0.48, 0.33. With only three variables there is no tetrad to check, so a one-factor model always fits three variables exactly: the exact fit is a consequence of counting, not of the data. Compared with PCA of the same matrix (Practical 7: 0.868, 0.821, 0.759), the factor loadings are smaller for \(X_2\) and \(X_3\): a common factor explains only the shared variance (here \(1.563333\), against \(\lambda_1 = 2.004458\)) and assigns the rest to \(\psi_j\) (worked also in Unit 4, Example 4.7).

Practical 12 — Path Coefficients and the Path Diagram

1. Problem

For the same three variables (\(r_{12} = 0.6\), \(r_{13} = 0.5\), \(r_{23} = 0.4\)), read as a model in which \(X_1\) and \(X_2\) are correlated causes of \(X_3\), find the path coefficients and split each correlation with \(X_3\) into its direct and indirect parts.

2. Aim

To estimate the path coefficients of a simple causal model and decompose correlations into direct and indirect effects.

3. Formula

In the path diagram, straight arrows are path coefficients and a double-headed line between \(X_1\) and \(X_2\) is an unanalysed correlation.

Applying it:

1. Compute the two path coefficients.
2. Compute \(R^{2}\) and the residual path.
3. Decompose \(r_{13}\) and \(r_{23}\).

4. Calculation

r12 <- .6; r13 <- .5; r23 <- .4
p31 <- (r13 - r12 * r23) / (1 - r12^2); p32 <- (r23 - r12 * r13) / (1 - r12^2)
R2 <- p31 * r13 + p32 * r23
round(c(p31 = p31, p32 = p32, R2 = R2, residual = sqrt(1 - R2),
        r13 = p31 + r12 * p32, r23 = p32 + r12 * p31), 6)

     p31      p32       R2 residual      r13      r23 
0.406250 0.156250 0.265625 0.856957 0.500000 0.400000 

5. Result

\(p_{31} = 0.406\), \(p_{32} = 0.156\), residual path 0.857. For \(X_1\), 81% of its correlation with \(X_3\) is direct (\(0.40625/0.5\)); for \(X_2\) only 39% is (\(0.15625/0.4\)) — most of its apparent association with \(X_3\) is borrowed from \(X_1\) (worked also, with the path diagram, in Unit 4, Example 4.8).

Practical 13 — Multidimensional Scaling

1. Problem

Three objects have distances \(d_{12} = 3\), \(d_{13} = 5\), \(d_{23} = 4\), and no coordinates. Recover a configuration of points by classical scaling.

2. Aim

To recover a configuration of points from their distances alone, by classical (metric) multidimensional scaling.

3. Formula

Check a scaling solution by rebuilding the distances from the recovered coordinates.

Applying it:

1. Square the distances; find the row, column and grand means.
2. Double-centre to get \(\mathbf{B}\); check that its rows sum to zero.
3. Find the eigenvalues and eigenvectors; scale the eigenvectors by \(\sqrt{\lambda}\) for the coordinates.
4. Rebuild the distances.

4. Calculation

Every row sums to zero \((52 - 2 - 50 = 0)\). The eigenvalues are \(12.964148\), \(3.702519\) and \(0\) (sum \(16.666667\) \(\checkmark\)), and the coordinates are:

and likewise \(d_{13} = 5.000000\) and \(d_{23} = 4.000000\). \(\checkmark\)

D <- matrix(c(0, 3, 5, 3, 0, 4, 5, 4, 0), 3)
fit <- cmdscale(D, k = 2, eig = TRUE)
round(fit$eig, 6)
round(fit$points, 6)
round(dist(fit$points), 6)               # the distances, rebuilt

[1] 12.964148  3.702519  0.000000
          [,1]      [,2]
[1,]  2.152311  1.070203
[2,]  0.658129 -1.531223
[3,] -2.810440  0.461020
  1 2
2 3  
3 5 4

5. Result

Two dimensions reproduce the three distances exactly (\(\lambda_3 = 0\), and the proportion \((\lambda_1 + \lambda_2)/\sum\lambda = 1\)). That had to be so: 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 the \(n(n-1)/2\) distances only approximately; the proportion then measures how good the picture is (worked also in Unit 4, Example 4.6).

Working These by Hand in an Examination

What the Practical Record Should Contain