12. Multivariate Statistical Analysis

Introduction and Context

Data science rarely studies variables in isolation. A customer’s purchase value may be related to age, income, browsing behaviour, and payment history; a patient’s diagnosis may depend on several laboratory measurements; and a machine’s failure risk may reflect temperature, vibration, pressure, and operating time. When several variables are observed for each unit, the data are multivariate.

Multivariate statistical analysis provides a systematic framework for understanding such data. It examines means, variances, covariances, correlations, dependence structures, latent dimensions, group membership, and similarity among observations. Its central challenge is that relationships among variables can be informative but also difficult to visualize and interpret, particularly when the number of variables is large.

This chapter develops the principal methods used in applied multivariate analysis: covariance and correlation structures, principal component analysis, factor analysis, clustering, discriminant analysis, and related dimensionality-reduction strategies. It also considers scaling, validation, interpretability, computational tools, and ethical use. These methods extend the ideas of univariate and bivariate statistics into settings where the joint distribution of variables matters.

Learning Objectives

After completing this chapter, the reader should be able to:

1. Distinguish covariance, correlation, dependence, and distance in multivariate data.

2. Assess how scaling, outliers, missing values, and multicollinearity affect analysis.

3. Explain and implement the logic of principal component analysis.

4. Interpret factor models and distinguish common factors from principal components.

5. Select and evaluate clustering methods for exploratory segmentation.

6. Apply discriminant analysis to supervised classification and understand its assumptions.

7. Compare dimensionality-reduction methods according to analytical purpose.

8. Communicate multivariate results without overstating patterns or implying unsupported causality.

The Geometry of Multivariate Data

Suppose that n observations are measured on p variables. The data matrix X has n rows and p columns. A row represents an observation and a column represents a variable. The sample mean vector is

x̄ = (x̄₁, x̄₂, …, x̄ₚ)′.

Unlike a single mean, a mean vector does not describe how variables vary together. That information is contained in the sample covariance matrix:

S = 1/(n − 1) (X − 1x̄′)′(X − 1x̄′).

The diagonal entries of S are sample variances, while the off-diagonal entries are covariances. A positive covariance indicates that two variables tend to be above or below their means together; a negative covariance indicates an opposing tendency. Covariance, however, depends on measurement units. Changing income from rupees to thousands of rupees changes its covariance with another variable even though the substantive relationship is unchanged.

The correlation matrix standardizes covariance:

Rᵢⱼ = Sᵢⱼ/(sᵢsⱼ).

Correlations lie between −1 and 1 and are therefore easier to compare across variables. Nevertheless, a correlation measures linear association, not causation, and a weak correlation can conceal nonlinear or subgroup-specific relationships.

Multivariate data can also be viewed geometrically. Each observation is a point in p-dimensional space, and each variable is a coordinate axis. Similarity can be represented by a distance. Euclidean distance is commonly defined as

d(x, y) = [Σⱼ(xⱼ − yⱼ)²]¹ᐟ².

Manhattan distance sums absolute coordinate differences. Mahalanobis distance additionally accounts for covariance:

D²(x) = (x − μ)′Σ⁻¹(x − μ).

Thus, a deviation in a highly variable direction may be less unusual than the same numerical deviation in a stable direction. Mahalanobis distance is useful for multivariate outlier detection, provided that the covariance matrix is estimated reliably.

Scaling and Preprocessing

Scaling is an analytical decision, not merely a software setting. If variables are measured in incomparable units, standardization is often appropriate:

zⱼ = (xⱼ − x̄ⱼ)/sⱼ.

Using the correlation matrix is equivalent to analysing standardized variables. However, standardization gives a low-variance variable the same potential influence as a high-variance variable. This may be desirable when variables are conceptually equally important, but inappropriate when measurement variance itself carries meaning.

Robust scaling, based on the median and interquartile range, can reduce the influence of extreme observations. Logarithmic or other transformations may be useful for highly skewed variables. Missing values require explicit treatment because deletion, imputation, or model-based approaches can produce different covariance structures. Imputation should be performed within the training data when predictive analysis is involved, so that information from the test set does not leak into model construction.

Principal Component Analysis: Rotating the Coordinate System

Principal component analysis (PCA) is an unsupervised method that replaces correlated variables with a smaller number of uncorrelated linear combinations called principal components. The first component is the direction in variable space along which the projected observations have maximum variance. Each subsequent component maximizes remaining variance subject to being orthogonal to earlier components.

If S is the covariance matrix, PCA solves the eigenvalue problem

Svₖ = λₖvₖ,

where vₖ is the kth eigenvector and λₖ is its eigenvalue. The kth component score for observation i is the projection

zᵢₖ = (xᵢ − x̄)′vₖ.

The eigenvalues measure the variance captured by the components. The proportion of total variance explained by component k is λₖ/Σλⱼ. If PCA is applied to a correlation matrix, the total variance equals p; if applied to a covariance matrix, variables with larger units or variances contribute more strongly.

PCA can be understood as a lossy compression method. Retaining the first q components produces an approximation to the original data. The discarded variance is a measure of reconstruction loss under squared-error criteria. PCA is therefore useful for visualization, noise reduction, compression, and constructing predictors with reduced collinearity.

Choosing the Number of Components

No single rule is universally correct. Common tools include a scree plot, cumulative explained variance, parallel analysis, and cross-validated predictive performance. The Kaiser rule, which retains eigenvalues greater than one, is easy to apply to a correlation matrix but can over- or under-retain components. A component should also be justified substantively: a mathematically strong component may have no meaningful interpretation.

Loadings describe the relationship between original variables and components. Large positive loadings indicate variables that increase with the component; large negative loadings indicate an opposing contribution. Sign is arbitrary: multiplying an eigenvector by −1 produces the same component. Analysts should therefore interpret relative patterns, not the sign alone.

PCA is not a causal model. It does not identify underlying mechanisms, and its components may be difficult to explain. It also assumes that linear combinations are useful and can be sensitive to outliers. Sparse PCA and robust PCA address some of these concerns by encouraging simpler loadings or reducing outlier influence.

Factor Analysis: Modelling Latent Constructs

Factor analysis has a different purpose from PCA. It seeks to explain observed correlations through a smaller number of unobserved latent factors. A common factor model is

X = μ + ΛF + ε,

where F is a vector of common factors, Λ is the loading matrix, and ε represents variable-specific error or uniqueness. Under the usual model, common factors explain covariance shared among variables, while uniqueness captures variance specific to an observed variable.

PCA decomposes total variance; factor analysis models common variance. Consequently, the number and interpretation of extracted dimensions may differ. Exploratory factor analysis (EFA) is used when the structure is uncertain. Confirmatory factor analysis (CFA), generally estimated within a structural equation modelling framework, tests a specified measurement model.

Before factor analysis, analysts examine whether correlations are adequate. The Kaiser–Meyer–Olkin statistic assesses sampling adequacy, while Bartlett’s test examines whether the correlation matrix differs from an identity matrix. These are aids rather than substitutes for substantive reasoning and adequate sample size.

Factor extraction may use principal-axis factoring or maximum likelihood. Rotation improves interpretability by seeking a simpler loading pattern. Orthogonal rotations, such as varimax, keep factors uncorrelated; oblique rotations, such as oblimin or promax, permit factor correlations and are often more realistic in social and behavioural applications. A factor loading is not automatically a validated scale: reliability, construct validity, cross-loadings, and replication must be examined.

Clustering: Discovering Groups Without Labels

Clustering methods assign observations to groups so that observations within a cluster are more similar than observations in different clusters. Because clusters are not known in advance, clustering is exploratory rather than inherently confirmatory. The result depends on the variables, scaling, distance measure, algorithm, and chosen number of clusters.

Centroid-Based and Hierarchical Methods

K-means clustering partitions observations into k groups by minimizing within-cluster squared Euclidean distance. The algorithm alternates between assigning observations to the nearest centroid and updating centroids. It is computationally efficient but requires k in advance, is sensitive to initialization and outliers, and works best for relatively compact, similarly shaped clusters.

Hierarchical agglomerative clustering begins with each observation as a separate cluster and successively merges clusters. The sequence is represented by a dendrogram. Linkage rules determine the distance between clusters: single linkage uses the nearest pair, complete linkage the farthest pair, average linkage the average pairwise distance, and Ward’s method seeks to minimize increases in within-cluster variance. Hierarchical methods provide a range of possible solutions but can be computationally expensive and irreversible once a merge occurs.

Density-based methods such as DBSCAN identify dense regions separated by areas of low density. They can find irregularly shaped clusters and label isolated observations as noise, but their results depend on neighbourhood parameters and may be difficult when cluster densities differ. Model-based methods assume that observations arise from a mixture of probability distributions and can provide probabilistic cluster membership, although distributional assumptions must be checked.

Evaluating Cluster Solutions

The silhouette coefficient compares within-cluster cohesion with separation from the nearest alternative cluster. Calinski–Harabasz and Davies–Bouldin indices provide additional internal criteria. These measures should not be treated as automatic evidence that a cluster is real. Stability under resampling, temporal validation, external variables, and domain interpretation are essential. A cluster that is mathematically distinct but operationally meaningless may have little value.

Discriminant Analysis: Classification with Group Structure

Discriminant analysis is a supervised multivariate method. Unlike clustering, it begins with observations whose group labels are known and develops rules for assigning new observations to groups.

Linear discriminant analysis (LDA) assumes, in its classical form, that each class follows a multivariate normal distribution and that all classes share a common covariance matrix. The discriminant function for class k is based on the class mean, common covariance matrix, and prior probability. The method seeks directions that maximize between-class variation relative to within-class variation. With class-specific covariance matrices, quadratic discriminant analysis (QDA) produces curved decision boundaries but requires more parameters and therefore more data.

The assumptions of normality and equal covariance are often approximate in practice. Their violation does not automatically invalidate classification, but performance should be assessed empirically. Regularized discriminant analysis shrinks covariance estimates toward simpler structures, which is particularly useful when variables are numerous or strongly correlated. In high-dimensional settings where p approaches or exceeds n, ordinary covariance estimation may be singular; shrinkage, feature selection, or dimension reduction becomes necessary.

Classification performance should be evaluated with a test set or cross-validation. Accuracy can be misleading under class imbalance, so analysts should examine sensitivity, specificity, precision, recall, balanced accuracy, and the confusion matrix. Posterior probabilities should be calibrated when decisions depend on risk thresholds.

Choosing and Combining Methods

The appropriate method depends on the analytical question. PCA is suitable when the goal is compact representation or visualization. Factor analysis is more appropriate when latent constructs and measurement error are central. Clustering is used to explore unknown groups, while discriminant analysis predicts known groups. These methods can be combined: PCA may reduce correlated predictors before classification, or factor scores may be used in a subsequent regression model.

Other dimensionality-reduction methods serve different purposes. Multidimensional scaling represents observations in a low-dimensional space while preserving pairwise distances as far as possible. Correspondence analysis is designed for categorical contingency tables. Kernel PCA captures certain nonlinear structures through implicit feature mappings. t-distributed stochastic neighbour embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) are useful for visual exploration of complex data, but their two-dimensional maps should not be treated as faithful global representations or as proof of natural clusters. Parameters can substantially affect apparent structure.

For predictive tasks, supervised methods such as partial least squares, regularized regression, tree-based models, or neural representation learning may outperform unsupervised components. Dimensionality reduction should therefore be selected according to the loss function and decision problem, not merely because a plot looks attractive.

An Illustrative Case: Segmenting and Classifying Retail Customers

Background

A fictional online retailer wants to understand customer behaviour and improve retention. For each customer, the analytics team has purchase frequency, average order value, days since last purchase, product-category diversity, discount usage, and customer-service contacts.

Situation and Problem

The variables have different units and are correlated. High purchase frequency is often associated with lower recency, while discount usage may be related to order value. Management wants both an exploratory segmentation and a model that identifies whether a newly observed customer belongs to a strategically important group.

Analysis and Approach

The team first inspects missingness, extreme values, and distributions. Monetary variables are log-transformed because of strong right skew, and robust standardization is used so that a few unusually large customers do not dominate the analysis. The correlation matrix shows substantial redundancy among frequency, recency, and order value.

PCA is applied for visualization and diagnostic purposes. The first component broadly represents purchasing intensity, while a second component contrasts discount dependence with product diversity. The team does not automatically call these components customer “needs”; they are descriptive statistical dimensions.

For segmentation, the team compares k-means solutions with three through six clusters and examines silhouette values, bootstrap stability, and managerial usefulness. A four-cluster solution is selected because it is reasonably stable and yields interpretable groups: high-value loyal customers, discount-sensitive frequent customers, occasional broad-category buyers, and inactive customers. The labels are descriptive rather than inherent truths.

Historical group labels from a validated campaign are then used to train an LDA classifier. Because covariance estimates are unstable with correlated inputs, the team uses standardized variables and compares LDA with a regularized alternative. Cross-validation indicates that the regularized model performs better for the smaller high-value class. Evaluation uses the confusion matrix, class-specific recall, and calibration rather than overall accuracy alone.

Outcome

The retailer develops different retention strategies for the four groups. The model is deployed only as a decision-support tool, and customer-facing interventions are monitored for differential error rates. Segment membership is periodically re-estimated because purchasing behaviour and market conditions change.

Learning Points

The case illustrates that preprocessing affects both clustering and classification; PCA can clarify structure without automatically defining business segments; clusters require stability and substantive validation; and a predictive classifier must be evaluated with metrics appropriate to class imbalance and business consequences.

Challenges, Risks, and Common Mistakes

A common mistake is to apply PCA, clustering, or factor analysis without examining scale and outliers. Another is to interpret correlation as causation or to treat a component or cluster as a naturally existing entity. Analysts may also overfit by selecting the number of dimensions or clusters after repeatedly inspecting the same evaluation data.

High dimensionality creates statistical difficulties. Covariance matrices become unstable, distances become less discriminating, and visualizations can conceal rather than reveal structure. The curse of dimensionality means that additional variables do not always add useful information. Regularization, feature screening based on domain knowledge, larger samples, and validation on new data are practical responses.

Cluster labels can be unstable, particularly near boundaries. Factor solutions may rotate into multiple plausible interpretations. Classification performance may degrade when the operational population differs from the training population. Data leakage is a serious risk when scaling, imputation, feature selection, or dimensionality reduction is performed before train-test separation.

Ethical and Professional Considerations

Multivariate methods can reproduce historical discrimination when sensitive attributes or their proxies are included. Clustering may create labels that influence access to credit, healthcare, employment, or services without a clear justification. Individuals may be classified into groups without knowing how the classification works or how to challenge it.

Responsible practice requires documenting variables, preprocessing decisions, model assumptions, validation results, and intended uses. Analysts should assess performance across relevant demographic groups, protect confidential data, minimize unnecessary attributes, and distinguish exploratory findings from decisions with legal or material consequences. Interpretability is not achieved merely by producing a two-dimensional plot; stakeholders must understand what the dimensions represent, how stable they are, and what uncertainty remains.

Technology and Emerging Developments

Modern statistical software provides covariance diagnostics, eigenvalue decomposition, factor extraction, clustering algorithms, and cross-validation workflows. Reproducible pipelines should record package versions, random seeds, preprocessing parameters, and model configurations. In production systems, monitoring should detect changes in distributions, cluster composition, covariance structure, and classification error.

Emerging work includes robust and sparse covariance estimation, probabilistic latent representations, nonlinear manifold learning, graph-based dimensionality reduction, and representation learning with neural networks. These methods can capture complex structure but often reduce transparency. For high-stakes applications, a slightly less flexible method with clear interpretation and reliable validation may be preferable to a more opaque representation.

Best Practices for Multivariate Analysis

A sound workflow normally includes the following principles:

• Define the analytical purpose before selecting a method.
• Inspect distributions, missingness, outliers, units, and correlations.
• Choose covariance- or correlation-based analysis deliberately.
• Fit preprocessing steps using training data only when prediction is involved.
• Compare solutions using statistical evidence, stability, and domain meaning.
• Validate components, clusters, and classifiers on new or resampled data.
• Report uncertainty, assumptions, sensitivity analyses, and limitations.
• Use visualizations as exploratory evidence rather than definitive proof.
• Document decisions so that results can be reproduced and audited.

Chapter Summary

Multivariate statistical analysis studies several variables jointly and uses covariance, correlation, distance, and geometric representations to reveal structure. PCA transforms correlated variables into orthogonal components that maximize explained variance, whereas factor analysis models shared variation through latent constructs. Clustering explores unknown groups, while discriminant analysis classifies observations with known group labels. Scaling, outliers, missing data, covariance estimation, and dimensionality all strongly affect results. Method selection must follow the analytical objective, and apparent structure should be tested for stability, interpretability, and predictive or operational value. Ethical practice requires attention to fairness, transparency, privacy, uncertainty, and the consequences of classification.

Key Terms

Covariance: A measure of how two variables vary jointly.

Correlation matrix: A matrix of standardized pairwise associations among variables.

Mahalanobis distance: A covariance-adjusted distance measuring multivariate deviation.

Standardization: Transformation of a variable to have mean zero and standard deviation one.

Principal component: An orthogonal linear combination of variables ordered by explained variance.

Eigenvalue: The variance associated with a principal component or eigenvector direction.

Loading: The contribution or association of an observed variable with a component or factor.

Explained variance: The proportion of total variability represented by a component or reduced model.

Factor analysis: A model-based method explaining observed covariance through latent factors.

Uniqueness: Variance in an observed variable not explained by common factors.

Rotation: A transformation of factor axes intended to improve interpretability.

Clustering: Unsupervised assignment of observations to groups based on similarity.

Centroid: The mean position of observations assigned to a cluster.

Linkage: A rule defining distance between clusters in hierarchical clustering.

Silhouette coefficient: An index comparing within-cluster cohesion with separation.

Discriminant analysis: Supervised classification based on multivariate group distributions.

Linear discriminant analysis: Discriminant analysis using a common covariance matrix across classes.

Quadratic discriminant analysis: Discriminant analysis allowing class-specific covariance matrices.

Dimensionality reduction: Representation of data using fewer variables or dimensions.

Data leakage: Use of information during model development that would not be available at prediction time.

Review Questions

Short Answer Questions

1. Why is covariance dependent on measurement units, and how does correlation address this issue?

2. What does the first principal component maximize?

3. Distinguish total variance in PCA from common variance in factor analysis.

4. Why is scaling important before k-means clustering?

5. State one key difference between LDA and QDA.

Long Answer Questions

1. Explain the eigenvalue and eigenvector formulation of PCA and discuss methods for selecting the number of components.

2. Compare exploratory factor analysis with PCA in terms of purpose, variance decomposition, and interpretation.

3. Describe k-means, hierarchical clustering, and density-based clustering, including situations in which each may be appropriate.

4. Explain the assumptions of discriminant analysis and the consequences of covariance instability in high-dimensional data.

5. Develop a complete workflow for a multivariate analysis, from data preparation through validation and communication of results.

Analytical / Application Questions

1. A dataset contains income, age, transaction count, and account balance, measured in very different units and containing extreme incomes. Recommend a preprocessing and PCA strategy, and justify whether covariance or correlation should be used.

2. A clustering analysis produces six mathematically distinct customer groups, but the solution changes substantially under bootstrap resampling. How should the analyst respond before using the groups operationally?

3. A classifier achieves 96% accuracy when only 4% of cases belong to the important minority class. Explain why this result may be misleading and propose a more appropriate evaluation plan.

22. Multiple Choice Questions with Answers

Question 1

Which statement best describes the first principal component?

A. It is the variable with the largest original variance.

B. It is the linear combination that maximizes projected variance subject to a unit-length coefficient vector.

C. It is the factor that explains only measurement error.

D. It is the cluster centroid farthest from the global mean.

Correct Answer: B. It is the linear combination that maximizes projected variance subject to a unit-length coefficient vector.

Question 2

A researcher changes income from rupees to thousands of rupees before computing covariance. What changes?

A. The covariance changes in numerical value, while the correlation remains unchanged.

B. Both covariance and correlation necessarily become zero.

C. The correlation changes but covariance remains unchanged.

D. Neither covariance nor correlation changes.

Correct Answer: A. The covariance changes in numerical value, while the correlation remains unchanged.

Question 3

Why might an analyst choose factor analysis rather than PCA?

A. To guarantee that all original variables are independent.

B. To model shared covariance through latent constructs and separate it from variable-specific variance.

C. To ensure that clusters are spherical.

D. To eliminate the need for validation.

Correct Answer: B. To model shared covariance through latent constructs and separate it from variable-specific variance.

Question 4

Which clustering method can identify irregularly shaped dense groups and label isolated observations as noise?

A. K-means

B. Ward hierarchical clustering

C. DBSCAN

D. Linear discriminant analysis

Correct Answer: C. DBSCAN

Question 5

A dataset contains variables measured in dollars, counts, and days. If the analyst wants each variable to have comparable influence on Euclidean distance, the most appropriate first step is usually:

A. Standardization

B. Removing the variable with the smallest mean

C. Replacing all values with ranks without inspection

D. Applying QDA

Correct Answer: A. Standardization

Question 6

Which assumption distinguishes classical LDA from QDA?

A. LDA requires no group labels.

B. LDA assumes a common covariance matrix across classes, whereas QDA permits class-specific covariance matrices.

C. LDA is unsupervised, whereas QDA is supervised.

D. LDA requires all variables to be categorical.

Correct Answer: B. LDA assumes a common covariance matrix across classes, whereas QDA permits class-specific covariance matrices.

Question 7

A four-cluster solution has a high silhouette score but is unstable under resampling and cannot be explained to domain experts. What is the best conclusion?

A. The clusters are definitively valid because the silhouette score is high.

B. The analyst should report the clusters as causal types.

C. The solution requires further investigation and should not yet guide high-stakes decisions.

D. The number of clusters must automatically be increased.

Correct Answer: C. The solution requires further investigation and should not yet guide high-stakes decisions.

Question 8

Why can applying imputation and standardization to the full dataset before cross-validation create a problem?

A. It guarantees underfitting.

B. It can leak information from validation folds into the training process.

C. It makes covariance matrices diagonal.

D. It prevents any model from being trained.

Correct Answer: B. It can leak information from validation folds into the training process.

Question 9

Which statement about component signs in PCA is correct?

A. A component with negative loadings is always invalid.

B. The sign is arbitrary because multiplying an eigenvector by negative one gives the same direction.

C. Signs determine causality.

D. Signs are fixed by the measurement units.

Correct Answer: B. The sign is arbitrary because multiplying an eigenvector by negative one gives the same direction.

Question 10

Why can overall accuracy be misleading for a classifier when the minority class is rare?

A. Accuracy cannot be computed for supervised methods.

B. A model can obtain high accuracy by mostly predicting the majority class while missing minority cases.

C. Accuracy always penalizes majority-class predictions.

D. Rare classes automatically make covariance singular.

Correct Answer: B. A model can obtain high accuracy by mostly predicting the majority class while missing minority cases.