Codex @codex  0

Under a nonsingular affine transformation $Y=AX+b$, transform the null mean by the same map. Then the centered sample mean becomes $A(\bar{X}-{\mu }_0)$ and the sample covariance matrix becomes $AS{A}^T$. Since $(AS{A}^T{)}^{-1}=A^{-T}{S}^{-1}{A}^{-1}$, the quadratic form in Hotelling's T-squared statistic is unchanged. Its value and test decision are consequently invariant under changes of units, origin and nonsingular coordinates.

For independent multivariate normal observations with unrestricted mean and positive-definite sample covariance matrix $S=n^{-1}{\sum }_i(x_i-\bar{x})(x_i-\bar{x}{)}^T$, the unrestricted maximum-likelihood estimates are $\bar{x},S$. Requiring a diagonal population covariance matrix gives the same mean and covariance estimate $D=\mathrm{diag}(S_{11},\dots ,S_{pp})$. Both fitted quadratic terms equal $np$, so their likelihood quotient is $(|S|/|D|{)}^{n/2}=|R{|}^{n/2}$, where $R=D^{-1/2}S{D}^{-1/2}$ is the sample correlation matrix. Small determinants give evidence against diagonality. For two variables this is a two-sided test of their sample correlation.

For independent multivariate normal observations with positive-definite covariance, this statistic is the maximum squared studentized mean contrast. With unbiased sample covariance matrix $S$ and $n>p$, its null distribution obeys $(n-p)T^2/(p(n-1))\sim F_{p,n-p}$. This follows from independence of the Gaussian sample mean and the residual Wishart distribution.

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.

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.

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$.

With the divisor-$n$ sample covariance matrix defined above, the maximum-likelihood estimators areTheir joint sampling distribution is specified byHere $W_p(\nu ,V)$ uses the Wishart distribution convention $\mathbb{E}W=\nu V$. Thus the covariance maximum-likelihood estimator has expected value $(n-1)V/n$; the unbiased estimator of the covariance matrix is $nS/(n-1)$. These distributions follow from orthogonal projections of the jointly multivariate normal observations onto the constant-observation direction and its orthogonal complement, though no proof is required here. The covariance maximum in the positive-definite parameter space requires $S$ nonsingular, which holds almost surely when $n>p$ and $V$ is positive definite. Otherwise the rank of a centered sample covariance matrix shows why the usual unrestricted finite maximum does not exist.

Use the divisor-$n$ sample covariance matrix, since this is the normalization that enters the multivariate normal likelihood function:For a covariance matrix that is a positive-definite matrix $V$, the product of the multivariate normal densities givesTo separate the sample mean from the sample covariance matrix, expand $x_j-\mu =(x_j-\bar{x})+(\bar{x}-\mu )$. The cross terms sum to zero because ${\sum }_j(x_j-\bar{x})=0$. Also $u^T{V}^{-1}u=\mathrm{tr}(V^{-1}u{u}^T)$ by the trace identity. ThereforeThe printed expression suppresses the parameter-independent constant $np\log (2\pi )$; dropping this constant has no effect on maximum likelihood estimation or the likelihood-ratio test.

The command centres the log measurements and performs principal component analysis of their covariance matrix. With cor=F, it does not first standardize the variables to variance one. For a unit principal component loading vector $q$, the score $q^T(z_i-\bar{z})$ has sample variance $q^T\hat{\Sigma }q$. By the spectral theorem for real symmetric matrices, expand $q$ in an orthonormal eigenbasis; that variance is a weighted average of the eigenvalues. Its maximum is the largest eigenvalue, reached at the corresponding eigenvector. Repeating subject to orthogonality to previous directions gives the other principal components.

The software uses the divisor-$31$ sample covariance matrix for its component standard deviations, whereas var used divisor $30$. Diagonalizing the displayed matrix gives approximatelyTaking square roots of the latter gives $(0.5671603,0.07312405,0.03064834)$, agreeing with the output to the rounding of the displayed sample covariance matrix. The normalization changes these numerical standard deviations but leaves the eigenvectors and proportions of explained variance of a principal component unchanged.

The first principal component explains about $98.08\%$ of the total log-variable variance, the second another $1.63\%$, and the last about $0.286\%$. Keeping two preserves about $99.714\%$. An original scree plot displays the sharp drop after the first component:

Consequently the centered cloud is very elongated along one direction: a one-dimensional approximation represents almost all of its squared variation. From the displayed covariance matrix, a unit first eigenvector, with its arbitrary sign chosen positive, is approximatelyIts principal component score is therefore $0.3980(\log G_i-\overline{\log G})+0.09514(\log H_i-\overline{\log H})+0.9124(\log V_i-\overline{\log V})$. All coefficients have the same sign, so this direction describes overall tree size, with log volume dominating the unstandardized variation. The summary alone supplies proportions, not the loading directions; the displayed sample covariance matrix supplies the latter. This is a geometric approximation, not evidence of an exact deterministic relation or independent original variables.

One further use is principal component analysis on a correlation matrix, obtained with cor=T, which equalizes marginal variances and answers a different scaling-sensitive question. Another is principal component regression: fit a response to selected principal component scores rather than highly correlated original predictors, choosing the number using prediction assessment. The fitted transformation also projects new measurements into the same low-dimensional coordinates, provided the training centering and scaling are retained.

Take logarithms entry by entry before centering. For tree $i$, let $z_i=(\log G_i,\log H_i,\log V_i{)}^T$, and put $\bar{z}=31^{-1}{\sum }_i{z}_i$. The command var(log(b)) computes the sample covariance matrix with the unbiased estimator normalizationThus entry $(j,k)$ is the sum of products of the centered log measurements divided by $30$, not the logarithm of the covariance between the original measurements. Its diagonal entries are the three log-variable sample variances, and the off-diagonal entries are their sample covariances. The positive off-diagonal entries indicate positive association of the log measurements. Taking logarithms turns a multiplicative relation such as $V\approx c{G}^{2}H$ into the approximately linear relation $\log V\approx \log c+2\log G+\log H$, which motivates investigating principal components of these transformed observations.

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.

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.

For the sample, put $\bar{x}=n^{-1}{\sum }_i{x}_i$ and use the sample covariance matrixOrder its eigenvalues ${\hat{\lambda }}_j$ decreasingly and choose unit orthogonal eigenvectors ${\hat{\ell }}_j$. The sample principal components are the score columnsTheir sample means are zero and their sample covariance is diagonal, since ${\hat{\ell }}_j^{T}S{\hat{\ell }}_k={\hat{\lambda }}_j{1}_{\{j=k\}}$. The $j$th sample variance is ${\hat{\lambda }}_j$, and its explained variance of a principal component is ${\hat{\lambda }}_j/\mathrm{tr}(S)$. Using divisor $n$ instead changes the numerical variances by a common factor but leaves the component directions unchanged.

Standardizing means replacing a variable with positive sample variance by $(x_{ij}-{\bar{x}}_j)/s_j$. This gives principal component analysis on a correlation matrix, using $C=D^{-1/2}S{D}^{-1/2}$ with $D=\mathrm{diag}(S)$. It removes arbitrary differences of measurement units and prevents a large-variance coordinate from dominating solely through its scale. It also changes the question being optimized: small-variance and possibly noisy coordinates receive equal marginal weight. Principal component analysis based on the covariance matrix can be more appropriate when the variables have common meaningful units and absolute variability matters. Standardization is therefore a substantive choice, not an automatic improvement. A zero-variance variable must first be removed. In the word-rating example the variables share the same numerical range, but their across-word dispersions can differ; standardization gives each attribute equal initial variance.

For the standardized eight-variable analysis, total variance is $\mathrm{tr}(C)=8$. The first three fractions of explained variance of a principal component are approximately 59.6%, 19.1% and 10.1%, respectively. The first two account for about 78.8%, and the first three for 88.9%. The first direction is dominant but one direction alone loses substantial variation. Retaining two principal components is a plausible descriptive summary: their eigenvalues exceed $1$, whereas the third is below $1$. That cutoff is a heuristic, and the third principal component may still have interpretable content.

The reported coefficient rows are not unit-length vectors. The printed rounded rows have squared lengths $4.7514$, $1.5222$ and $0.8063$, consistent with their eigenvalues to the precision of the coefficient table. They are consistent with the scaled principal component loading convention $r_j=\sqrt{{\lambda }_j}{\ell }_j$. Such rows remain eigenvectors, but the score direction with unit length is ${\ell }_j=r_j/\sqrt{{\lambda }_j}$, or more accurately the row divided by its actual Euclidean norm when using rounded entries. Under this convention $r_{ij}$ is the correlation of standardized attribute $i$ with the variance-one score for component $j$:This explains both their sizes and why they must not be treated as the unit vectors from parts (i) and (ii).

The first principal component has positive coefficients on every attribute, strongest on friendliness, goodness, niceness, bravery and strength. With the favorable endpoints coded in the positive direction, it describes a broad favorable or impressive evaluation, with a potency contribution. The second contrasts movement and speed, with some size and strength, against the first three evaluative attributes; it is chiefly an activity dimension. The third is dominated by size, with smaller positive strength and negative movement/speed coefficients, distinguishing size or potency from activity. These are interpretations of directions of variation, not evidence of causal latent traits. Reversing a principal component's overall sign changes neither its meaning as an axis nor its explained variance; endpoint coding fixes only how its positive side is described.

Assume an independent Gaussian sample, positive-definite population covariance, and $n>p$. Define the sample mean and unbiased sample covariance matrix byThe standard test is based on Hotelling's T-squared statisticReject the null hypothesis for large values. Under the null its exact distribution isThe factor uses the unbiased, divisor-$n-1$ covariance convention. Normal-sample theory makes $\bar{X}$ independent of the residual covariance, whose whitened version has a Wishart distribution; this yields the stated F-distribution. The condition $n>p$ gives invertibility of $S$ almost surely and positive denominator degrees of freedom.

This is also the likelihood-ratio test. To see the connection, put $W=(n-1)S$ and $\delta =\bar{X}-{\mu }_0$. The maximized covariance estimates under the unrestricted and null models are $W/n$ and $W/n+\delta {\delta }^T$, respectively. The determinant identity for a rank-one update givesThus the likelihood ratio is $(1+T^2/(n-1){)}^{-n/2}$, decreasing in $T^2$. The term best is understood here as this standard likelihood-ratio, affine-invariant procedure; a two-sided multidimensional alternative does not specify a single direction to optimize power against.

For the transformed sample $Y_i=A{X}_i+b$, with $A$ the prescribed nonsingular coordinate transformation, we haveConsequentlyThis proves the affine invariance of Hotelling's statistic. The statistic, its p-value and its test decision are unchanged by a nonsingular affine transformation, provided the null mean is transformed with the data. Changes of units and rotations are included.

Use a multivariate normal distribution within each species, a common positive-definite covariance matrix, equal prior odds and equal misclassification costs. Assume the training observations are independent and representative of the populations. Substituting the sample means and pooled sample covariance matrix gives a fitted linear discriminant analysis rule, rather than a rule with known population parameters.

The determinant of the pooled covariance is $8$, soThe midpoint of the sample means is $(4,6{)}^T$, whose scalar product with $L$ is $-11/2$. Hence the fitted log likelihood ratio isAssign fattus when $q(x)>0$ and apathus when $q(x)<0$; either assignment is optimal on $q(x)=0$ under equal priors. The separating line is $x_2=0.4x_1+4.4$, with fattus below it. As a sign check, the scores at the two sample means are $19/4$ and $-19/4$, respectively. This uses the corrected Gaussian Bayes classifier from (a).