Past exam of the mathematics course of the University of Cambridge / 2001 / iii / Paper 30 (Experimental Design and Applied Multivariate Analysis)

Let each observation have dimension $p$, and define the sample mean and the divisor-$n$ sample covariance matrix byThe multivariate normal distribution gives the log-likelihoodWriting $y_i-\mu =(y_i-\bar{y})+(\bar{y}-\mu )$ makes the cross terms vanish, and the last sum becomes $n\mathrm{tr}(V^{-1}S)+n(\bar{y}-\mu {)}^T{V}^{-1}(\bar{y}-\mu )$. Since $V$ is a positive-definite matrix, its minimum over the mean is at $\hat{\mu }=\bar{y}$. ConsequentlyWhen $S$ is a positive-definite matrix, $\log \det (V^{-1}S)=\log \det S-\log \det V$. Thuswhere $C^{\prime }$ depends on the data but not on $V$.

For the maximum-likelihood estimator, consider the positive-definite matrix $Q=S^{1/2}{V}^{-1}{S}^{1/2}$, with positive eigenvalues $a_1,\dots ,a_p$. The part of the profile likelihood depending on $V$ is $\frac{n}{2}{\sum }_j(\log a_j-a_j)$. The elementary inequality $\log a-a\le -1$, with equality exactly at $a=1$, shows that the unique maximizing matrix is $Q=I$. HenceThe divisor is $n$, not the unbiased covariance divisor $n-1$. The regularity qualification matters: if $S$ is singular, the displayed logarithm of its determinant is not finite, and no positive-definite covariance maximizes the likelihood. Sending the fitted variance in a null direction of $S$ to zero makes the likelihood unbounded. For a nonsingular normal population, $S$ is positive definite almost surely when $n>p$, which is the setting of the large-sample test.

Under the diagonal covariance restriction, the profile log-likelihood separates into $p$ terms,Differentiating each term gives the restricted maximum-likelihood estimator ${\hat{V}}_0=D=\mathrm{diag}(S_{11},\dots ,S_{pp})$. Both fitted trace terms equal $p$, so the likelihood-ratio test statistic iswhere $R=D^{-1/2}S{D}^{-1/2}$ is the sample correlation matrix. The statistic is nonnegative by the Hadamard inequality; large values indicate dependence. Equivalently, the likelihood ratio is $(\det S/\det D{)}^{n/2}$ and rejection is for small values of this ratio.

For fixed $p$ and a positive-definite true covariance, the Wilks theorem says that twice the maximized log-likelihood difference for nested regular models converges under the null to a chi-squared distribution with degrees of freedom equal to the dimension difference. The unrestricted covariance has $p(p+1)/2$ parameters and the restricted covariance has $p$; the unknown mean has the same $p$ parameters in both. Therefore the Gaussian covariance diagonality likelihood-ratio test of asymptotic level $\alpha$ isWithin the multivariate normal distribution, diagonal covariance also means independent coordinates.

Treat the displayed entries as the lower triangle of a symmetric sample covariance matrix $S$, in the order mechanics, vectors, algebra, analysis, statistics. Compute the precision matrix $K=S^{-1}$. For a multivariate normal distribution, conditioning on all remaining coordinates givesThese identities follow by fixing the other coordinates in the quadratic form of the normal density and completing the square. For the second identity, the conditional precision matrix for the retained pair is its $2\times 2$ principal block; inverting that block gives the negative off-diagonal entry divided by the square root of the two diagonal entries.

Thus the two quantities to evaluate areAn equivalent procedure uses the Schur complement covariance: regress mechanics on the other four variables, or regress analysis and statistics separately on the remaining three. The first residual variance is $S_{11}-S_{1,-1}{S}_{-1,-1}^{-1}{S}_{-1,1}$; the second result is the partial correlation of the two residuals. No arithmetic is required to specify this method.

The model qualification is important: a covariance matrix alone determines linear-regression residual variances and partial correlations. Identifying them with conditional quantities independent of the observed conditioning values uses the multivariate normal model, or another model with the same conditional moments; it does not follow for arbitrary distributions.

Principal component analysis replaces correlated measurements by orthogonal linear coordinates ordered by their variance. For centered observations collected in an $n\times p$ matrix $Z$, write the sample covariance matrix as $S=Z^{T}Z/(n-1)$. A unit loading vector $a$ gives score vector $Za$ with sample variance $a^{T}Sa$. The Rayleigh quotient is maximized by a unit eigenvector of $S$ for its largest eigenvalue.

The subsequent principal components maximize the same quadratic form subject to orthogonality to the earlier loading vectors. By the spectral theorem, if $S=Q\Lambda Q^T$ with descending eigenvalues, the scores $ZQ$ have diagonal covariance $\Lambda$. Their fractions of explained variance are ${\lambda }_j/{\sum }_i{\lambda }_i$. Repeated eigenvalues identify an eigenspace rather than a unique axis.

Keeping the first $q$ principal components gives the rank-$q$ reconstruction $\hat{Z}=Z{Q}_q{Q}_q^T$. Its squared reconstruction error isFor any rank-$q$ orthogonal projection $P$, retained variance is $\mathrm{tr}(SP)={\sum }_j{\lambda }_j{q}_j^{T}P{q}_j$. The weights $q_j^{T}P{q}_j$ lie in $[0,1]$ and sum to $q$, so this is at most the sum of the largest $q$ eigenvalues. Projecting onto any candidate reconstruction subspace is its best least-squares reconstruction; this proves the minimum-error property and agrees with the singular value decomposition. Thus PCA finds the linear subspace preserving the greatest variance, equivalently minimizing squared reconstruction error. The upper-left sketch shows the first axis aligned with the elongated cloud; the second is perpendicular. A scree plot and substantive interpretability help choose $q$.

These methods depend on measurement scale. A variable measured in large numerical units can dominate covariance-based principal component analysis; standardizing variables gives principal component analysis on a correlation matrix. Scores summarize patterns, but the method neither establishes causation nor guarantees that the largest-variance directions best predict a separate response.

Classical multidimensional scaling begins with a symmetric dissimilarity matrix $D$, rather than necessarily with measured coordinates. Its aim is to find a low-dimensional Euclidean configuration representing those dissimilarities. Square the entries, not the matrix product, and defineFor centered coordinates $X$, squared Euclidean distances satisfy $D_{ij}^2=B_{ii}+B_{jj}-2B_{ij}$ with $B=X{X}^T$. Multiplication on both sides by $J$ annihilates the two single-index terms, proving the double-centering formula.

If $B$ is positive semidefinite, use its spectral decomposition $B=U\Lambda U^T$ to constructKeeping all positive eigenvalues reproduces the original Euclidean distances exactly; their number is the minimal embedding dimension. Keeping fewer gives a low-rank approximation to the Gram matrix. It minimizes squared Gram-matrix error, often called strain, rather than necessarily the sum of squared errors in raw distances. Negative eigenvalues signal that the dissimilarities cannot be represented exactly by Euclidean coordinates; discarding them produces an approximation whose adequacy should be inspected.

In the upper-right sketch, the supplied distances come from four points in a rectangle, and the two-dimensional reconstruction recovers that shape. The overall translation, rotation and reflection of the configuration are not identifiable from distances. Classical scaling reconstructs centered geometry from pairwise distances. If the input distances are computed from centered multivariate observations, classical scaling and principal component analysis give the same score configuration up to an orthogonal transformation: their matrices $Z{Z}^T$ and $Z^{T}Z$ share the nonzero squared singular values.

Agglomerative hierarchical clustering starts with each observation as a singleton cluster in cluster analysis, repeatedly joins the closest pair of current clusters, and continues until one remains. The input is a dissimilarity matrix, so choosing variable scales and a scientifically sensible dissimilarity is part of the analysis.

Different linkages define closeness differently:Single-linkage clustering can connect elongated chains through nearest neighbors; complete-linkage clustering emphasizes compact clusters; average-linkage clustering averages all cross-pair distances. For Euclidean observations, Ward minimum-variance clustering instead chooses the smallest increase in within-cluster squared error,This follows by expanding squared deviations around the merged mean and using the zero sums of within-cluster residuals.

A dendrogram plots merge heights. In the lower-left sketch, two close pairs first merge at small heights, then join across a much larger separation; a horizontal cut between these heights yields two groups. The hierarchy gives nested partitions, not a uniquely established number of populations. Examine sensitivity to scaling, linkage and ties; once an agglomerative merge is made, the method does not undo it. The absence of a probability model also means that a visually attractive dendrogram is not itself a significance test.

Multivariate analysis of variance tests differences between group mean vectors when each experimental unit has several correlated responses. For independent $p$-dimensional observations in $g$ groups, use a common positive-definite covariance matrix $\Sigma$ and the multivariate normal model $Y_{ri}\sim N_p({\mu }_r,\Sigma )$. The null hypothesis is equality of all $g$ means.

Write $n={\sum }_r{n}_r$, group means ${\bar{y}}_r$, and the size-weighted grand mean $\bar{y}$. The within-group and between-group scatter matrices areExpanding deviations about the group means yields total scatter $T=E+H$, with zero cross terms. Under the null, orthogonal normal projections give independent Wishart distributions $E\sim W_p(\Sigma ,n-g)$ and $H\sim W_p(\Sigma ,g-1)$.

With nonsingular $E$, the Wilks lambda statistic iswhere the nonnegative ${\theta }_j$ are eigenvalues of $E^{-1/2}H{E}^{-1/2}$. Maximizing the common covariance under the two mean models gives the likelihood ratio ${\Lambda }^{n/2}$, so small $\Lambda$ indicates large group separation relative to within-group spread. For fixed $p,g$ and large group sizes, the Wilks theorem gives $-n\log \Lambda \Rightarrow {\chi }_{p(g-1)}^2$ under the null; a finite-sample reference distribution can also be used.

The lower-right sketch shows two elongated groups whose centroid difference lies nearly across, rather than along, their long axes. Their correlation makes this multivariate separation much clearer than either marginal mean shift alone. Nonsingular linear rescaling changes $E$ and $H$ by congruence, leaving the determinant ratio unchanged. MANOVA measures mean separation against the full within-group covariance geometry. Its interpretation requires independent observations and an adequate common-covariance model; separate coordinatewise tests discard that geometry and raise multiple-testing issues.

A balanced incomplete block design has $t$ treatments and $b$ blocks, each containing $k$ distinct treatments with $1<k<t$. Every treatment occurs in $r$ blocks and every distinct treatment pair occurs together in $\lambda >0$ blocks. A symmetric balanced incomplete block design has $b=t$, which implies $r=k$ by the first counting identity below.

For the separate even-order condition, let $N$ be the square treatment-by-block incidence matrix of a set system of a symmetric design. Its diagonal overlap counts are $r$ and its off-diagonal counts are $\lambda$, givingThe counting identities proved in the next slots show that $r=k$ and $r+(t-1)\lambda =r^2$. Consequently the eigenvalues of $N{N}^T$ are $r^2$ on the constant vector and $r-\lambda$ on its $(t-1)$-dimensional orthogonal complement. Taking determinants givesIf $b=t$ is even, the exponent $t-1$ is odd. The exponent of every prime in $\det (N{)}^2$ and in $r^2$ is even; therefore its exponent in the positive integer $r-\lambda$ must be even. This proves the even-order symmetric design square obstruction:It is a necessary condition, not a sufficiency claim.

Count incidences of a treatment with a block containing it. Summing over blocks gives $bk$, while summing over treatments gives $tr$. Hence the balanced incomplete block design satisfies

Fix one treatment. Each of its $r$ blocks supplies $k-1$ companion treatments, giving $r(k-1)$ incidences. Alternatively, each of the other $t-1$ treatments accompanies it in exactly $\lambda$ blocks. Thus

Let $N$ be the $t\times b$ incidence matrix of a set system. The balanced incomplete block design overlap counts give $N{N}^T=(r-\lambda )I+\lambda J$. On the zero-sum subspace its eigenvalue is $r-\lambda =r(t-k)/(t-1)>0$; on the constant vector its eigenvalue is $r+(t-1)\lambda >0$. Therefore $N{N}^T$ has matrix rank $t$. Since a product cannot have larger rank than either factor,This proves Fisher's inequality for block designs directly. The incomplete, nontrivial hypotheses ensure the positive eigenvalue needed for this argument.

Symmetry forces $r=k$, and the pair-counting identity becomes $k(k-1)=2\cdot 45=90$. Its only positive integral solution is $k=r=10$. Since the number of blocks is even, the even-order symmetric design square obstruction requires $r-\lambda =8$ to be a square. It is not. No design with these parameters exists.

The identities require $r=k$ and $k(k-1)=2\cdot 3=6$, so $k=r=3$. Take treatments $1,2,3,4$ and the four blocksEach treatment is absent from exactly one block, hence occurs three times. Each pair is absent from the two blocks missing one of its members, hence occurs together in the other two. This is the required symmetric balanced incomplete block design; here $r-\lambda =1$ is indeed a square.

Code the two levels of each factor as $\pm 1$. For a factor subset $A$, its factorial contrast column is the character ${\chi }_A(x)={\prod }_{j\in A}{x}_j$. Products of columns correspond to symmetric differences of subsets, because $x_j^2=1$.

Let $G_1,\dots ,G_k$ be defining words whose incidence vectors are linearly independent over ${\mathbb{F}}_2$. Write $x_j=(-1{)}^{u_j}$, with $u_j\in {\mathbb{F}}_2$. Specifying signs ${\chi }_{G_i}(x)=s_i$ is then a rank-$k$ system of binary linear equations. Every one of its $2^k$ right-hand sides has $2^{n-k}$ solutions. Thus $k$ independent defining contrasts split the full $2^n$ design into $2^k$ equal fractions. Independence is necessary: dependent defining words do not yield $2^k$ distinct sets.

The defining contrast subgroup has $2^k$ words. On a selected fraction its word $G$ has a constant sign $s_G$, soHence the coefficients of $A$ and $AG$ cannot be estimated separately without assumptions or extra runs: this is aliasing in a fractional factorial design. Conversely, averaging a character over the affine binary solution space vanishes unless its word is in the defining subgroup, so these signed subgroup cosets give all aliases. Here a lowercase treatment label records the factors at their high levels, and $(1)$ puts all five at their low levels; signs must be retained when listing aliases.

On these eight runs, $W=M$ and $SCT=-1$. The two independent defining words therefore give the full signed defining relationMultiplying each desired factorial contrast by its four words gives all aliases. In the following equalities a minus sign means opposite contrast columns:Thus two main effects, $W$ and $M$, are inseparable even when every unwanted interaction vanishes. In addition, the important interaction $ST$ aliases with the main effect $C$. Design (a) is unsuitable. The seven-column design matrix for the intercept, the five main effects and $ST$ has rank five, confirming the two independent estimability failures.

For this fraction, $SCW=-1$ and $SCMT=1$, so their product is $WMT=-1$. The signed defining relation isThe complete factorial contrast alias classes areThe five main effects and $ST$ occupy six distinct alias classes, and none aliases with the intercept. Every undesired two-factor term in these classes is assumed negligible, and all terms of order at least three are also assumed negligible. Thus design (b) is suitable under the stated interaction assumptions. Its intercept and six requested columns are mutually orthogonal, with $X^{T}X=8I_7$, so each coded regression coefficient is ${\hat{\beta }}_A=8^{-1}{\sum }_{i=1}^8{\chi }_A(x_i)Y_i$. The corresponding high-minus-low main effect is $2{\hat{\beta }}_A$.

There is one residual degree of freedom, but no replicated setting and therefore no separate pure error estimate. Suitability here means estimability under the stated model; it does not protect estimates from active interactions that the assumptions exclude.

For the general normal linear model, minimizing $\Vert Y-X\beta {\Vert }^2$ gives the least-squares normal equations $X^{T}X\hat{\beta }=X^{T}Y$. Full column rank makes $X^{T}X$ invertible, soThe mean and covariance follow by substituting $Y=X\beta +ϵ$ and using $\mathbb{E}ϵ=0$, $\mathrm{Cov}(ϵ)={\sigma }^{2}I$. Normal errors also give the exact multivariate normal distribution of this estimator.

Use coded experimental variablesTheir origin is the middle of the allowed rectangle, and the stated endpoints become $\pm 1$. The four corners and $n_0$ center runs give mutually orthogonal columns $1,x_1,x_2$, withWriting ${\bar{Y}}_0=n_0^{-1}{\sum }_{i=5}^{4+n_0}{Y}_i$ when $n_0>0$, the ordinary least squares estimates areTheir variances are ${\sigma }^2/(4+n_0)$, ${\sigma }^2/4$, ${\sigma }^2/4$, with zero covariances.

The repeated center runs estimate pure error without requiring the fitted mean surface to be correct:They also permit a center-point curvature contrast. Let ${\bar{Y}}_F=(Y_1+Y_2+Y_3+Y_4)/4$ and $\hat{\gamma }={\bar{Y}}_F-{\bar{Y}}_0$. Under the first-order model, $\mathbb{E}\hat{\gamma }=0$ and $\mathrm{Var}(\hat{\gamma })={\sigma }^2(1/4+1/n_0)$. Under a quadratic surface, its expectation is ${\beta }_{11}+{\beta }_{22}$, since the linear and $x_1{x}_2$ terms average to zero. Thusunder the first-order normal model. With only one center run there is no pure-error degree of freedom. A zero center contrast does not exclude all curvature, because the two quadratic coefficients can cancel. The corners also estimate an interaction contrast; the full Lack-of-fit F-test compares the first-order fit with the five setting means, with two lack-of-fit and $n_0-1$ pure-error degrees of freedom.

If the first-order model is adequate and $\hat{b}=({\hat{\beta }}_1,{\hat{\beta }}_2{)}^T\ne 0$, steepest ascent in response surface methodology uses the unit coded direction $d=\hat{b}/\Vert \hat{b}\Vert$. The Cauchy-Schwarz inequality proves that it maximizes the predicted directional derivative. Run sequential settings $x=hd$ for increasing feasible $h$, corresponding to physical increments $(5h{d}_1,2.5h{d}_2)$ from the center. Stop or reduce the steps when measured yield ceases to improve, and refit locally near promising settings. Respect the allowed rectangle; if a boundary is reached, investigate feasible boundary directions rather than extrapolating outside it. A near-zero gradient offers no first-order ascent direction and calls for further modeling.

If curvature is indicated, add the four axial settings $(1,0),(-1,0),(0,1),(0,-1)$, together with further center replication. This face-centered central composite design stays inside the allowed physical ranges and separates the two squared coordinates, which the corners alone cannot do. Fit the full quadratic regression surfaceWrite its fitted quadratic part as $x^{T}Bx$, where $B_{11}={\hat{\beta }}_{11}$, $B_{22}={\hat{\beta }}_{22}$ and $B_{12}=B_{21}={\hat{\beta }}_{12}/2$. Its gradient and Hessian are $\hat{b}+2Bx$ and $2B$. If $B$ is nonsingular, the stationary setting isIf $B$ is negative definite and $x_{\ast }$ is feasible, it is the fitted maximum. Otherwise compare feasible stationary settings, optima on each edge obtained from the restricted one-variable quadratic, and the corners. A singular or poorly determined $B$ calls for attention to ridge directions and more runs. Convert any chosen coded setting back using ${\xi }_1=95+5x_1$, ${\xi }_2=37.5+2.5x_2$, and confirm its yield experimentally. The fitted local surface guides the search; it does not establish an untested physical optimum.