Past exam of the mathematics course of the University of Cambridge / 2002 / iii / Paper 39 (Experimental Design and Multivariate Analysis)

Work with a nonsingular multivariate normal distribution, so $V$ is positive-definite. The claimed finite maximum-likelihood estimates also require the centered sample covariance matrix to be nonsingular; the exceptional case is discussed below.

Independence gives, up to a constant, the log-likelihoodLet $\bar{Y}=n^{-1}{\sum }_i{Y}_i$ and $S=n^{-1}{\sum }_i(Y_i-\bar{Y})(Y_i-\bar{Y}{)}^T$. The centered sum is zero, and expansion therefore givesFor every fixed positive-definite matrix $V$, the extra quadratic term in the log-likelihood is $-n(\bar{Y}-\mu {)}^T{V}^{-1}(\bar{Y}-\mu )/2$, with unique maximum at $\mu =\bar{Y}$.

Assume $S$ is positive-definite. Maximizing over $V$ now means minimizing $\log \det V+\mathrm{tr}(V^{-1}S)$. Put $A=S^{1/2}{V}^{-1}{S}^{1/2}$, another positive-definite matrix. ThenIf $a_1,\dots ,a_p>0$ are its eigenvalues, the variable part is ${\sum }_j(a_j-\log a_j)$. The scalar inequality $a-\log a\ge 1$, with equality only at $a=1$, shows that this is minimized uniquely at $A=I$, or $V=S$. ThusThe divisor is $n$ because this is maximum likelihood estimation, not the unbiased covariance estimate with divisor $n-1$.

The nonsingularity condition matters. The centered scatter has rank at most $n-1$. If $S$ is singular, choose $\mu =\bar{Y}$, retain positive covariance eigenvalues on its range and let a covariance eigenvalue on its null space tend to zero. The quadratic likelihood term remains finite, while $-\log \det V\to +\infty$. This is singular covariance and nonexistence of a Gaussian maximum likelihood estimate: there is no finite maximizer over nonsingular $V$. Under a nonsingular Gaussian population, $S$ is positive definite almost surely when $n>p$; for $n\le p$ it is necessarily singular. The asserted estimate and the inverse in the next part have their usual meaning under $n>p$.

By the spectral theorem, write $V=U\mathrm{diag}(v_1,\dots ,v_p)U^T$, where $U$ is orthogonal and each $v_j>0$. Choose the whitening transformationA linear image of a multivariate normal vector is multivariate normal, and $LV{L}^T=I$. Thus $Z_i=L(Y_i-\mu )$ are independent $N_p(0,I)$ random vectors. Their whole joint distribution is independent of $\mu ,V$.

Define $\bar{Z}=n^{-1}{\sum }_i{Z}_i$ and $S_Z=n^{-1}{\sum }_i(Z_i-\bar{Z})(Z_i-\bar{Z}{)}^T$. These satisfyThereforeThe right side is a function of independent standard multivariate normal observations only, so its distribution depends on $n,p$ but not on $\mu ,V$. This is the same affine invariance of Hotelling's statistic that permits multivariate mean tests with unknown covariance. The proof concerns the jointly transformed mean and covariance; whitening the mean alone would not establish the conclusion. Inverses exist almost surely under $n>p$. If $V$ is singular, no square matrix can whiten it to $I$, so nonsingularity of the population covariance is necessary as well.

The calculation is principal component analysis on a correlation matrix, rather than on an unscaled covariance matrix. For each of the eleven rating variables, subtract its column mean and divide by its column standard deviation. If $Z$ is the resulting $43\times 11$ matrix, use the divisor convention of the fitted object to writeHere the scaling standard deviations use divisor $43$, consistently with this convention. Using sample standard deviations with divisor $42$ gives the same correlation matrix and directions, with the corresponding score scale. Every diagonal entry of $R$ is one, so ${\sum }_k{\lambda }_k=\mathrm{tr}R=11$.

A principal component is the linear combination of the standardized ratings along $a_k$. Its principal component score for judge $i$ is $t_{ik}=z_i^T{a}_k$, and its variance is ${\lambda }_k$. The displayed component standard deviations are $\sqrt{{\lambda }_k}$; the explained variance of a principal component is ${\lambda }_k/11$. In particular,The first component explains about $92.122\%$ of the standardized variation, and the first two jointly explain $95.982\%$. Components three, four and five add about $2.321\%$, $0.830\%$ and $0.339\%$. The cumulative line adds these proportions successively; it is not the proportion explained by that component alone. Retaining two components discards about $4.018\%$ of the total standardized squared variation.

The displayed principal component loadings are the unit eigenvector coefficients $a_{jk}$. The first direction has positive, broadly similar weights for all ratings. It measures an overall favorable-rating tendency: a judge high on most ratings has a high first score. The second direction compares integrity and demeanour, with large positive coefficients, against such attributes as decision promptness, case-flow management and physical ability, with negative coefficients. Thus among similarly rated judges, a positive second score indicates relatively stronger integrity/demeanour compared with those operational ratings. This is a description of a data direction, not a causal interpretation. Multiplying a whole eigenvector and its scores by minus one reflects the display but changes nothing substantive.

These principal component loadings are not themselves variable-component correlations. With standardized variables,For example, the integrity correlation with component one is approximately $3.1833\times 0.289\simeq 0.920$. Blank loadings in the printed second column indicate suppressed small entries, not exact zeros. The complete unrounded loading vectors in the fitted object must be used to calculate scores.

A two-dimensional representation plots the 43 pairswith each point labeled by its judge. In the fitted software object these are the first two columns of the score matrix. Nearby points have similar projections of their standardized ratings. The rank-two reconstruction ${\hat{z}}_i=t_{i1}{a}_1+t_{i2}{a}_2$ is the orthogonal projection onto the leading principal plane and minimizes the total squared reconstruction error over all two-dimensional linear subspaces. To see this, for any orthonormal pair of directions the captured variance is the sum of its two Rayleigh quotients, at most ${\lambda }_1+{\lambda }_2$ by diagonalizing $R$; the remaining variance is the total minus that captured sum. Distances in the plot omit the smaller components, so a close projected pair need not be equally close in every original rating. The actual 43 coordinates cannot be recovered from the summary alone; the supplied data matrix or fitted score matrix is needed.

The printed descriptive list omits INTG, although the loading table includes it. INTG denotes integrity and supplies the eleventh variable. The standard data set also contains a contacts variable; the stated eleven-rating analysis uses the ratings, not contacts.

The response is the species categorical variable, while the four measurements are candidate predictors. The dot in the formula means all other columns of the data frame. The fitted object is a classification tree: it recursively partitions predictor space by binary threshold rules and assigns a species and class-probability vector to each terminal region. Sepal width was available but was not selected in the displayed tree.

At a node $t$, let $n_{tk}$ be the number of observations of class $k$ and $n_t={\sum }_k{n}_{tk}$. The multinomial likelihood is maximized by ${\hat{p}}_{tk}=n_{tk}/n_t$. The fitted class is one attaining the largest proportion, with a software convention resolving ties. The classification-tree deviance isA pure node has zero deviance. At each nonterminal node, consider candidate split points between distinct observed predictor values, and choose the split with the largest decrease $D(t)-D(t_L)-D(t_R)$. Repeating this local maximization constructs a binary decision tree; it does not solve a global optimization over all trees. Growth stops when nodes are pure or size/deviance controls prohibit further splitting. Mixed leaves with five or six observations can therefore remain. The node number $j$ has children $2j$ and $2j+1$, and a terminal-node marker means that no further split was retained.

The probability triples are in class order $(c,s,v)$, not $(s,c,v)$. For example, the pure left leaf has probability vector $(0,1,0)$ and predicts $s$, while the 54-observation node has counts $(49,0,5)$ and predicts $c$. The root counts are $(50,50,50)$, giving $D=300\log 3\simeq 329.584$ and a tied majority class. Its reported choice of $c$ does not indicate a more common class.

The complete derived decision rule is shown below. Each leaf includes its fitted class and the observed class counts, so both the branching logic and its remaining errors are visible. For the observed data, the printed strict inequalities leave no equality cases because the thresholds fall between measurement values. A consistent right-branch convention at equality extends the rule to new specimens.

The six terminal nodes contain the class countsThus the training misclassification error isThe terminal classification-tree deviance is insteadThe displayed residual mean deviance divides this by $150-6=144$, yielding approximately $0.1253$. This is the summary's reporting convention, not a normal-error residual mean square or an automatic chi-squared calibration with 144 statistical degrees of freedom. The two summary measures answer different questions: deviance depends on the fitted probability of each observed class, while misclassification counts only whether the majority-class prediction is correct. Both are training quantities. Cross-validation or independent test data, rather than the displayed training error, are needed to assess prediction on new specimens; cost-complexity tree pruning can select a less elaborate tree if predictive performance warrants it.

A randomized complete block design groups similar experimental units into experimental blocks and places every treatment once in each experimental block, independently randomizing treatment allocation within experimental blocks. Blocking removes background between-block variation from within-block treatment comparisons. A balanced incomplete block design has $v$ treatments and $b$ experimental blocks of size $k<v$, with each treatment occurring in $r$ experimental blocks and every distinct treatment pair together in $\lambda$ experimental blocks. Counting treatment occurrences and co-occurrences gives $vr=bk$ and $r(k-1)=\lambda (v-1)$.

A Latin square of order $q$ places each of $q$ treatment symbols once in every row and every column of a $q\times q$ array. This controls two nuisance directions. A Graeco-Latin square superposes two mutually orthogonal Latin squares: each square is Latin, and each ordered pair of symbols from the two squares occurs exactly once. The two symbol factors can then be fitted together with row and column effects. Random permutations of row labels, column labels and treatment symbols provide a suitable randomized allocation without destroying these balances.

For the operator experiment, fit the additive normal linear modelHere $r_i,c_j$ account for the two fabric directions, and ${\tau }_o$ is the operator effect. Additivity, independent equal-variance errors and the normal assumption supply the exact F-test; there is only one observation per cell, so an unrestricted interaction model is not available.

In the analysis of variance for a Latin square, the three centered factor spaces each have dimension $4-1=3$. Together with the intercept the model has rank ten, leaving $16-10=6$ residual statistical degrees of freedom. The corrected total has $16-1=15$. The completed calculation isUnder the null hypothesis ${\tau }_1=\cdots ={\tau }_4=0$, the operator F-statistic isThe upper five-percent critical value for $F_{3,6}$ is $4.76$, so the operator effect is not significant at five percent in this model.

The decision is independent of the order of fitting rows, columns and operators. To prove the required orthogonality, let $u_i$ be any centered row coefficients and $w_o$ any centered operator coefficients. Their observation-space inner product is ${\sum }_{i,j}{u}_i{w}_{o(i,j)}={\sum }_i{u}_i{\sum }_o{w}_o=0$, since every operator occurs once in each row. The same calculation applies to columns and operators, and all row-column pairs occur once as well. Thus the three centered spaces are mutually orthogonal; their projections commute and sequential sums of squares equal adjusted sums of squares in every fitting order.

With machines included, the layout is a Graeco-Latin square: the machine symbols are Latin, and all sixteen operator-machine pairs occur exactly once. Fitin addition to the earlier constraints. The machine space contributes three new statistical degrees of freedom and is orthogonal to the previous factor spaces, so the operator sum of squares remains $18.25$. The residual sum of squares and statistical degrees of freedom becomeHenceThis exceeds $F_{3,3}(0.05)=9.28$ but not $F_{3,3}(0.01)=29.46$. Therefore the operator effect is now significant at five percent, but not at one percent. Adjustment for machines has removed substantial residual variation without altering the operator numerator. This illustrates why a balanced nuisance factor can change precision and significance even when it does not change the estimated treatment contrasts.

Code each factor level by $x_j\in \{-1,1\}$. A nonconstant factorial contrast indexed by a nonempty factor set $A$ has sign column $g_A(x)={\prod }_{j\in A}{x}_j$. Flipping any coordinate in $A$ reverses that sign and pairs the treatment combinations, so exactly $2^{m-1}$ have each sign. Allocating the two signs to two experimental blocks therefore divides the full factorial design into equal experimental blocks. The contrast is confounded with experimental blocks because its column is constant within each experimental block: changing its coefficient and compensating by the two experimental block intercepts leaves every fitted value unchanged. It has no information separate from experimental block differences. Randomize the allocation of the two groups to days and the run order within days.

Write $A,B,C$ for temperature, pressure and catalyst, respectively, with low levels and catalyst one coded minus. In the supplied allocation, $ABC=-1$ on day one and $ABC=1$ on day two. Thus the three-factor interaction $ABC$ is confounded with days. Every other nonconstant factorial column is orthogonal to the day contrast, because its product with $ABC$ is a nonconstant sign column on the full cube and consequently sums to zero.

If all interactions are negligible, the fitted model contains an intercept, one day contrast and the three main effects. The corrected total statistical degrees of freedom partition asThe three residual statistical degrees of freedom are the unused two-factor interaction columns $AB,AC,BC$, pooled as error under the stated negligible-interaction assumption. There is no independent pure-error replication of treatment combinations.

For four days with two runs each, all three main effects can be estimated. One regular allocation pairs each treatment combination with its complete opposite:The block-generating columns can be $AB$ and $AC$; their product $BC$ is also constant in experimental blocks. Thus $AB,AC,BC$ span the three experimental block contrasts, while $A,B,C$ and $ABC$ vary within experimental blocks. The main effects remain mutually orthogonal and separate from the experimental block effects; with an additive main-effect model there is one residual degree of freedom, represented by $ABC$.

For the further request about $AC$, the distinction between regular confounding and unrestricted pairing is important. A regular four-block allocation cannot retain all three main effects and $AC$. Its three confounded columns form a rank-two defining subgroup. To preserve main effects, these columns must lie in $\{AB,AC,BC,ABC\}$. If the subgroup includes $ABC$ and a two-factor column, their product is a main-effect column, which is forbidden. Hence the only permissible nonidentity subgroup is $\{AB,AC,BC\}$, and it necessarily confounds $AC$. This proves impossibility in the usual regular factorial-blocking framework, rather than relying on a count of available statistical degrees of freedom.

The literal request does not explicitly restrict arrangements to regular experimental blocks. If arbitrary pairings are permitted, a nonregular construction can estimate $A,B,C,AC$ together, assuming the other interactions are omitted from the model. An explicit allocation isEvery treatment combination occurs exactly once. Fit one unrestricted intercept per day and the four coefficients ${\beta }_A,{\beta }_B,{\beta }_C,{\beta }_{AC}$. By estimability from within-block differences, subtracting the second run from the first in each experimental block removes all day effects and yields the coefficient matrixThus all four coefficients are uniquely determined from the four day differences. In this broader design class the answer is yes, under the reduced four-effect model. The resulting fit is saturated: four experimental block intercepts plus four treatment coefficients use all eight observations, leaving no residual estimate of error variance. Replication or external error information would be needed for inference. This does not contradict the impossibility for regular orthogonal confounding, and it does not claim simultaneous unbiased estimation when arbitrary additional interactions are retained.

Let $f(x)=(f_1(x),\dots ,f_p(x){)}^T$ and let an approximate experimental design $\xi$ specify the proportions of observations allocated to points of the design region. Its normalized information matrix of an experimental design isFor $N$ independent equal-variance observations realizing those proportions, $X^{T}X=NM(\xi )$, and the least-squares estimator has covariance matrix ${\sigma }^{2}M(\xi {)}^{-1}/N$. A D-optimal design maximizes $\det M(\xi )$, equivalently minimizing the determinant of this covariance at fixed $N$. A G-optimal design minimizes the largest variance of a fitted regression mean over the design region, equivalentlywhere $d$ is the design sensitivity function. This is the variance of the estimated mean multiplied by $N/{\sigma }^2$, not the variance of a new noisy observation.

For continuous regressors spanning a $p$-dimensional parameter space on a compact design region, the general equivalence theorem for optimal design states that, among approximate designs with nonsingular information matrices, the following are equivalent: being a D-optimal design, being a G-optimal design, andAt every support point of an optimal design, $d=p$. To see the criterion, the sensitivity average is $\int dd\xi =\mathrm{tr}(M^{-1}M)=p$, so the supremum is at least $p$. Adding an infinitesimal allocation at $x$ has directional derivative $d(x,\xi )-p$ for $\log \det M$. At a D-optimum each derivative is nonpositive. Conversely, if all sensitivities are at most $p$, the tangent bound for the concave function $\log \det M$ shows that no competing design increases it. A nonsingular D-optimum exists by the spanning assumption and compactness of the information-matrix set, so the smallest possible G-value is also $p$. Equality at support points follows from the average. The theorem refers to the full approximate-design class; constrained integer allocations or mandated augmentations require separate comparison.

For the linear mean ${\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2$, take $f(x)=(1,x_1,x_2{)}^T$. Equal replication at the four corners givesAll odd corner sums and the sum of $x_1{x}_2$ vanish, while each squared regressor sum is $4m$. If ${\bar{Y}}_{ab}$ is the mean of the $m$ observations at $(a,b)$, $a,b\in \{-1,1\}$, the least-squares normal equations giveIn particular the intercept is the average of the four cell means, and each slope is the average of its sign-weighted cell means. Independence and equal error variance suffice for this covariance formula; normal errors are not needed for it.

The normalized corner design has $M=I_3$, so its sensitivity iswith equality at all four corners. Since $p=3$, the general equivalence theorem for optimal design proves that the equally weighted corner design is D-optimal and G-optimal.

For the augmentation, all $m$ added observations share the regressor $f=(1,x_1,x_2{)}^T$. The enlarged design matrix satisfiesBy the matrix determinant lemma,For the permitted nine points, the determinant is $80m^3$ at the center, $96m^3$ at an edge midpoint and $112m^3$ at a corner. ThereforeAll four corners tie for the best mandated single-point augmentation. This is optimality among the specified augmented designs, not unrestricted D-optimality with $5m$ observations: the resulting normalized determinant is $112/125<1$, whereas equal corner proportions have determinant one. The restriction that the extra $m$ runs all occur at one point is what prevents recovering equal proportions.