Past exam of the mathematics course of the University of Cambridge / 2004 / iii / Paper 43 (Applied Multivariate Analysis)

Write the expected value as ${\mu }^T=({\mu }_1^T,{\mu }_2^T)$ and assume the relevant covariance matrix $V_{22}$ is invertible. Remove the linear predictor from the first block by definingA linear image of a multivariate normal vector is again multivariate normal. Direct covariance calculation givesThus $Y$ and $X_2$ are independent: uncorrelated jointly Gaussian variables are independent. Expanding the covariance matrix of $Y$, using $V_{21}=V_{12}^T$, givesConditioning on $X_2=x_2$ translates $Y$ by the fixed vector ${\mu }_1+V_{12}{V}_{22}^{-1}(x_2-{\mu }_2)$ and leaves its covariance matrix unchanged. Consequently the conditional multivariate normal distribution isThe conditional covariance is the Schur complement and does not depend on $x_2$.

Put $K=V^{-1}$, the precision matrix, with $X_1$ scalar and $X_2$ containing all remaining coordinates. Holding $X_2=x_2$ in the multivariate normal density, the terms involving $x_1$ in its exponent areCompleting the square gives a normal distribution with mean ${\mu }_1-K_{11}^{-1}{K}_{12}(x_2-{\mu }_2)$ andEquivalently, the block matrix inverse identity gives $K_{11}=(V_{11}-V_{12}{V}_{22}^{-1}{V}_{21}{)}^{-1}$. The Gaussian conditional variance from precision is a reciprocal diagonal entry of the precision matrix, rather than the marginal variance $V_{11}$.

Write ${\sigma }_i^2=V_{ii}$ and $V_{ij}={\rho }_{ij}{\sigma }_i{\sigma }_j$. Apply the Schur complement covariance to the first two coordinates, conditioning on the third. The conditional covariance and conditional variances areDivide the conditional covariance by the product of the conditional standard deviations. The conditional correlation is thereforeThis is the partial correlation after removing the linear effect of $X_3$. For a nonsingular multivariate normal distribution, both residual variances are positive and the answer is independent of $x_3$. If a residual variance vanishes, the requested correlation coefficient is undefined, even though the residual covariance formula remains valid.

Use the usual sample mean and the divisor-$n-1$ sample covariance matrix,The residual sum is zero, so expanding $x_j-\mu =(x_j-\bar{x})+(\bar{x}-\mu )$ givesMultiplying the independent multivariate normal densities and taking logarithms yields the log-likelihoodHere $|V|$ is the determinant and $\mathrm{tr}$ denotes the trace. The maximum-likelihood estimators areIf $S$ instead denotes the divisor-$n$ sample covariance matrix, replace $(n-1)S$ by $nS$ and then $\hat{V}=S$. These formulas describe a nonsingular maximum-likelihood estimator when the empirical scatter is positive-definite; if it is singular, there is no maximizer over positive-definite covariance matrices.

Let $n={\sum }_{\nu =1}^g{n}_{\nu }$, ${\bar{x}}_{\nu }=n_{\nu }^{-1}{\sum }_j{x}_j^{(\nu )}$ and $\bar{x}=n^{-1}{\sum }_{\nu }{n}_{\nu }{\bar{x}}_{\nu }$. Define the within-group and between-group scatter matrices byThe within-group and between-group scatter decomposition gives total scatter $T=W+B$: in each group the residuals about ${\bar{x}}_{\nu }$ sum to zero, so all cross terms vanish. Under unrestricted group means, the maximum-likelihood estimators are ${\hat{\mu }}_{\nu }={\bar{x}}_{\nu }$ and $\hat{V}=W/n$. Under the common-mean null hypothesis, they are $\hat{\mu }=\bar{x}$ and ${\hat{V}}_0=(W+B)/n$.

For either fitted scatter matrix $A$, substitution into the log-likelihood givessince $\mathrm{tr}[(A/n{)}^{-1}A]=np$. Thus the maximized likelihood ratio isA generalized likelihood-ratio test rejects for small ${\Lambda }_{\mathrm{LR}}$, equivalentlyThe cutoff $c$ is chosen for the desired significance level. The determinant ratio $|W|/|W+B|$ itself is the Wilks lambda statistic. The derivation assumes $W$ is positive-definite, which for continuous normal samples requires enough within-group degrees of freedom, $n-g\ge p$; singular $W$ does not support this ordinary determinant-based test without a modified model.

The within-group and between-group scatter decomposition separates noise from group separation: $W$ measures variation about each group sample mean, and $B$ measures separation of those means about the overall sample mean. In linear discriminant analysis, seek a projection $a^{T}x$ maximizing between-group relative to within-group variation,Normalize by $a^{T}Wa=1$; a Lagrange multiplier calculation gives the generalized eigenvalue problem $Ba=\lambda Wa$. The leading canonical discriminant directions have the largest eigenvalues, and subsequent directions are orthogonal in the $W$ inner product. There are at most $\min (p,g-1)$ nonzero separating directions. Projecting observations onto these directions gives a low-dimensional display of group separation.

For classification, estimate the common covariance matrix by $S_p=W/(n-g)$ and use the linear discriminant analysis scoreswhere ${\pi }_{\nu }$ are class prior probabilities. Assign to the largest score. The same within-group scatter sets the noise scale, while between-group scatter selects the directions in which the group means separate most strongly.

The first command reads a rectangular data file into a; header=T treats the first line as variable names. There are twelve observations, identified by student names, and nine binary attributes. A one records presence of an attribute and a zero records absence. This is a description of the encoding, not a claim that every presence is equally important for a scientific application.

dist(a, metric="binary") compares rows using Jaccard distance. If $n_{11}$ counts shared presences and $n_{10},n_{01}$ count mismatched presences, thenJoint absences do not enter the denominator, unlike the simple matching coefficient or normalized Hamming distance. dist2full(d) expands the stored pairwise distances into a symmetric dissimilarity matrix; round(...,2) changes only its displayed precision. The clustering command still receives the original unrounded distances.

For example, Tom and Martin have six shared presences and only one mismatch, giving $d=1/7=0.142857\dots$. Taeko and Luitgard have three shared presences and four mismatches, giving $d=4/7=0.571428\dots$. LingChen and Mohammad share only the meat attribute and differ on four others, giving $d=4/5=0.8$. The diagonal is zero, and the entries are symmetric because exchanging the two observations exchanges $n_{10}$ and $n_{01}$.

hclust(d, method="compact") constructs an agglomerative hierarchical clustering object by complete-linkage clustering. Its distance between two current groups isStart with singleton groups, merge a pair with smallest intergroup distance, update these maximum cross-pair distances, and repeat until one group remains. The object stores the merge structure, merge heights and a leaf order suitable for a dendrogram; assigning it to h does not itself display a graph. The historical S-Plus help identifies compact with complete linkage and describes the plotting object.

The first merge is Tom with Martin at $1/7$. At the next height $1/5$, either LinYee with HuiChin or HuiChin with Meg may merge: the distances are tied. The following dendrogram shows one valid resolution, using original row order to resolve ties. In this tree, LinYee joins HuiChin at $1/5$, Pio joins Meg at $1/4$, and Nicolas joins Tom–Martin at $3/8$. The latter height is $\max \{d(\text{Nicolas},\text{Tom}),d(\text{Nicolas},\text{Martin})\}=\max \{3/8,2/7\}$, illustrating the complete-linkage update.

For the displayed tree, a cut at height $0.65$ gives three groups: Taeko–Luitgard; Alet–Tom–Martin–Nicolas–Mohammad; and LinYee–Pio–LingChen–HuiChin–Meg. The last two merges occur at $5/7$ and $4/5$. The hierarchy is a description under the chosen binary distance and linkage rule; it neither supplies an intrinsically correct number of groups nor guarantees a unique tree in the presence of ties. Changing the treatment of joint absences, or the weighting of attributes, can change that description.

princomp(): principal component analysis. Supply observations in rows and numeric variables in columns. After centering, let $Z$ denote the resulting matrix and $C=Z^{T}Z/n$ its sample covariance matrix, with an explicitly stated normalization. WriteThe first direction maximizes the Rayleigh quotient $q^{T}Cq$ over unit vectors; successive directions maximize it subject to orthogonality to the earlier directions. The principal component scores are $z_{ij}=q_j^T(x_i-\bar{x})$, with variances ${\lambda }_j$ and zero pairwise sample covariance. Zero covariance alone does not imply independence. Each principal component loading describes a direction in variable space; its sign is arbitrary, and repeated eigenvalues leave the basis inside the corresponding eigenspace arbitrary.

The function returns principal component loadings, component standard deviations sdev and, when requested, principal component scores. By default it uses a covariance matrix; cor=T instead performs principal component analysis on a correlation matrix, standardizing nonconstant variables to variance one. This choice matters when variables have different measurement units. These interface details are documented in the S-Plus help for princomp.

A scree plot shows the ordered eigenvalues; a plot of the first two principal component scores can reveal groups, unusual observations or a dominant gradient. The first $k$ components explain ${\sum }_{j\le k}{\lambda }_j/{\sum }_j{\lambda }_j$ of the total variance. Reconstructing each centered observation by its orthogonal projection onto these directions minimizes average squared reconstruction error among $k$-dimensional linear subspaces. Thus PCA gives an optimal linear compression for squared reconstruction error, with scaling and the retained dimension chosen explicitly. It is an unsupervised learning method, rather than a classification rule fitted to supplied labels.

tree(): classification trees and regression trees. A call such as tree(response ~ predictors, data=dat) fits a recursive partitioning model. At a current node, test numerical predictors against candidate thresholds, or divide the levels of a categorical predictor into two sets. Select a split that most improves the fitting criterion, then repeat within its children. Numerical splits give axis-aligned regions. A categorical response produces a classification tree; a numeric response produces a regression tree, as specified in the S-Plus tree documentation.

For a regression tree, the constant prediction in a terminal region is its sample mean, minimizing the sum of squared residuals there. For a classification tree, estimate class probabilities by the observed proportions in a leaf and predict a most probable class under equal misclassification costs. The classification-tree deviance of a node iswith $0\log 0=0$; splitting maximizes the decrease from parent to children. This fitting criterion differs from the number of wrong classifications. The middle sketches show the same three-region partition as a binary tree: first split on variable one, then split the right-hand branch on variable two.

Growing to small leaves can overfit. Use cost-complexity tree pruning, balancing fitted loss $R(T)$ against leaf count,and select complexity by cross-validation. Plotting the fitted tree and annotating nodes displays the rules; prediction routes each new observation through the fitted splits. A tree gives an interpretable, piecewise-constant prediction rule; training fit alone does not measure prediction accuracy.

The rpart() alternative. This also fits classification and regression trees by recursive partitioning, with squared-error splitting for ordinary numeric responses and, by default for classification, Gini impurity $1-{\sum }_r{p}_r^2$. Classification settings can instead use information-based splitting and specify prior probabilities and a loss matrix. It supports surrogate splits to approximate a node's primary left/right assignment when the primary predictor is missing. Its complexity-parameter path and cross-validated errors support cost-complexity tree pruning; for example, select a complexity parameter and apply prune(fit, cp=...). The maintainer's rpart documentation describes the criteria and the distinction in handling surrogate predictors. The partition and tree sketches apply to either implementation; the fitted tree need not be the same because splitting criteria and controls can differ.

cmdscale(): classical multidimensional scaling. Its input is a dissimilarity matrix or stored distance object. It outputs coordinates for the observations in a chosen dimension, usually two for visualization. Let $D^{(2)}$ contain the squared entries of $D$, and putFor centered Euclidean coordinates $X$, the distance identity $d_{ij}^2=\Vert x_i{\Vert }^2+\Vert x_j{\Vert }^2-2x_i^T{x}_j$ gives $G=X{X}^T$: multiplying by $H$ removes the row and column norm terms. Diagonalize this Gram matrix and retain its $k$ largest positive eigenvalues. The coordinates areThus classical scaling exactly recovers Euclidean distances when all positive eigenvalues are retained. Reducing dimension approximates the centered Gram matrix; it does not in general minimize squared discrepancies in the original distances. Negative eigenvalues identify a failure of the Euclidean distance matrix criterion, so their size matters when interpreting a display. The historical S-Plus cmdscale help describes the dimension argument, eigenvalue output and optional additive adjustment to dissimilarities.

Plot the two coordinate columns on equally scaled axes and label observations. Nearby plotted points suggest small dissimilarities only insofar as the dimension reduction fits the input well. The bottom sketches give an exact Euclidean example: a configuration and its independently recovered two-dimensional coordinates have identical pairwise distances, despite different coordinate axes. Classical scaling converts distances into geometry; translations, rotations and reflections do not change the represented distances. For Euclidean distances between observations, its coordinates agree with the corresponding principal component scores up to these harmless choices of axes.

Under equal costs for either wrong assignment, the Bayes classifier chooses the class of greatest posterior probability. With positive prior probabilities, it assigns class one when ${\pi }_1{f}_1(x)>{\pi }_2{f}_2(x)$; ties can be resolved arbitrarily. For the two multivariate normal densities with common positive-definite covariance matrix $V$, their normalizing factors cancel. The log posterior odds areTherefore the Gaussian Bayes classifier has the required linear form, withA larger prior probability for class one decreases $b$, expanding its decision region. The prior-dependent Gaussian discriminant boundary is the hyperplane $a^{T}x=b$.

For equal prior probabilities, $b=\frac{1}{2}{a}^T({\mu }_1+{\mu }_2)$. The projected score is normal under either class because it is a linear image of a multivariate normal vector, and its variance isMoreover $a^T{\mu }_1-a^T{\mu }_2={\delta }^2$, so the midpoint threshold lies ${\delta }^2/2$ from either projected mean. HenceStandardizing the two normal distributions, with $\Phi$ the standard normal distribution function, givesThus the equal-covariance Gaussian classification error is $p(\delta )=\Phi (-\delta /2)$. The Mahalanobis distance $\delta$ between the class means is their separation in units of within-class noise. It yields $p\to 1/2$ as $\delta \downarrow 0$ and $p\to 0$ as $\delta \to \infty$. The stipulated $\delta >0$ means the two class means are distinct.