Notes for Applied Multivariate Statistical Analysis

lectured by Fang Yao∗ LATEXed by T.-Y. Li†

2019 Fall

Abstract
The main goal of this course is to study multivariate statistical techniques, their conventional and
modern ideas, innovative statistical methods, and novel computational tools, as well as new applica-
tions. The topics covered include Multivariate Modeling and Inferences, Model Assessment, Principal
Components, Dimension Reduction, Classification, Cluster Analysis, if time permits, Nonparametric
Methods, Support Vector Machines, Latent Variable Models, Artificial Neural Networks and so on.

All rights reserved. Please feel free to leave a message at [Link]

(under construction). Any comments are welcome. f

∗ Homepage: [Link] E-Mail: fyao@[Link]

† E-Mail: kellty@[Link]

i
 Notes for Applied Multivariate Statistical Analysis

Contents
§1 Review of Linear Algebra (2019/9/10) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 1
§2 Review of Linear Algebra Cont’d (2019/9/17) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 1
§3 Prediction (2019/9/19) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 3
§4 Prediction Cont’d & Multivariate Data (2019/9/24) · · · · · · · · · · · · · · · · · · · · · · · · · · 3
§5 Descriptive Statistics & Normal Distribution (2019/10/8) · · · · · · · · · · · · · · · · · · · · · · 4
§6 Multivariate Normal Distribution Cont’d & MLE (2019/10/15) · · · · · · · · · · · · · · · · · 6
§7 Maximum Likelihood Estimation Cont’d (2019/10/17) · · · · · · · · · · · · · · · · · · · · · · · · · 6
§8 Diagnostics for Normality & Matrix Distribution (2019/10/22) · · · · · · · · · · · · · · · · · 7
§9 Cochran’s Theorem & Hypothesis Testing (2019/10/29) · · · · · · · · · · · · · · · · · · · · · · · 9
§10 Hypothesis Testing Cont’d (2019/10/31) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 10
§11 Inferences about Multi-Population Means (2019/11/5) · · · · · · · · · · · · · · · · · · · · · · · · 11
§12 Multivariate Analysis of Variance (2019/11/12) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 13
§13 Linear Regression (2019/11/14) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 16
§14 Multivariate Linear Regression Cont’d (2019/11/19) · · · · · · · · · · · · · · · · · · · · · · · · · · 20
§15 Model Diagnostics and Selection (2019/11/26) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 22
§16 Model Selection Cont’d & Principal Components (2019/11/28) · · · · · · · · · · · · · · · · · 24
§17 Principal Component Analysis Cont’d (2019/12/3) · · · · · · · · · · · · · · · · · · · · · · · · · · · 26
§18 Canonical Correlation & Factor Analysis (2019/12/10) · · · · · · · · · · · · · · · · · · · · · · · · 27
§19 Factor Analysis Cont’d & Classification (2019/12/12) · · · · · · · · · · · · · · · · · · · · · · · · · 29
§20 Classification Cont’d (2019/12/17) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 30
§21 Classification Cont’d (2019/12/24) · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 32
§22 Classification Cont’d & Final Review (2019/12/26) · · · · · · · · · · · · · · · · · · · · · · · · · · · 34

References
[JW] Richard A. Johnson, & Dean W. Wichern (2007). Applied Multivariate Statistical Analysis (6th ed.).
Pearson Prentice Hall. [Link]
[I] Alan J. Izenman (2008). Modern Multivariate Statistical Techniques: Regression, Classification, and
Manifold Learning. Springer-Verlag. [Link]
[JWHT] Gareth James, Daniela Witten, Trevor Hastie, & Robert Tibshirani (2013). An Introduction to
Statistical Learning, with Applications in R. Springer Science+Business Media. [Link]
1007/978-1-4614-7138-7
[G] Hui-Xuan Gao (2005). Applied Multivariate Statistical Analysis (in Chinese). Peking University Press.
[HS] Wolfgang K. Härdle, & Léopold Simar (2019). Applied Multivariate Statistical Analysis (5th ed.). Springer
Nature. [Link]
[K] Inge Koch (2013). Analysis of Multivariate and High-Dimensional Data. Cambridge University Press. https:
//[Link]/10.1017/CBO9781139025805
[V] Matthias Vallentin (Dec 19th , 2017). Probability and Statistics Cookbook. [Link]
[PP] Kaare B. Petersen, & Michael S. Pedersen (Nov 15th , 2012). The Matrix Cookbook. [Link]
[Link]/viewdoc/summary?doi=[Link].6244

[Z] Topics in Matrix Theory (in Chinese). [Link]

ii
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

List of Figures
1 Difference between statistics and probability . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Insignificant interaction vs. significant interaction . . . . . . . . . . . . . . . . . . . . . . . . 15
3 Profile plot of a two-way ANOVA with an interaction effect . . . . . . . . . . . . . . . . . . 15
4 Projection in linear regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
5 PCA vs. OLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
6 Orthogonal rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
7 Oblique rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
8 Maximum likelihood rule for univariate normal distributions . . . . . . . . . . . . . . . . . . 30
9 LDA decision boundary for bivariate normal distributions . . . . . . . . . . . . . . . . . . . 31
10 Fisher’s linear discriminant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

List of Tables
1 One-way MANOVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2 Distribution of Wilks’ Lambda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Two-way balanced MANOVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4 Confusion matrix in binary classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Get your facts first, and then you can distort them as much as you please. — Mark Twain

iii
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§1 Review of Linear Algebra (2019/9/10)

For any x, y ∈ Rn , their inner product is
X
n
hx, yi = x′ y = xk yk .
k=1

The projection of x on y 6= 0n is
hx, yi
projy x = y.
hy, yi
Consider the linear model
Y = Xβ + ε,
where X = (x1 , . . . , xp ) ∈ R n×p
is called the design matrix. Usually, n is # of observations and p is # of
variables. We are often faced with trade-off between parsimony and goodness of fit.
Choosing the maximal linearly independent system of xk ’s leads to dimension reduction. Recall
that rank (X) is the maximum # of linearly independent columns/rows, i.e., the dimension of the column
space Col(X) or the row space Col(X ′ ).
Given a symmetric A ∈ Rn×n , the corresponding quadratic form is

x ∈ Rn 7→ x′ Ax ∈ R.

The matrix A is said to be positive semi-definite if x′ Ax ≥ 0, ∀x ∈ Rn ; furthermore, A is said to be

positive definite if x′ Ax > 0 whenever x 6= 0n . Clearly, X ′ X is non-negative definite.
The eigen-decomposition of symmetric A ∈ Rn×n is
X
n
A= λi ei e′i = (e1 , . . . , en ) diag(λ1 , . . . , λn ) (e1 , . . . , en )′ ,
i=1

where λi ’s are eigenvalues and ei ’s are orthonormal eigenvectorsi) . Note that ei e′i = projei .
We say that Q ∈ Rn×n is an orthogonal matrix if Q′ Q = In . Any orthogonal matrix induces a map
preserving inner product, i.e.,
hQx, Qyi = hx, yi , ∀x, y ∈ Rn .
Write A = QΛQ′ with Q = (e1 , . . . , en ) and Λ = diag(λ1 , . . . , λn ). If y = Q′ x, then

X
n
x′ Ax = y ′ Λy = λk yk2 ,
k=1

which is closely related to ellipsoid when A is positive definite. For any analytic function f : R → R, we
can define
Xn
f (A) := f (λi )ei e′i .
i=1

Any orthogonal matrix of order 2 has the form

 
cos θ − sin θ
Q= , θ ∈ [0, 2π).
sin θ cos θ

One can justify this fact by a few calculations using orthonormal relations.

§2 Review of Linear Algebra Cont’d (2019/9/17)

Given a positive definite A ∈ R2×2 , write its spectral decomposition as A = QΛQ′ , where
   
λ1 cos θ − sin θ
Λ = diag(λ1 , λ2 ) = , Q = (e1 , e2 ) = .
λ2 sin θ cos θ
i) cf. [Link]

1
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

If y = Q′ x, then x′ Ax = y ′ Λy. Note that

   
kxk cos(θx ) kxk cos(θx − θ)
x= =⇒ y = ,
kxk sin(θx ) kxk sin(θx − θ)
which is explained as a rotation. Henceforth,

x ∈ R2 : x′ Ax = c2
is an ellipse with the orientation defined by e1 , e2 and length √cλ , √cλ .
1 2
Next, we explore orthogonalization for X ∈ Rn×p , i.e., to find Z ∈ Rn×p with orthogonal columns
such that Col(Z) = Col(X). The so called Gram-Schmidt process makes itii) . To be explicit,
z1 = x1 , z2 = x2 − projz1 x2 ,
..
.
zk = xk − projspan(z1 ,...,zk−1 ) xk , k ≤ p.
| {z }
∑
= k−1
i=1 projz xk i

It follows immediately that

X
k
xk = projzi xk , 1 ≤ k ≤ p,
i=1
which gives rise to X = ZΓ with upper-triangular Γ = (γkℓ ) ∈ Rp×p , where


 1, k = ℓ;
γkℓ = projzk xℓ , k < ℓ;


0, k > ℓ.
Let D = diag(kz1 k , . . . , kzp k), then
Q = ZD−1 , R = DΓ =⇒ X = QR.
The QR decomposition involves Q ∈ Rn×p with orthonormal columns and upper-triangular R ∈ Rp×p .
In linear regression (see Lecture §13), the hat matrix
PX = X(X ′ X)−1 X ′
is actually the orthogonal projection onto Col(X). Using the QR decomposition, we obtain
PX = QR(R′ Q′ QR)−1 R′ Q′ = QQ′ .
Write Q = (q1 , . . . , qp ), then
X
p X
p
ŷ = PX y = qk qk′ y = hy, qk i qk .
k=1 k=1
The above results can be generalized for any separable Hilbert space, where we need apply some cut-off to
pre-specified basis (e.g., Fourier basis, B-splines, wavelets, and other data-driven ones).
The singular value decomposition (SVD) is, sayiii) ,
rank(X)
X
′
X = U ΣV = σk uk vk′
k=1

such that Xvk = σk uk and X ′ uk = σk vk , where σk > 0 decreases in k. Note that

rank(X)
X rank(X)
X
X ′ X = V Σ2 V ′ = σk2 vk vk′ =⇒ PX = U U ′ =⇒ ŷ = PX y = hy, uk i uk .
k=1 k=1

From the perspective of matrix completion, SVD provides the best low-rank approximation in that
X
m
σk uk vk′ = arg min tr ((X − Y )(X − Y )′ ) , ∀m ≤ rank (X) .
k=1 Y ∈Rn×p : rank(Y )≤m
ii) cf. [Link]
iii) cf. [Link]

2
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§3 Prediction (2019/9/19)
Early statistics focused on inference. Nowadays, people are more concerned about prediction accuracy. The
AIC and BIC are frequently used for model selectioniv) . Generally, given data (X, Y ), we are longing for

min EL(Ŷ , Y ) = E[E[L(Ŷ , Y )|X]],

where L stands for loss function such as L2 loss (risk) and 0-1 loss (cost).
• Let Ŷ = f (X) and L(ŷ, y) = |ŷ − y|2 , then
   
2
E[L(Ŷ , Y )|X] = E (Y − f (X))2 X = E (Y − E[Y |X])2 X + E[Y |X] − f (X) .
MSE Var Bias

The best predictor

fˆ(x) = E[Y |X = x]
for squared loss is appropriate for both fixed design and random/conditional design.
• In classification problem, Y ∈ {1, 2, · · · , k} is categorical and we often adopt 0-1 loss L(ŷ, y) = 1[ŷ̸=y] .
Let Ŷ = g(X), then
X
E[L(Ŷ , Y )|X = x] = P(Y = y|X = x) = 1 − P(Y = g(x)|X = x).
y̸=g(x)

Maximizing P(Y = g(x)|X = x) gives the naïve Bayesian classifier ĝ(x) — the most probable class
given x. We may treat y as a vector of dummy variables

yj = 1[y=j] , j ∈ {1, · · · , k},

which add up to 1 and thus the degree of freedom is k − 1. Henceforth,

ĝ(x) = arg max E[Yj |X = x].

j∈{1,··· ,k}

 (

 continuous

numerical
 discrete (counting)
Most data fall into one of two groups: (

 nominal


categorical
ordinal
Given data (xi , yi )1≤i≤n , one may estimate
1 X
E[Y |X = x] ≈ yi ,
#{i : xi = x} x =x
i

where {1 ≤ i ≤ n : xi = x} is the effective sample. If we replace it with the neighborhood Nk (x), referred
to as k-nearest neighbor (kNN), the kernel estimation (or kNN classifier) will come up. The case k = 1
results in perfectv) fit, where no data reduction exists, and hereby the variance tends to explode even though
the bias may vanish. The case when k is large will see larger bias and smaller variance. Anyway, we are
desired to obtain suitable effective sample size, which could be rather difficult in real world.

§4 Prediction Cont’d & Multivariate Data (2019/9/24)

The curse of dimensionality turns out to be more
P significant when
R the model considered is nonparametric
rather than parametric. The kNN estimator k1 xi ∈Nk (x) yi has Nk (x) fX (u) du ∝ fX (x)hx as its effective
sample size, which often takes the form A
nhx + Bh4x of MSE = Var +Bias2 due to the smoothness of the
iv) cf. [Link] and Lecture §15
v) ironic?

3
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

density function fX ∈ C 2 , where hx is the neighborhood size (bandwidth). Trade-off between variance and
bias provides the minimizer h∗ ∝ n−1/5 and the minimum MSE(h∗ ) ∝ n−4/5 . Note that the parametric
model has MSE ∝ n−1 . Generally, one may encounter
A
MSE(hx ) = + Bhbx ,
nhdx
where d denotes the dimension and b is twice the order of differentiability of fX . The optimal results are
h∗ ∝ n−1/(d+b) and MSE(h∗ ) ∝ n−b/(d+b) . For efficiency, we should always require h → 0 and nhd → ∞ as
the sample size n tends to infinity. See also exercise 5 in homework 1 for example.
The object to be studied is data matrix, which consists of observations and variables, denoted by
 ′
X1
 .. 
X = (xij )1≤i≤n
1≤j≤p =  .  = (X(1) , . . . , X(p) ).
Xn′
To check data quality, one may plot histogram, boxplot, pairwise scatter plots, etc. The concept of outlier
is subtle, which suggests some alert. Exploratory data analysis helps finding structures in data —
pattern, describe, visualize, PCA, cluster, etc. These are nonparametric (model-free) approaches, also
called unsupervised learning. Confirmatory data analysis is model-based, including statistical inference
and prediction. In other words, we assume

data ∼ model + error.

The model is wished both parsimony (as simple as possible, causing larger bias) and fidelity (adequate to
describe data, causing larger variance). It’s worth noting thatvi)
all models are wrong, but some are useful. — George E. P. Box

§5 Descriptive Statistics & Normal Distribution (2019/10/8)

Given data matrix X = (xij )1≤i≤n
1≤j≤p , the sample mean is defined to be

1X
n
X̄ = (x̄1 , . . . , x̄p )′ , x̄j = xij for 1 ≤ j ≤ p.
n i=1

The sample covariance is defined to be

1 X
n
S = (sjk )1≤j,k≤p , sjk = (xij − x̄j )(xik − x̄k ).
n − 1 i=1

The sample correlation is defined to be

√
R = (rjk )1≤j,k≤p , rjk = sjk / sjj skk ∈ [−1, 1].
√
Consider standardized data x∗ij = (xij − x̄j )/ sjj , its sample covariance is exactly S ∗ = R = D−1 SD−1 ,
√ √

where D = diag s11 , . . . , spp . The sample correlation characterizes the linear association between
parameters. In linear regression with one regressor, the sample covariance appears naturally. Recall that
the Cauchy-Schwarz inequality states that
√ √
2
a′ A Bb ≤ (a′ Aa)(b′ Bb).

For sampling statistics,

Plet X1 , . . . , Xn be i.i.d.
Pn random vectors with mean µ and covariance
Pn matrix Σ. The
n
sample mean X̄ = n1 i=1 Xi has mean n1 i=1 EXi = µ and covariance matrix n12 i=1 Var(Xi ) = n1 Σ.
The denominator (n − 1) of the sample covariance

1 X
n
S= (Xi − X̄)(Xi − X̄)′
n − 1 i=1
vi) see [Link] for more fun.

4
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

guarantees the unbiasedness ES = Σ. PnTo see this, one can plug EXi Xi′ = Σ + µµ′ and EX̄ X̄ ′ = n1 Σ + µµ′
′ ′
into the expectation of (n − 1)S = P i=1 Xi Xi − nX̄ X̄ , and thus the corresponding degree of freedom
n
(d.o.f.) is n − 1 since one constraint i=1 (Xi − X̄) = 0 is imposed. Also, the MLE (see Lecture §7)
1
Σ̂n = (Xi − X̄)(Xi − X̄)′
n
for Σ in Gaussian setting is biased. In matrix language, X̄ = n1 X ′ 1n and (n − 1)S = X ′ (In − n1 1n 1′n )X.
Denote the centering matrix by
H = In − n1 1n 1′n ⊥ n1 1n 1′n .
Clearly H = H ′ = H 2 is symmetric and idempotent, and thus rank (H) = tr (H) = n − 1.
Probability theory deduces the underlying structure of observed data/samples. Statistical inference
aims to learn unknown (or ad hoc) population/distribution/parameters from the observations.

Figure 1: Difference between statistics and probabilityvii)

In multivariate framework, we often specify notations in this manner:

• EX = µ = (µj )1≤j≤p ;
• Var(X) = Σ = (σjk )1≤j,k≤p ;
√
• ρ = (ρjk )1≤j,k≤p with ρjk = σjk / σjj σkk ; and
• Cov(X, Y ) = E(X − EX)(Y − EY )′ = ΣXY . Note that Cov(AX + a, BY + b) = AΣXY B ′ .
To justify the normal distribution, the central limit √theorem (CLT) states that if X1 , X2 , . . . are
i.i.d. random variables with mean µ and variance σ 2 , then n(X̄n − µ)/σ → N (0, 1) in distribution. For
non-degenerate multivariate normal distribution Np (µ, Σ), the probability density function (p.d.f.) is
 
1 1 ′ −1
f (x) = exp − (x − µ) Σ (x − µ) , x ∈ Rp .
(2π)p/2 (det Σ)1/2 2

The crucial part is Σ−1 , called the precision matrix, which is important in graphical models. Generally,
independence implies uncorrelatedness, but not vice versa. If random p-vector X1 and random q-vector X2
are jointly normal, then Cov(X1 , X2 ) = 0 if and only if X1 ⊥⊥ X2 . The conditional distribution is more
complicated, involving some matrix holing tricksviii) . Some frequently used properties are as follows:
i) X ∼ Np (µ, Σ) =⇒ AX + b ∼ Nq (Aµ + b, AΣA′ ) for A ∈ Rq×p and b ∈ Rq ;
ii) if a′ X ∼ N (a′ µ, a′ Σa) for any a ∈ Rp , then X ∼ Np (µ, Σ);

√
iii) X ∼ Np (µ, Σ) ⇐⇒ E exp(it′ X) = exp it′ µ − 12 t′ Σt , t ∈ Rp , where i = −1.
vii) from Herman Bennett’s 14.30 Introduction to Statistical Method in Economics. Spring 2006. [Link]
viii) cf. [Link]

5
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§6 Multivariate Normal Distribution Cont’d & MLE (2019/10/15)

The distribution of any random vector X corresponds uniquely to all one-dimensional
√ distributions of its
linear functions a′ X. To see this,nthe characteristic
o function ϕ X (t) = E exp( −1t ′
X) applies.
|X−µ|
For inference, we consider P ≤c
for random variable X ∼ N (µ, σ 2 ) and constant c ∈ R+ .
σ

If instead X ∼ Np (µ, Σ) is a random vector, the region considered  becomes (X − µ)′ Σ−1 (X − µ) ≤ c2 ,
which correspond to contours with constant density x ∈ R : fX (x) ≥ c20 . Write Σ = QΛQ′ with
p

Q = (e1 , . . . , ep ) orthogonal and Λ = diag(λ1 , . . . , λp ). Let Y = (Y1 , . . . , Yp )′ = Q′ (X − µ) ∼ Np (0p , Λ),

then X
(X − µ)′ Σ−1 (X − µ) = Y ′ Λ−1 Y = Yj2 /λj ∼ χ2p ,
geometricallycalling
 on 
some ellipsoid.
 
X1 µ1 Σ11 Σ12
Suppose ∼N , , then X1 − Σ12 Σ−122 X2 is independent of X2 , and thus
X2 µ2 Σ21 Σ22


X1 |X2 ∼ N µ1 + Σ12 Σ−1 −1
22 (X2 − µ2 ), Σ11 − Σ12 Σ22 Σ21 .

Note that E[X1 |X2 ] = µ1 + Σ12 Σ−1 22 (X2 − µ2 ) minimizes E|Z − X1 | among all Z = g(X2 ), i.e., the
2

conditional expectation is the population version of the best linear prediction, compared to the sample
−1
version Ŷ = Ȳ + SY X SXX (X − X̄) in linear regressionix) (see Lecture §13).
We now explore the maximum likelihood estimation (MLE) for θ ∈ Θ given a sample x1 , . . . , xn . The
likelihood function is the joint density
Y
n
Ln (θ|x1 , . . . , xn ) = f (x1 , . . . , xn |θ) = f (xi |θ).
i=1

Hence, the log-likelihood function is

X
n
ℓn (θ) = log Ln (θ|x1 , . . . , xn ) = log f (xi |θ).
i=1

To justify the MLE, suppose Xi ’s are i.i.d. with p.d.f. f0 = f (·|θ0 ). Shannon’s entropy inequality states
that
E log f (X) ≤ E log f0 (X)
for any p.d.f. f = f (·|θ), with equality holding if and only if f = f0 almost everywhere. Note that the
identifiablity means
θ1 6= θ2 ⇐⇒ f (·|θ1 ) 6= f (·|θ2 ).
By SLLN, n1 ℓn (θ) → Elog f (X|θ) almost surely as n → ∞, so the maximizer θ̂nMLE of ℓn (·) seems reasonable
to estimate the true parameter θ0 . Under some mild assumptions, θ̂nMLE is consistent, as desired.

§7 Maximum Likelihood Estimation Cont’d (2019/10/17)

In multivariate normal setting,
Yn   
1 1 ′ −1
Ln (µ, Σ) = exp − (xi − µ) Σ (xi − µ)
i=1
(2π)p/2 (det Σ)1/2 2
( )
1 1X
n n n o
′ −1 ′ −1
= exp − (x i − x̄) Σ (x i − x̄) exp − (x̄ − µ) Σ (x̄ − µ) ,
(2π)np/2 (det Σ)n/2 2 i=1 2
Pn
so µ̂MLE = x̄ and Σ̂MLE = n1 i=1 (xi − x̄)(xi − x̄)′ (exercise 2 in homework 2) with Bias(Σ̂MLE ) = − n1 Σ.
They are asymptotically consistent, minimal sufficient and completex) .
ix) regression toward the mean is the phenomenon that arises if a random variable is extreme on its first measurement but

closer to the mean or average on its second measurement and if it is extreme on its second measurement but closer to the
average on its first.
x) a statistic T = T (X) is said to be complete for θ ∈ Θ if, E f (T ) = 0 for all θ ∈ Θ implies f (T ) = 0 a.s., and thus T
θ
appears to be successful in reducing the data.

6
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

To derive MLEs, we can use matrix differentiationxi) as well. For symmetric positive definite A = (aij ),

• ∂(x′ Ax)/∂x = 2Ax;

• ∂ tr (AB) /∂A = 2B − diag(B);
(
2Aij / det(A), i 6= j

• ∂ log det(A)/∂aij = and thus ∂ log det(A)/∂A = 2A−1 − diag A−1 .
Aii / det(A), i=j

Note that X̄ ∼ Np (µ, n1 Σ) is independent of nΣ̂MLE ∼ Wp (Σ, n − 1). See Lecture §9 for the Wishart
distribution. As the sample size n → ∞, the following large sample properties hold:

• X̄ →P µ, Σ̂ →P Σ;
√
• n(X̄ − µ) →d Np (0p , Σ);

• n(X̄ − µ)′ Σ̂−1 (X̄ − µ) →d χ2p .

Generally, the derivative of ℓn (θ) is called the score function, denoted by Sn (θ). Under some regularity
conditions,
  Z Z
∂ 1∂f (x|θ)  ∂
Eθ S(θ) = Eθ log f (X|θ) =  
f (x|θ) dx = f (x|θ) dx = 0p ,
∂θ f (x|θ)
 ∂θ ∂θ
| {z }
=1

and (check!)  
∂
−Eθ Sn (θ) = Eθ [Sn (θ)Sn (θ)′ ] = Varθ (Sn (θ))
∂θ′
is exactly the Fisher information In (θ). Note that In (θ) = nI1 (θ). As n → ∞,
√

• n(θ̂MLE − θ) →d Np 0p , I1 (θ)−1 ;

• n(θ̂MLE − θ)′ I1 (θ̂MLE )(θ̂MLE − θ) →d χ2p .

∂
Sketch of Proof. Denote by Hn (θ) = ∂θ ′ Sn (θ) the Hessian matrix of ℓn (θ). Then
√
−1
0p = Sn (θ̂MLE ) ≈ Sn (θ) + Hn (θ)(θ̂MLE − θ) =⇒ n(θ̂MLE − θ) ≈ − n1 Hn (θ) √1 Sn (θ),
n

where − n1 Hn (θ) →P I1 (θ) and √1 Sn (θ)

n
→d Np (0p , I1 (θ)).

§8 Diagnostics for Normality & Matrix Distribution (2019/10/22)

For univariate data, ordered as x(1) ≤ · · · ≤ x(n) called sample quantiles, the most useful diagnostic tool is
the Q-Q plot. Let q(1) , . . . , q(n) be theoretical quantiles of N (0, 1), i.e.,

i − 21
Φ(q(i) ) = , 1 ≤ i ≤ n.
n
Plot x(i) against q(i) , and then compare it to a line. If Xi ∼ Np (µ, Σ) are multivariate, then

(Xi − µ)′ Σ−1 (Xi − µ) ∼ χ2p .

By large sample property X̄ →P µ and S →P Σ, as n  p,

◦
d2i := (Xi − X̄)′ S −1 (Xi − X̄) ∼ χ2p .

Sort d2i to obtain sample quantiles d2(i) , and plot them against theoretical quantiles of χ2p .
xi) cf. [Link]

7
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

To detect outliers and clean data, one may plot histogram or dot plot for each variable, and pairwise
scatterplot for each pair. Sometimes we prefer standardized variable
xik − x̄k
zik = √ , 1 ≤ i ≤ n, 1 ≤ k ≤ p.
skk

Note that judging outliers should take both n and p into consideration. The empirical 3σ rulexii) is not
always appropriate.
Transforms are used to normalize univariate data, including
√
• y ∈ N 7→ y for counts;
 
• logit : y ∈ [0, 1] 7→ log 1−yy
for proportions;
 
• r ∈ [−1, 1] 7→ log 1−r 1+r
for correlation coefficients;
(
(xλ − 1)/λ, λ > 0
• Box-Cox transform x 7→ x = (λ)
, where λ is a tuning parameter to maximize
log(x), λ=0
 P  Pn
n (λ) 2
ℓ(λ) = − n2 log n1 i=1 xi − x(λ) + (λ − 1) i=1 log(xi ), or equivalently, to minimize the sample
(λ) xλi −1
variance of yi := ∏
λ( n (λ−1)/n . For multivariate data, we often choose λk for Xik separately,
j=1 xj )
pretending the components are independent.
i.i.d.
Next, multivariate data should be tackled. If X1 , . . . , Xn ∼ Np (µ, Σ), then it is said that
 ′
X1
 .. 
X =  .  = (X(1) , . . . , X(p) ) ∼ MNn×p (1n µ′ , In , Σ)
Xn′
is a normal matrix with p.d.f.
 
−n/2 1 −1 ′ ′ ′


f (X) = det(2πΣ) exp − tr Σ (X − 1n µ ) (X − 1n µ ) .
2

Using vectorization and Kronecker productxiii) , we have

 
X(1)
 
vec(X) =  ...  ∼ Nnp (µ ⊗ 1n , Σ ⊗ In ).
X(p)
Several frequently used properties are linearity, mixed-product property, and vectorization formula
vex(AXB) = (B ′ ⊗ A) vec(X).
It follows that (see exercise 4 in homework 2)
vex(AXB) ∼ Nmq ((B ′ µ) ⊗ (A1n ), (B ′ ΣB) ⊗ (AA′ ))
for A ∈ Rm×n and B ∈ Rp×q . If µ = 0p and Q ∈ Rn×n is orthogonal, then
vec(QX) =d vec(X) ∼ Nnp (0p ⊗ 1n , Σ ⊗ In ).
Let Y = AXB and Z = CXD, then
Cov(vec(Y ), vec(Z)) = (B ′ ΣD) ⊗ (AC ′ ),
so Y ⊥⊥ Z if and only if AC ′ = 0 or B ′ ΣD = 0. As an application, to see µ̂ = 1 ′
n X 1n independent of
Σ̂ = n1 X ′ (In − n1 1n 1′n )X, denoting H = In − n1 1n 1′n , it suffices that
Cov(vec(1′n X), vec(HX)) = Σ ⊗ (1′n H) = 0.
xii) cf. [Link]
xiii) cf. [Link]

8
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§9 Cochran’s Theorem & Hypothesis Testing (2019/10/29)

The celebrated Cochran’s theoremxiv) concerns simultaneous diagonalizationxv) to idempotent matrices.
Pn Pk
Let X = (X1 , . . . , Xn )′ ∼ Nn (0n , σ 2 In ). Suppose i=1 Xi2 = j=1 Qj with Qj = X ′ Aj X
Pk Pk
for some symmetric Aj ∈ Rn×n of rank rj . That is, In = j=1 Aj . If j=1 rj = n, then there
exists some orthogonal C ∈ Rn×n such that C ′ X = (Y1 , . . . , Yn )′ satisfies
X
r1 +···+rj−1 +rj
Qj = Yi2 , 1 ≤ j ≤ k.
i=1+r1 +···+rj−1

Note that Qj ∼ σ 2 χ2rj (j = 1, · · · , k) are independent.

For a multivariate version, some acquaintance of the Wishart distribution is prerequisite. Consider the
i.i.d. Pn
data matrix X = (X1 , . . . , Xn )′ with Xi ∼ Np (0p , Σ), it is said that X ′ X = i=1 Xi Xi′ ∼ Wp (Σ, n).
• W1 (σ 2 , n) = σ 2 χ2n when p = 1;
• BM B ′ ∼ Wq (BΣB ′ , n) for M ∼ Wp (Σ, n) and B ∈ Rq×p ;
• a′ M a ∼ (a′ Σa)χ2n for M ∼ Wp (Σ, n) and a ∈ Rp ;
• M1 + M2 ∼ Wp (Σ, n1 + n2 ) for independent Mi ∼ Wp (Σ, ni ), i = 1, 2.
i.i.d.
Let X = (X1 , . . . , Xn )′ with Xi ∼ Np (0p , Σ), and A1 , . . . , Ak , A, B ∈ Rn×n be symmetric, then
(a) X ′ AX is a λ(A)-weighted sum of independent Wp (Σ, 1) matrices;
(b) X ′ AX ∼ Wp (Σ, r) if and only if A is idempotent of rank r;
(c) X ′ AX ⊥⊥ X ′ BX if and only if AB = 0.
Pk Pk Pk
If X ′ X = j=1 X ′ Aj X (i.e., In = j=1 Aj ) and j=1 rj = n, where rj = rank (Aj ) for 1 ≤ j ≤ k, then
X ′ Aj X are independent and distributed as Wp (Σ, rj ) matrices, respectively.
Since H = In − n1 1n 1′n has a spectral decomposition H = CΛC ′ such that Λ = diag(In−1 , 0) and
C = (e1 , . . . , en−1 , en ) with en = √1n 1n , denoting Y = (Y1 , . . . , Yn )′ = C ′ X, it follows that

X
n−1
nΣ̂ = (n − 1)S = X ′ HX = Y ′ ΛY = Yi Yi′ ∼ Wp (Σ, n − 1).
i=1

From the perspective of Cochran’s theorem,

X ′ X = X ′ (In − n1 1n 1′n )X + X ′ ( n1 1n 1′n )X,
where In − n1 1n 1′n and 1 ′
n 1n 1n are idempotent and of rank n − 1 and 1, respectively.
i.i.d.
Next, we explore inferences for the mean of normal distribution. For univariate Xi ∼ N (µ, σ 2 ) with
H0
σ > 0 unknown, T = X̄−µ √ 0 is often employed to test H0 : µ = µ0 , since T ∼ tn−1 no matter which value
s/ n
µ0 takes.
h To reject H0 at level α, one can adopt
i either the criterion {|T | > tn−1 (α/2)}, or equivalently, that
H
µ0 ∈
/ X̄ − √s tn−1 (α/2), X̄
n
+ √s tn−1 (α/2)
n
. Note that T 2 = n(X̄ − µ0 )(s2 )−1 (X̄ − µ0 ) ∼0 F1,n−1 .
In multivariate setting, if d ∼ Np (0p , Ip ) and M ∼ Wp (Ip , m) are independent, then
mp
md′ M−1 d ∼ T 2 (p, m) = Fp,m−p+1 ,
m−p+1
i.i.d. √
called Hotelling’s T-squared distribution. Let Xi ∼ Np (µ, Σ). Now that n(X̄ − µ) ∼ Np (0p , Σ)
and (n − 1)S ∼ Wp (Σ, n − 1) are independent, it’s immediate that
√ −1 √
T 2 = n(X̄ − µ)′ S−1 (X̄ − µ) = (n − 1) n(X̄ − µ)′ ((n − 1)S) n(X̄ − µ) ∼ T 2 (p, n − 1).
Hotelling’s T 2 statistic is clearly a generalizationnof Student’s t-statistic for testing H0 : µ = µ0o
. Hence,
′ −1 (n−1)p
a 100(1 − α)% confidence region/ellipsoid for µ is µ : n(X̄ − µ) S (X̄ − µ) ≤ n−p Fp,n−p (α) . For a
◦
large sample (n  p), we have n(X̄ − µ)′ S−1 (X̄ − µ) ∼ χ2p , so (n−1)p
n−p Fp,n−p (α) can be replaced by χ2p (α).
xiv) cf. [Link]
xv) cf. [Link]

9
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§10 Hypothesis Testing Cont’d (2019/10/31)

Generally, the likelihood-ratio test statistic for H0 : θ ∈ Θ0 vs. H1 : θ ∈ Θ \ Θ0 is
maxθ∈Θ0 L(θ)
Λ= ,
maxθ∈Θ L(θ)
with d.o.f.= dim(Θ) − dim(Θ0 ). Wilks’ theorem asserts that

−2 log(Λ) →d χ2dim(Θ)−dim(Θ0 ) .

The rejection region has the form {Λ < c} for some c.

i.i.d.
Consider Xi ∼ Np (µ, Σ) and H0 : µ = µ0 ↔ H1 : µ 6= µ0 . The likelihood function is
 
1

L(µ, Σ) = (2π)−np/2 det(Σ)−n/2 exp − tr Σ−1 (X − 1n µ′ )′ (X − 1n µ′ ) .
2

Hence (see exercise 2 in homework 2)

max L(µ, Σ) = (2π)−np/2 det(Σ̂)−n/2 e−np/2 ,

µ,Σ

Pn
where Σ̂ = 1
n i=1 (Xi − X̄)(Xi − X̄)′ , and

max L(µ0 , Σ) = (2π)−np/2 det(Σ̂0 )−n/2 e−np/2 ,

Σ
Pn
where Σ̂0 = 1
n i=1 (Xi − µ0 )(Xi − µ0 )′ is the constrained/restricted MLE. It follows that
Pn ′


det(Σ̂) det i=1 (Xi − X̄)(Xi − X̄)
Λ 2/n
= = Pn
det(Σ̂0 ) det ( i=1 (Xi − µ0 )(Xi − µ0 )′ )

is equivalent to Hotelling’s T 2 in the sense of monotone bijection, since

X
n X
n
(Xi − µ0 )(Xi − µ0 )′ = (Xi − X̄)(Xi − X̄)′ + n(X̄ − µ0 )(X̄ − µ0 )′ ,
i=1 i=1

the decomposition of the sums of squares, implies (see exercise 1 in homework 1)

−1 1
Λ−2/n = 1 + n(X̄ − µ0 )′ ((n − 1)S) (X̄ − µ0 ) = 1 + T 2.
n−1

To test H0 : θ = θ0 in other ways, denoting by θ̂ = arg maxθ∈Θ L(θ) the unconstrained MLE,
• the Wald test employs

d θ̂)−1 (θ̂ − θ0 ) = (θ̂ − θ0 )′ In (θ̂)(θ̂ − θ0 ) −−→ χ2

W = (θ̂ − θ0 )′ Var(
d
dim(Θ) ,
H0

noting that Varθ (θ̂) ≈ In (θ)−1 (see Lecture §7); and

• the score test (a.k.a. the Lagrange multiplier test, particularly in econometrics) uses

Un (θ0 )′ In (θ0 )−1 Un (θ0 ) −−→ χ2dim(Θ) ,

d
H0

∂ log L(θ)
where U (θ) = ∂θ is the score.
The three tests above are equivalent for large samples. An advantage of the Wald test over the other two
is that it only requires the estimation of the unrestricted model, which lowers the computational burden.
The score test only requires the estimation of the likelihood function under the null hypothesis, and thus
it is less specific than the other two tests about the precise nature of the alternative hypothesisxvi) .
xvi) a.k.a. the research hypothesis, being more important and meaningful than the null hypothesis (benchmark).

10
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

The multiple testing problem occurs when one considers a set of statistical inferences simultaneously.
i.i.d.
Let X1 , . . . , Xn ∼ Np (µ, Σ) with Σ unknown. For any single a ∈ Rp ,
√
n(a′ X̄ − a′ µ)
T (a) = √ ∼ tn−1
a′ Sa
q
′
gives a 100(1 − α)% confidence interval a′ X̄ ± tn−1 ( α2 ) a nSa for a′ µ. To obtain some c > 0 such that
n o
2
P |T (a)| ≤ c2 , ∀a ∈ Rp = 1 − α,

note that
2
n a′ (X̄ − µ)
= n(X̄ − µ)′ S −1 (X̄ − µ) = T 2
2
max |T (a)| = max
a a a′ Sa
by Cauchy-Schwarz inequality, so

\  q !
′ ′ (n−1)p a′ Sa
P a µ ∈ a X̄ ± n−p Fp,n−p (α) n ≥ 1 − α.
a∈Rp

p q
(n−1)p
In other words, replacing tn−1 ( α2 ) = F1,n−1 (α) with n−p Fp,n−p (α) yields simultaneous confidence
intervals for all possible linear functions of µ.
When a small number of statements are of interest, the Bonferroni correction can be used to adjust
confidence intervals. If P(Ck ) = 1 − αk for 1 ≤ k ≤ m, then
X
P(∩Ck ) = 1 − P(∪Ck{ ) ≥ 1 − αk ,

and thus P(∩Ck ) ≥ 1 − α as long as αk = α/m. For example,

\m  q !
′ ′ a′ Sa
P ak µ ∈ ak X̄ ± tn−1 ( 2m )
α
n ≥ 1 − α, ∀a1 , . . . , am ∈ Rp .
k=1

§11 Inferences about Multi-Population Means (2019/11/5)

To compare two treatments, the same units
 i = 1, 2, · · · 
 , n are assigned
 to both
 so that the comparison
X1i i.i.d. µ1 Σ11 Σ12
makes sense. Formally, consider Xi = ∼ N2p , , and H0 : µ1 = µ2 is of
X2i µ2 Σ21 Σ22
interest. Let Di = X1i − X2i and µd = µ1 − µ2 . It suffices to test H0 : µd = 0p using

−1 (n − 1)p
= nD̄′ SD
H
2
TD D̄ ∼0 T 2 (p, n − 1) = Fp,n−p .
n−p

If C ∈ Rq×p has full row rank q ≤ p, then the statistic used to test the constraint

H0 : Cµ = ν0 (= 0q )
i.i.d.
on the population mean µ ∈ Rp of Xi ∼ Np (µ, Σ) is

(n − 1)q
T 2 = n(C X̄ − ν0 )′ (CSC ′ )−1 (C X̄ − ν0 ) ∼0
H
Fq,n−q ,
n−q

since CXi ∼ Nq (Cµ, CΣC ′ ). Note that T 2 does not depend on the particular choice of C, i.e., T 2 computed
with C̃ instead of C yields the same result provided that C̃ = QC for some invertible Q ∈ Rq×q .
A contrast is a linear combination of variables whose coefficients add up to zero. Attention is usually
centered on contrasts when comparing different treatments. The example at the beginning introduces
C = (Ip , −Ip ) as a contrast matrix xvii) , from which Di = CXi , D̄ = C X̄ and SD = CSX C ′ are derived.
xvii) cf. [Link]

11
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

In repeated measurements, subjects i are assumed independent, but treatments j within each subject are
dependent, which differs from that in one-way univariate ANOVA (see also Lecture §12). The responses
are
xij = µj + εij , i = 1, · · · , n; j = 1, · · · , p.
i.i.d.
Here εi = (εi1 , . . . , εip )′ ∼ Np (0p , Σ) are random errors. In order to test H0 : µ1 = µ2 = · · · = µp , one
may take    
1 −1 1 −1
 1 −1  1 −1 
   
C= .. ..  or  . ..  ∈ R(p−1)×p ,
 . .   .. . 
1 −1 1 −1
each of which gives
(n − 1)(p − 1)
T 2 = n(C X̄)′ (CSC ′ )−1 (C X̄) ∼0
H
Fp−1,n−p+1 .
n−p+1
If we write µj = µ + αj , then (µ, α1 , . . . , αp )′ ∈ Rp+1 is not identifiable/estimable,
Pp but H0 : α1 = · · · = αp
is indeed estimable and thus can be tested. It’s often imposed that j=1 αj = 0.
i.i.d.
Next, we will compare means from two populations. Suppose X11 , . . . , X1n1 ∼ Np (µ1 , Σ1 ) and
i.i.d.
X21 , . . . , X2n2 ∼ Np (µ2 , Σ2 ) are two independent samples. The goal is to test H0 : µ1 = µ2 . As-
sume Σ1 = Σ2 . Note that the first-order information is meaningful only if the second-order information
keeps almost the same. Denote
1 X
nk
X̄k = Xki
nk i=1
and
1 X
nk

Sk = (Xki − X̄k )(Xki − X̄k )′

nk − 1 i=1
for k = 1, 2. Let
(n1 − 1)S1 + (n2 − 1)S2
Spooled = ,
n1 + n2 − 2
whose d.o.f. is n1 + n2 − 2 since two sample means are computed separately. It holds that
  −1

′ 1 1
(n1 + n2 − 2)p
X̄1 − X̄2 − (µ1 − µ2 ) + Spooled X̄1 − X̄2 − (µ1 − µ2 ) ∼ Fp,n1 +n2 −p−1 .
n1 n2 n1 + n2 − p − 1
For large samples, having no need of Σ1 = Σ2 or normality,
 −1

′ 1 1
◦
X̄1 − X̄2 − (µ1 − µ2 ) S1 + S2 X̄1 − X̄2 − (µ1 − µ2 ) ∼ χ2p .
n1 n2
Profile analysis pertains to situations in which a series of treatments are administered to two or more
groups of subjects. It is assumed that (i) all measurements are expressed in similar scales; and (ii) all
subjects are independent. Three questions on the population means µ1 , µ2 ∈ Rp are formulated stepwise:
(1) Are the profiles parallel in the sense that µ1 − µ2 is constant?
(2) If (1) holds, are the profiles the same in the sense that µ1 = µ2 ?
(3) If (2) holds, are the treatments identical in the sense that µ,1 = · · · = µ,p ?
To that end, denoting the contrast matrix by
 
−1 1
 −1 1
 
C= ..  ∈ R
(p−1)×p
.. ,
 . . 
−1 1
we can carry out the sequential tests as follows:

12
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

(1) Test H01 : Cµ1 = Cµ2 using

  −1

′ 1 1
H (n1 + n2 − 2)(p − 1)
X̄1 − X̄2 C ′ + CSpooled C ′ C X̄1 − X̄2 ∼0 Fp−1,n1 +n2 −p .
n1 n2 n1 + n2 − p
If H01 is rejected, stop; otherwise, continue.
(2) Test H02 : 1′p µ1 = 1′p µ2 using
  −1

1 1
H
1′p X̄1 − X̄2 + 1′p Spooled 1p 1′p X̄1 − X̄2 ∼0 F1,n1 +n2 −2 = t2n1 +n2 −2 .
n1 n2
If H02 is rejected, stop; otherwise, continue.
(3) Test H03 : Cµ = 0p−1 using

(n1 + n2 − 1)(p − 1)
(n1 + n2 )X̄ ′ C ′ (CSC ′ )−1 C X̄ ∼0
H
Fp−1,n1 +n2 −p+1 ,
n1 + n2 − p + 1
where X̄ and S are based on all n1 + n2 observations.

§12 Multivariate Analysis of Variance (2019/11/12)

Usually, more than two treatments (e.g. drugs) are administered. In the one-way (single-factor) MANOVA
(multivariate analysis of variance) model,

Xℓr = µ + τℓ + εℓr (r = 1, · · · , nℓ ; ℓ = 1, · · · , g)
(observation) (overall mean) (the ℓth treatment effect) (error)

represents the p-dimensional response under the ℓth treatment in the rth replication, where the constraint
Pg i.i.d.
ℓ=1 nℓ τℓ = 0p is specified. The errors are assumed to be εℓr ∼ Np (0p , Σ), which only influences the
distribution of the test statistic.
ANOVA is a special case of linear models (see Lecture §13), and here
   
X11 ε11
 ..   .. 
 .       . 
  1n1 ⊗ Ip 1n1 ⊗ Ip µ  
X1n1  ε1n1 
  1n2 ⊗ Ip 1n2 ⊗ Ip     
 ..    τ1   .. 
 .  = .. ..   .  +   ,
   . .   ..   . 
 Xg1   εg1 
  1ng ⊗ Ip 1ng ⊗ Ip np×(g+1)p τg (g+1)p×1  
 .   . 
 ..   .. 
Xgng np×1 εgng np×1
Pg
where n = ℓ=1 nℓ is the total number of observations. Note that

Xℓr = X̄·· + (X̄ℓ· − X̄·· ) + (Xℓr − X̄ℓ· ),

=µ̂ =τ̂ℓ =ε̂ℓr
Pnℓ Pg Pnℓ
where X̄ℓ· = n1ℓ r=1 Xℓr and X̄·· = 1
n ℓ=1 r=1 Xℓr . Thereby, the sum of squares and cross products
(SSP) breaks up into
X
g X
nℓ X
g X
g X
nℓ
(Xℓr − X̄·· )(Xℓr − X̄·· )′ = nℓ (X̄ℓ· − X̄·· )(X̄ℓ· − X̄·· )′ + (Xℓr − X̄ℓ· )(Xℓr − X̄ℓ· )′
ℓ=1 r=1 ℓ=1 ℓ=1 r=1

since the cross terms add up to

X
g X
nℓ
(Xℓr − X̄ℓ· )(X̄ℓ· − X̄·· )′ = 0p×p .
ℓ=1 r=1
| {z }
=0p

13
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Source SSP d.o.f.

P
g
Treatment/between-group B= nℓ (X̄ℓ· − X̄·· )(X̄ℓ· − X̄·· )′ g−1
ℓ=1
P
g P nℓ
Residual/within-group W = (Xℓr − X̄ℓ· )(Xℓr − X̄ℓ· )′ n−g
ℓ=1 r=1
Pg Pnℓ
Total B+W = (Xℓr − X̄·· )(Xℓr − X̄·· )′ n−1
ℓ=1 r=1

Table 1: One-way MANOVA

Pg
We summarizePthe variations in table 1. The within SSP can be expressed as W = ℓ=1 (nℓ − 1)Sℓ ,
(Xℓr − X̄ℓ· )(Xℓr − X̄ℓ· )′ is the sample covariance matrix for the ℓth treatment, and
nℓ
where Sℓ = nℓ1−1 r=1
Pg
thus the corresponding d.o.f. is ℓ=1 (nℓ − 1) = n − g.
The hypothesis of no treatment effects, H0 : τ1 = · · · = τg = 0p , can be tested using Wilks’ Lambda
det(W )
Λ∗ = ,
det(B + W )
 n/2
det(Σ̂) 1
which is equivalent to Λ = det(Σ̂0 )
in the likelihood-ratio test, where Σ̂ = nW is the unrestricted
1 ∗
MLE and Σ̂0 = n (B + W ) is the restricted MLE. Therefore, we reject H0 if Λ < c for some c.

dim of variables no. of groups Sampling distribution for normal data

∗ H
p=1 g≥2 n−g 1−Λ
g−1 Λ∗ ∼0 Fg−1,n−g
√
∗ H
p=2 g≥2 n−g−1 1−
g−1
√ Λ
Λ∗
∼0 F2(g−1),2(n−g−1)
n−p−1 1−Λ∗ H
p≥1 g=2 p Λ∗ ∼0 Fp,n−p−1
√
∗ H
p≥1 g=3 n−p−2 1−
p
√ Λ
Λ∗
∼0 F2p,2(n−p−2)
∗ d
n  pg (large sample) −(n − 1 − p+g
2 ) log Λ −−→ χ2(g−1)p
H0

Table 2: Distribution of Λ∗ = det(W )/ det(B + W )

The s-way ANOVA model with s ∈ N may be handled in a similar manner. Suppose that there are two
factors under consideration. Proceeding by analogy, the two-way balanced MANOVA model reads
Xℓkr = µ + αℓ + βk + γℓk + εℓkr , ℓ = 1, · · · , g; k = 1, · · · , b; r = 1, · · · , n,
Pg Pb Pg Pb i.i.d.
where ℓ=1 αℓ = k=1 βk = ℓ=1 γℓk = k=1 γℓk = 0p , and εℓkr ∼ Np (0p , Σ). Note that

Xℓkr = X̄··· + (X̄ℓ·· − X̄··· ) + (X̄·k· − X̄··· ) + (X̄ℓk· − X̄ℓ·· − X̄·k· + X̄··· ) + (Xℓkr − X̄ℓk· ),
=µ̂ =α̂ℓ =β̂k =γ̂ℓk =ε̂ℓkr

where a dot is used to denote averaging over the indicated subscript. Accordingly, the partition of sums of
squares and cross products (SSP) is
∑
g
∑
b ∑
n ∑
g
(Xℓkr − X̄··· )(Xℓkr − X̄··· )′ = bn(X̄ℓ·· − X̄··· )(X̄ℓ·· − X̄··· )′
ℓ=1 k=1 r=1 ℓ=1

∑
b
+ gn(X̄·k· − X̄··· )(X̄·k· − X̄··· )′
k=1

∑
g
∑
b
+ n(X̄ℓk· − X̄ℓ·· − X̄·k· + X̄··· )(X̄ℓk· − X̄ℓ·· − X̄·k· + X̄··· )′
ℓ=1 k=1

∑
g
∑
b ∑
n
+ (Xℓkr − X̄ℓk· )(Xℓkr − X̄ℓk· )′ .
ℓ=1 k=1 r=1

14
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

The variations are summarized in table 3.

Source SSP d.o.f.

P
g
Factor 1 SSPfac1 = bn(X̄ℓ·· − X̄··· )(X̄ℓ·· − X̄··· )′ g−1
ℓ=1
Pb
Factor 2 SSPfac2 = gn(X̄·k· − X̄··· )(X̄·k· − X̄··· )′ b−1
k=1
P
g P
b
Interaction SSPint = n(X̄ℓk· − X̄ℓ·· − X̄·k· + X̄··· )(X̄ℓk· − X̄ℓ·· − X̄·k· + X̄··· )′ (g − 1)(b − 1)
ℓ=1 k=1
P
g P
b P
n
Residual SSPres = (Xℓkr − X̄ℓk· )(Xℓkr − X̄ℓk· )′ gb(n − 1)
ℓ=1 k=1 r=1
Pg P b P n
Total SSPT = (Xℓkr − X̄··· )(Xℓkr − X̄··· )′ gbn − 1
ℓ=1 k=1 r=1

Table 3: Two-way balanced MANOVA

Ordinarily, the test for interaction γℓk is conducted before the tests for main effects αℓ , βk . If interaction
effects exist, the factor effects will not have a clear interpretation. An interaction effect means that the

Figure 2: Insignificant interaction (upper half) vs. significant interaction (lower half)

effect of one factor depends on the other factor, and it’s indicated by the lines in our profile plot not running
parallel. Figure 3 shows 4(medicine) × 2(gender) = 8 observed mean BDI scores. In this case, the effect

Figure 3: Profile plot of a two-way ANOVA with an interaction effectxviii)

for medicine interacts with gender, i.e., medicine affects females differently from males. Roughly speaking,
the red line (Female) descends quite steeply from “None” to “Pharmaceutical” whereas the blue line (Male)
is much more horizontal. Depending on gender, the main effect of medicine obscures –rather than clarifies–
how medicine really affects the BDI scores, and thus should be ignored even if it’s statistically significant.
xviii) from [Link]

15
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

First, test H0 : γℓk = 0p , ∀ℓ = 1, · · · , g; k = 1, · · · , b using

det(SSPres )
Λ∗ =
det(SSPint + SSPres )
with h i
log Λ∗ −−→ χ2(g−1)(b−1)p .
p+1−(g−1)(b−1) d
− gb(n − 1) − 2 H0

If Λ∗ < c, reject H0 ; otherwise, continue.

Once the full model X = Xmain +Xint could be restricted to X = Xmain , its deserved orthogonal designxix)
would allow us to test main effects of Factor 1 and 2 separately.
• H01 : α1 = · · · = αg = 0p may be tested using

det(SSPres )
Λ∗1 =
det(SSPfac1 + SSPres )
with h i
log Λ∗1 −−→ χ2(g−1)p .
p+1−(g−1) d
− gb(n − 1) − 2 H0

• H02 : β1 = · · · = βb = 0p may be tested using

det(SSPres )
Λ∗2 =
det(SSPfac2 + SSPres )
with h i
log Λ∗2 −−→ χ2(b−1)p .
p+1−(b−1) d
− gb(n − 1) − 2 H0

§13 Linear Regression (2019/11/14)

To begin with, we discuss the univariate regression model

yi = β ′ xi + ε i , i = 1, · · · , n,

where the response yi is one-dimensional and β ∈ Rp consists of unknown parameters β1 , . . . , βp . When the
explanatory variables xi1 , . . . , xip are more than one (i.e., p > 1), the model is called multiple regression,
compared to simple linear regression yi = a + bzi + εi . These n equations are often stacked together and
written in matrix notation as
y = Xβ + ε,
where         
y1 x′1 x11 x12 ... x1p β1 ε1
 y2   x′2   x21 x22 ... x2p   β2   ε2 
         
y =  . , X= . = . .. ..  , β =  . , ε =  . .
 ..   ..   .. . .   ..   .. 
yn x′n xn1 xn2 ... xnp βp εn
Besides, the errors are assumed to satisfy that Eεi = 0, Var(εi ) = σ for some unknown σ 2 > 0, and
2

Cov(εi , εj ) = 0 for i 6= j. That is to say,

Eε = 0n & Var(ε) = σ 2 In .

The method of least squares yields

β̂ = arg min(y − Xβ)′ (y − Xβ) = (X ′ X)− X ′ y.

β∈Rp

Note that the generalized inverse (X ′ X)− of X ′ X ∈ Rp×p may not be well-defined. Nevertheless, the
(ordinary) least squares estimator of β refers unambiguously to

β̂OLS = (X ′ X)+ X ′ y = X + y,
xix) cf. [Link]

16
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

taking advantage of the unique Moore–Penrose inversexx) . Denote rank (X) = r ≤ p, then X = X∗ C for
some X∗ ∈ Rn×r and C ∈ Rr×p . Let β∗ = Cβ ∈ Col(C) = Rr . It follows that

Xβ = X∗ Cβ = X∗ β∗ .

Note that β̂∗ = (X∗′ X∗ )−1 X∗′ y is uniquely defined.

In many applications, we are interested in estimating some linear functions θ = c′ β of β, where c ∈ Rp .
Recall that ϑ is said to be estimable if and only if there exists an unbiased estimator of ϑ.
Theorem (estimability in linear models). The following facts suggest that c′ β is estimable if c ∈ Col(X ′ ),
which is also necessary for c′ β to be estimable under the assumption that ε is normally distributed.
(1) If c ∈ Col(X ′ ), then c′ β̂ is unique and unbiased for c′ β.
(2) If c′ β is estimable and ε ∼ Nn (0n , σ 2 In ), then c ∈ Col(X ′ ).
Proof . (1) If c = X ′ l for some l ∈ Rn , then c′ β̂ = l′ X(X ′ X)− X ′ y. It suffices to show
Lemma. The matrix PX = X(X ′ X)− X ′ is unique. Indeed, PX is the orthogonal projection onto the
column space Col(X) of the design matrix X ∈ Rn×p .
∈ Rn , write v = x + w for x ∈ Col(X) and w ∈ Col(X)⊥ , which exist and are
Proof of Lemma. For any vL
unique since R = Col(X) Col(X)⊥ . It follows from X ′ w = 0p that PX w = 0n and PX v = PX x. To
n

prove PX x = x, it is equivalent to show X(X ′ X)− X ′ X = X , or u′ X(X ′ X)− X ′ X = u′ X, ∀u ∈ Rn . Now

that X ′ X(X ′ X)− X ′ X = X ′ X by definition, it suffices to point out that u′ X = z ′ X ′ X for some z ∈ Rp .
What we need is exactly that Col(X ′ ) = Col(X ′ X) . Clearly Col(X ′ X) ⊂ Col(X ′ ), so the well-known
relation rank (X ′ X) = rank (X) = rank (X ′ ) completes the proof.

Therefore, c′ β̂ = l′ PX y is unique, and

E[c′ β̂] = l′ PX Ey = l′ PX Xβ = l′ Xβ = c′ β.

Note that c = X ′ l ∈ Col(X ′ ) =⇒ E[l′ y] = l′ Xβ = c′ β.

Remark. We say that c′ β is linearly estimable if, there exists some l ∈ Rn such that

E[l′ y] = c′ β, ∀β ∈ Rp .

(2) If there is an estimator T (y, X) unbiased for c′ β, then

Z  
′ 1 1 2
cβ= T (y, X) exp − 2 ky − Xβk dy.
Rn (2πσ 2 )n/2 2σ
Differentiation with respect to β yields
Z  
′ y − Xβ 1 2
c=X T (y, X) exp − 2 ky − Xβk dy,
Rn (2πσ 2 )n/2 σ 2 2σ
and thus c ∈ Col(X ′ ). 

In summary, Aβ with A ∈ Rq×p is estimable iff Col(A′ ) ⊂ Col(X ′ ), i.e., A = A∗ X for some A∗ ∈ Rq×n .
Particularly, β is estimable if and only if X has full column rank.
Sometimes denoted by H as well, the projection matrix PX = X(X ′ X)− X ′ is also named hat matrix
as it “puts a hat on y”, mapping the vector of response values to the vector of fitted values

ŷ = X β̂ = Hy.

The vector of residuals is

ε̂ = y − ŷ = (In − H)y.
Since HX = X, X ′ H = X ′ and
H ′ = H ′ H = (H ′ H)′ = H
(the check is left to the reader), we have X ′ ε̂ = 0p and ŷ ′ ε̂ = 0.
xx) cf. [Link]

17
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Figure 4: Projection in linear regression

We may view the design matrix from the SVD (see Lecture §2)

X
r
X = U D V′ = dj uj vj′
(n×p) (n×n)(n×p)(p×p)
j=1
 
diag(d1 , . . . , dr , 0, · · · , 0)
= (u1 , . . . , un ) (v1 , . . . , vp )′ .
0(n−p)×p

The row space of X is Col(X ′ ) = Col(X ′ X) = span (v1 , . . . , vr ), since

X
r
X ′X = d2j vj vj′ .
j=1

The column space of X is Col(X) = Col(XX ′ ) = span (u1 , . . . , ur ), since

X
r
XX ′ = d2j uj u′j .
j=1

In addition, Col(X)⊥ = span (ur+1 , . . . , un ). Moreover,

X
r
PX = X(X ′ X)− X ′ = uj u′j = Ur Ur′ ,
j=1

and thus
X
r X
n
PX y = hy, uj i uj , & (In − PX )y = (U U ′ − Ur Ur′ )y = hy, uj i uj .
j=1 j=r+1

The above results are extendable to infinite-dimensional separable Hilbert spaces, e.g., Fourier expansions.
In multivariate regression,
′
Yi = β Zi + Ui , i = 1, · · · , n,
(p×1) (r×p) (r×1) (p×1)

i.i.d.
where Ui ∼ (0p , Σ) for some unknown Σ = (σjk ) ∈ Rp×p . Denote
 ′  ′  
Y1 Z1 U1′
 ..   .. 
 
Y =  .  = (Y(1) , . . . , Y(p) ), Z =  .  , β = β(1) , . . ., β(p) , U =  ...  = (U(1) , . . . , U(p) ).
...
Yn′ Zn′ (r×1) (r×1)
Un′

18
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Henceforth,
Y = Z β + U ,
(n×p) (n×r)(r×p) (n×p)

where U is a data matrix from (0p , Σ). Rewrite

Y(j) = Z β(j) + U(j) , j = 1, · · · , p,
(n×r)
(n×1) (r×1) (n×1)

where Var(U(j) ) = σjj In and Cov(U(j) , U(k) ) = σjk In . The method of least squares yields β̂(j) = (Z ′ Z)− Z ′ Y(j) ,
so

β̂ = β̂(1) , . . . , β̂(p) = (Z ′ Z)− Z ′ Y
minimizes the total variation
X
p
β 7→ tr ((Y − Zβ)′ (Y − Zβ)) = (Y(j) − Zβ(j) )′ (Y(j) − Zβ(j) )
j=1

of all observations. The fitted values are

Ŷ = Z β̂ = Z(Z ′ Z)− Z ′ Y = PZ Y.
The residuals are
Û = Y − Ŷ = (In − PZ )Y.
It’s easily seen that
Z ′ Û = 0r×p , i.e., Z ⊥ Û(j) ;
and
Ŷ ′ Û = 0p×p , i.e., Ŷ(j) ⊥ Û(k)
for all 1 ≤ j, k ≤ p. Furthermore, we have the decomposition of sum of squares and cross products
Y ′Y = Ŷ ′ Ŷ + Û ′ Û .
(total) (regression) (residual)

For instance, Ŷ = n1 1n 1′n Y = Ȳ 1n if Z = 1n , and all formulae above may come into play.
Assume that the design matrix Z ∈ Rn×r has full column rank, in which case β = (β(1) , . . . , β(p) ) ∈ Rr×p
is estimable. We now investigate properties of β̂ = (Z ′ Z)−1 Z ′ Y , the least squares estimator. First,
Eβ̂ = (Z ′ Z)−1 Z ′ EY = (Z ′ Z)−1 Z ′ Zβ = β,
which demonstrates that β̂ is unbiased for β, and thus EÛ = 0n×p . Second,
Cov(β̂(j) , β̂(k) ) = (Z ′ Z)−1 Z ′ Cov(U(j) , U(k) )Z(Z ′ Z)−1 = σjk (Z ′ Z)−1
for 1 ≤ j, k ≤ p. Third,
Cov(β̂(j) , Û(k) ) = (Z ′ Z)−1 Z ′ Cov(U(j) , U(k) )(In − Z(Z ′ Z)−1 Z ′ ) = σjk (Z ′ Z)−1 (Z ′ − Z ′ ) = 0r×n .
When it comes to error terms, note that
 
′ ′ ′
E[Û(j) Û(k) ] = E[U(j) (In − PZ )U(k) ] = tr (In − PZ )E[U(k) U(j) ] = σkj tr (In − PZ ) = (n − rank (Z))σjk .
1
Thus, n−rank(Z) E[Û ′ Û ] = Σ, and Σ̂ = n−rank(Z)
1
Û ′ Û is called the within-sample MSExxi) of the predictor.
Prediction of a new observation is more uncertain than estimating the expected value. See also exercise 7
in homework 1. Consider
Y0 = β ′ z0 + U0
(new response) (expected value) (new error)

at z0 ∈ R . Suppose that Z ∈ R
r n×r
has full column rank, then z0′ β̂ is an unbiased predictor of z0′ β. The
forecast error is
Y0 − β̂ ′ z0 = (β − β̂)′ z0 + U0 ,
whose covariance matrix consists of
Cov(z0′ (β(j) − β̂(j) ), z0′ (β(k) − β̂(k) )) + Cov(U0,j , U0,k ) = σjk [1 + z0′ (Z ′ Z)−1 z0 ], 1 ≤ j, k ≤ p.
In short, Var(Y0 − β̂ ′ z0 ) = [1 + z0′ (Z ′ Z)−1 z0 ]Σ, where z0′ (Z ′ Z)−1 z0 ≥ 0.
xxi) the definition of the mean squared error differs according to whether one is describing a predictor or an estimator (see

Lecture §3 for the latter).

19
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§14 Multivariate Linear Regression Cont’d (2019/11/19)

From now on, we assume, WLOG, that Z, the design matrix in the linear model

Y = Z β + U ,
(n×p) (n×r)(r×p) (n×p)

has full column rank so that β is estimable (see Lecture §13); otherwise, we may pick out some maximal
linearly independent columns of Z, resulting in dimensionality reduction of covariates. Denote by

β̂ = (Z ′ Z)−1 Z ′ Y = β + (Z ′ Z)−1 Z ′ U

the least squares estimator (LSE) of β. Note that the rows of the error term U are uncorrelated zero-mean
random vectors, whose covariance matrices are identically Σ ∈ Rp×p , an unknown nuisance parameter.
The Gauss–Markov theorem claims that β̂ is the best linear unbiased estimator (BLUE) of β, in
the sense that for any c′ β (c ∈ Rr ), of the estimators l′ Y (l ∈ Rn ) that are unbiased, c′ β̂ has minimum
covariance matrix in Loewner order.
Proof. Since c′ β = E[l′ Y ] = l′ Zβ for all β, we have c = Z ′ l. To show that vec(c′ β̂) = vec(l′ PZ Y ) =
(Ip ⊗ (l′ PZ )) vec(Y ) has smaller covariance matrix than vec(l′ Y ) = (Ip ⊗ l′ ) vec(Y ), direct calculations yield

Var(vec(c′ β̂)) = (Ip ⊗ (l′ PZ ))(Σ ⊗ In )(Ip ⊗ (l′ PZ ))′ = (l′ PZ l)Σ

and Var(vec(l′ Y )) = (l′ l)Σ. It suffices to observe that l′ (In − PZ )l ≥ 0.

Hereafter, assume further that the rows of U are i.i.d. random vectors drawn from Np (0p , Σ), i.e.,
vec(U ) ∼ Nnp (0np , Σ ⊗ In ). In such a Gaussian framework, the MLE of β is identical to its LSE. Actually,
the likelihood function
( !)
1 Xn
−np/2 −n/2 −1 ′ ′ ′
Ln (β, Σ|Y1 , . . . , Yn ) = (2π) det(Σ) exp − tr Σ (Yi − β Zi )(Yi − β Zi )
2 i=1

attains its maximum at

(β̂ MLE , Σ̂MLE ) = (β̂, n1 Û ′ Û ),

where Û = Y − Ŷ = (In − PZ )Y has Ûi = Yi − β̂ ′ Zi as its rows. To see this, notice that the SSP of

Yi − β ′ Zi = Ûi + (β̂ − β)′ Zi

is
X
n X
n X
n
(Yi − β ′ Zi )(Yi − β ′ Zi )′ = Ûi Ûi′ + (β̂ − β)′ Zi Zi′ (β̂ − β) = Û ′ Û + (β̂ − β)′ Z ′ Z(β̂ − β),
i=1 i=1 i=1

since the cross terms add up to

X
n
Ûi Zi′ (β̂ − β) = Û ′ Z(β̂ − β) = 0p×r (β̂ − β) = 0p×p .
i=1

Combined with exercise 2 in homework 2, the decomposition of SSP makes the expressions of β̂ MLE and
Σ̂MLE clear. Since Σ̂MLE is biased, the unbiased Σ̂MoM = n−r
1
Û ′ Û (method of moments estimator) prevails
in linear regression. Note that

nΣ̂MLE = Y ′ (In − Z(Z ′ Z)−1 Z ′ )Y = (Y − Zβ)′ (In − PZ )(Y − Zβ) = U ′ (In − PZ )U ∼ Wp (Σ, n − r),

not depending on β.  
β1
After estimation, let’s get down to inferences, say, the likelihood-ratio test (LRT). Write β = and
β2
Z = (Z1 , Z2 ) such that β1 ∈ Rq×p , β2 ∈ R(r−q)×p , Z1 ∈ Rn×q , and Z2 ∈ Rn×(r−q) for some q < r. Then,

Y = Z1 β1 + Z2 β2 + U.

20
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

We are interested in H0 : β2 = 0(r−q)×p . Since Y = Z1 β1 +U under H0 , similar arguments give constrained

MLEs
β̂1 = (Z′1 Z1 )−1 Z′1 Y
and
Σ̂MLE
1 = 1
n (Y − Z1 β̂1 )′ (Y − Z1 β̂1 ) = 1 ′
n Y (In − PZ1 )Y.
The LRT statistic is
!n/2
Ln (β̂1 , Σ̂MLE
1 ) (2π)−np/2 det(Σ̂MLE
1 )−n/2 exp(− np
2 ) det(Σ̂MLE )
Λ= = = .
Ln (β̂, Σ̂MLE ) (2π)−np/2 det(Σ̂MLE )−n/2 exp(− np
2 ) det(Σ̂MLE
1 )
Also, Wilks’ Lambda
det(Σ̂MLE )
Λ∗ = Λ2/n =
det(Σ̂MLE
1 )
is commonly used. Here, the SSP decomposition is
(Y − Z1 β̂1 )′ (Y − Z1 β̂1 ) = (Y − Zβ̂)′ (Y − Zβ̂) + (Zβ̂ − Z1 β̂1 )′ (Zβ̂ − Z1 β̂1 ),
=nΣ̂MLE
1 =nΣ̂MLE =n(Σ̂MLE
1 −Σ̂MLE )

where the cross terms add up to

(Y − Zβ̂)′ (Zβ̂ − Z1 β̂1 ) = Y′ (In − PZ )(PZ − PZ1 )Y = 0p×p ,
since Col(Z1 ) ⊂ Col(Z). From the perspective of Cochran’s theorem (see Lecture §9),
In = PZ1 + (PZ − PZ1 ) + (In − PZ )
decomposes into mutually orthogonal projections onto Col(Z1 ), Col(Z) ∩ Col(Z1 )⊥ , and Col(Z)⊥ , respec-
tively. Therefore,
H H
nΣ̂MLE
1 ∼0 Wp (Σ, n − q), & n(Σ̂MLE
1 − Σ̂MLE ) ∼0 Wp (Σ, r − q).
Moreover, the large sample approximation is
!
  det(Σ̂MLE ) ◦
−2 log Λ ≈ − n − r − 1 − 21 (p − r + q + 1) log ∼ χ2(r−q)p .
det(Σ̂MLE
1 )
Generally, to test H0 : C β= A , we can use the fact that
(r−q)×r (r−q)×p

− Σ̂MLE ) = (C β̂ − A)′ [C(Z ′ Z)−1 C ′ ]−1 (C β̂ − A) ∼0 Wp (Σ, r − q),

H
n(Σ̂MLE
0

H
and nΣ̂MLE
0 ∼0 Wp (Σ, n − q).
Recall that
M mp
T 2 = d′ ( )−1 d ∼ T 2 (p, m) = Fp,m−p+1
m m−p+1
when d ∼ Np (0p , Σ) ⊥⊥ M ∼ Wp (Σ, m). For any z0 ∈ Rr , we have

β̂ ′ z0 ∼ Np (β ′ z0 , z0′ (Z ′ Z)−1 z0 Σ) ⊥⊥ nΣ̂MLE ∼ Wp (Σ, n − r),

and thus
!′ !−1 !
β̂ ′ z − β ′ z0 nΣ̂MLE β̂ ′ z − β ′ z0 (n − r)p
p 0 p 0 ∼ T 2 (p, n − r) = Fp,n−r−p+1 .
z0′ (Z ′ Z)−1 z0 n−r z0′ (Z ′ Z)−1 z0 n−r−p+1

In the course of predicting/forecasting Y0 = β ′ z0 + U0 , we have

Y0 − β̂ ′ z0 = (β − β̂)′ z0 + U0 ∼ Np (0p , [1 + z0′ (Z ′ Z)−1 z0 ]Σ) ⊥⊥ nΣ̂MLE ∼ Wp (Σ, n − r),

and thus
!−1
nΣ̂MLE (n − r)p
(Y0 − β̂ ′ z0 )′ (Y0 − β̂ ′ z0 ) ∼ [1 + z0′ (Z ′ Z)−1 z0 ] Fp,n−r−p+1 .
n−r n−r−p+1

21
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§15 Model Diagnostics and Selection (2019/11/26)

Recall the univariate linear model (see Lecture §13)

y = X β + ε .
(n×1) (n×p)(p×1) (n×1)

The vector of residuals is

ε̂ = y − ŷ = (In − X(X ′ X)− X ′ )y.
Denote
H = X(X ′ X)− X ′ = (hij )1≤i,j≤n .
Since H = H ′ = H 2 , we have
X
n X
n
hij = hik hkj = hik hjk .
k=1 k=1

The diagonal elements of H are called the leverages (i.e., influential points). Note that

X
n
hii = tr (H) = rank (H) = rank (X) ,
i=1

and X
hii (1 − hii ) = h2ij ≥ 0 =⇒ hii ∈ [0, 1].
j̸=i

If hii % 1, then X
ŷi = hii yi + hij yj
j̸=i

reflects more contribution from yi . However, for prediction, we hope that yj (j 6= i) contribute to ŷi as
much as possible, while yi contributes to ŷi as even as possible.
It is imperative to examine the adequacy of the model before making decisions. Since

Var(ε̂) = (In − H) Var(y)(In − H)′ = (In − H)(σ 2 In )(In − H) = σ 2 (In − H),

we often prefer studentized residuals

ε̂i
ε̂∗i = p , i = 1, · · · , n
s2 (1 − hii )
1
Pn 2
using the mean squared residual s2 = n−p 2
i=1 ε̂i as an estimate of σ . Also, we may studentize ε̂i using
the delete-one estimated variance s2(−i) , computed from the regression with the ith observation dropped.
The following graphical diagnostics are based on these studentized residuals:
• Q-Q plots and histograms for normality.
• Scatter plots of ε̂∗i vs. i (indexes) for temporal/spatial correlation.
• Scatter plots of ε̂∗i vs. ŷi for homoscedasticity.
• Scatter plots of ε̂∗i vs. Xi for flexible fit (e.g., need for more terms in the model).
Note that there exist some subjective judgements. A good sense comes from practices!
In spite of most realistic situations that the true model is not in our candidate class, we shall try our
best to make our model useful, i.e., to approximate the true model, which often involves nonparametric
(infinite-dimensional) methods. If, coincidentally, the true model belongs to the candidate class, then our
goal is simply to find it out, in which case we often utilize parametric (finite-dimensional) methods and pay
attention to consistency. These are the so called model selection. For example, we often use penalties
or shrinkage in high-dimensional statistics (p  n), where models rest on a few tuning parameters. For
the purpose of selection, some criteria are needed. Note that models, offering ways of inference for the
population feature, ultimately serve the prediction.

22
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Suppose the true model is

y = µ + ε, µ ∈ Rn , ε ∼ (0n , σ 2 In ).
The candidate class consists of stepwise regression models

y = Xp βp + ε, p = 1, · · · , m

where Xp ∈ Rn×p form a nested sequence X1 ⊂ X2 ⊂ · · · ⊂ Xm , e.g., splines with more and more knots,
or Fourier expansions of higher and higher order. The method of least squares yields β̂p = (X′p Xp )−1 X′p y,
and thereby
µ̂p = Xp β̂p = Xp (X′p Xp )−1 X′p y = Pp y.
Clearly, the bias
bp = µ − Eµ̂p = (In − Pp )µ
decreases in p. Particularly, if µ = Xp∗ βp∗ for some p∗ , i.e., the true model belongs to the candidate class,
then bp = 0n for all p ≥ p∗ . Besides, the mean squared error of µ̂p ,
2 2 2 2 2
dp := E kµ̂p − µk = E kPp ε − (In − Pp )µk = E kPp εk + k(In − Pp )µk = pσ 2 + kbp k ,
2
decomposes into its variance pσ 2 and its squared bias kbp k . Generally, the optimal p can be estimated by

p̂ = arg min dˆp ,

1≤p≤m

where the estimate dˆp of dp should be specified. Since the sum of squared errors for the candidate model
with p covariates,
2 2
SSEp = ky − Xp β̂p k2 = k(In − Pp )(µ + ε)k = kbp + (In − Pp )εk ,

satisfies that
2 2 2
ESSEp = E k(In − Pp )εk + kbp k = (n − p)σ 2 + kbp k ,
we can define
dˆp = SSEp − (n − 2p)σ̂ 2 ,
where σ̂ 2 is unbiased so that Edˆp = dp . A similar criterion is Akaike’s final prediction error (FPE).
Imagine a vector of new data y∗ = µ + ε∗ . The expectation of the sum of squared prediction errors is

dp := EPEp = E ky∗ − µ̂p k = nσ 2 + E kµ̂p − µk = (n + p)σ 2 + kbp k .

2 2 2

Hence, we may take

dˆp = FPEp = SSEp + 2pσ̂ 2 ,
and then we are faced with the trade-off between SSEp and 2pσ 2 . So, how to estimate σ 2 ?
1) If we choose σ̂ 2 = MSEp = SSEp /(n − p), then

dˆp = SSEp + 2pMSEp = n+p

n−p SSEp .

For underfit models, the MSEp can be dramatically larger than σ 2 . Thus, a large p̂ will be encouraged.
2) If we choose σ̂ 2 = MSEm typically, using the largest model, then

dˆp = SSEp + 2p SSE

n−m ,
m

which, however, cannot capture consistent models.

One of the most popular criterion,
Cp = 1
n (SSEp + 2pσ̂ 2 ),
is called Mallow’s Cp (with σ̂ 2 = MSEm ), sometimes defined as

Cp′ = SSEp /σ̂ 2 − (n − 2p) = nCp /σ̂ 2 − n.

Other criteria to be minimized among 1 ≤ p ≤ m include

23
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

• Akaike information criterion

AIC = −2 log(L̂n ) + 2p,
where L̂n is the maximum value of the likelihood function of the model with p covariates. In the case
of (Gaussian) linear regression, ignoring the constant terms allows us to conveniently take

AIC = n log(SSEp /n) + 2p,

which is asymptotically equivalent to Mallow’s Cp .

• Bayesian information criterion (developed by Schwarz who gave a Bayesian argument)

BIC = −2 log(L̂n ) + log(n)p,

where L̂n is the maximum value of the likelihood function of the model with p covariates. In the case
of (Gaussian) linear regression, ignoring the constant terms allows us to conveniently take

BIC = n log(SSEp /n) + log(n)p.

When fitting models, it is possible to increase the likelihood by adding parameters, but doing so may result
in overfitting. Both BIC and AIC attempt to resolve this problem by introducing a penalty term for the
number of parameters in the model; the penalty term is larger in BIC than in AIC.

§16 Model Selection Cont’d & Principal Components (2019/11/28)

Let p̂M be the model chosen by Mallow’s Cp , which is asymptotically equivalent to AIC, and let p̂S be the
model chosen by Schwarz’s BIC. That is,

p̂M = arg min Cp , and p̂S = arg min BICp .

p∈{1,··· ,m} p∈{1,··· ,m}

Their optimalities are as follows. Assume that m → ∞ as n → ∞.

• If µ 6= Xp βp for p = 1, · · · , m, i.e., the true model is not in the candidate class, then

kµ − µ̂p̂M k P dp̂M P
−
→ 1, and −
→ 1.
inf1≤p≤m kµ − µ̂p k inf1≤p≤m dp

• If µ = Xp∗ βp∗ for some p∗ , i.e., the true model belongs to the candidate class, then

P{p̂S = p∗ } → 1, while P{p̂M > p̂S } > 0.

A widely used method for assessing the performance of prediction models is cross-validation (CV).
In K-fold cross-validation, the original sample (xi , yi )1≤i≤n is randomly partitioned into K equal-sized
subsamples. Let κ : {1, · · · , n} → {1, · · · , K} be an indexing function that indicates the partition to
which observations are allocated by the randomization, respectively. Denote by fˆ−k (·) the fitted function,
computed with the k th part of the data removed. Then the estimate of prediction error is

1 X 1
K X h i2 1 Xh
n i2
CV = yi − fˆ−k (xi ) = yi − fˆ−κ(i) (xi ) .
K n/K n i=1
k=1 κ(i)=k

In practice, one typically performs K-fold CV using K = 5 or K = 10, which is computationally efficient.
The case K = n is known as leave-one-out cross-validation, where

1 Xh i2
n
LOOCV = yi − fˆ(−i) (xi ) .
n i=1

Consider the linear regression model with p regressors, where f (x) = x′ β and y −f (x) = ε ∼ (0, σ 2 ). Clearly
we have fˆ(x) = x′ β̂ and fˆ(−i) (x) = x′ β̂(−i) by the method of least squares. Forecasting a new response y ∗
 2
at x∗ gives prediction error PE = y ∗ − fˆ(x∗ ) . Then it holds that ELOOCV = EPE + O(p/n2 ).

24
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Generally, LOOCV can be easily calculated for ridge-type regression, where the hat matrix takes the
form
H = X(X ′ X + λJ)−1 X = (hij )1≤i,j≤n .
The goal is
" #2
1 X yi − fˆ(xi )
n
LOOCV = ,
n i=1 1 − hii
so it suffices to show that
yi − fˆ(xi )
yi − fˆ(−i) (xi ) = .
1 − hii
 ′
Proof. Let ỹ = (ỹj )1≤j≤n = y1 , . . . , yi−1 , fˆ(−i) (xi ), yi+1 , . . . , yn . Then the minimizer fˆ(−i) (·) of the
P 2 Pn
2
leave-one-out penalized sum of squares [yj − f (xj )] + λJ (f ) also minimizes [ỹj − f (xj )] + λJ (f ).
j̸=i j=1
  P n Pn
It follows that fˆ(−i) (xj ) = H ỹ. Just combine fˆ(−i) (xi ) = hij ỹj with fˆ(xi ) = hij yj .
1≤j≤n j=1 j=1

It’s immediate that LOOCV is not stable if some hii % 1. To alleviate the tendency to undersmooth,
a feasible approximation is generalized cross-validation, given as
" #2
1X
n
yi − fˆ(xi )
GCV = .
n i=1 1 − tr (H) /n

Using the fact that (1 − p/n)−2 ≈ (n + p)/(n − p) when tr (H) = p  n, we have

1 p −2 X h i2
n
1 n+p 1
GCV = 1− yi − fˆ(xi ) ≈ SSE = FPE.
n n i=1
nn−p n
   
AIC MLE 2
The above discussions can be extended to multivariate cases, e.g., = n log(det Σ̂d )+ pd
BIC log(n)
for Gaussian linear regression of p-dimensional responses with d covariates.
For computational and statistical reasons, best subset selection may suffer from an enormous search
space. Thus, stepwise methods are attractive alternatives. Hybrid versions of forward and backward
stepwise selection attempts to more closely mimic best subset selection while retaining the computational
advantages of forward and backward stepwise selection. Rather than subset selection methods, we may fit
a model containing all predictors using a technique that regularizes the coefficient estimates, i.e., we may
introduce a shrinkage penalty and select the tuning parameter.
In multivariate statistics, data reduction plays a central role, which transforms the predictors and
reproduces the variability. Principal components analysis (PCA) is commonly used for deriving a low-
dimensional set of features from a large set of variables. As for intuition, (xi , 0)1≤i≤n should be extracted
from its perturbed realization (xi , ui )1≤i≤n , where |ui |  1. For any data cloud, we should find some
directions onto which the projection of the data shows the most information/variation.
From the point of view of population, let X be a random p-vector and suppose that EX = 0p for
simplicity. We are to find some unit vector ℓ1 ∈ Rp such that Y1 = ℓ′1 X captures most of the variability.
Indeed,
ℓ1 = arg max Var(ℓ′ X) = arg max ℓ′ Σℓ
ℓ∈Rp :∥ℓ∥=1 ℓ∈Rp :∥ℓ∥=1

turns out to be the first eigenvector of Var(X) = Σ. Sequentially, the j th principal component of X
is Yj = ℓ′j X, where ℓj = arg maxℓ∈Rp ℓ′ Σℓ s.t. kℓk = 1 and ℓ′ ℓk = 0 for 1 ≤ k < j. Using the eigen-
decomposition
X
p
Σ = (e1 , . . . , ep ) diag(λ1 , . . . , λp ) (e1 , . . . , ep )′ = λj ej e′j ,
j=1

where λ1 > · · · > λp xxii)

, we have xxiii)
Yj = e′j X and thus Var((Yj )1≤j≤p ) = diag(λ1 , . . . , λp ).
xxii) the case when λj = λk for some j ̸= k (multiplicity) is usually not considered in statistical analysis, since observations are
assumed to be drawn randomly.
xxiii) cf. [Link]

25
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

§17 Principal Component Analysis Cont’d (2019/12/3)

For unknown Σ = Var(X), we can estimate it by the sample covariance matrix S or the MLE Σ̂. If n  p,
then the sample principal components (PCs) are consistent in the sense that

{êj , λ̂j } →P {ej , λj }.

One may consider the random effect model (linear mixed-effects model)
Yi = Xi′ β + Zi′ γi + εi ,
mean (fixed) subject-specific (random)

where the PCs seem inclined to approximate these fixed effects.

Back to the population view: let
Q = (e1 , . . . , ep ),
then
Y = (Y1 , . . . , Yp )′ = Q′ X,
and
X
p X
p
X = (X1 , . . . , Xp )′ = QY = Yj ej = hX, ej i ej .
j=1 j=1

Generally, for a (non-centered) random element X taking values in some separable Hilbert space with
EX = µ, we have the so called Karhunen-Loève expansionxxiv)
∞
X
X =µ+ hX − µ, ej i ej .
j=1

Note that the total variance is invariant in that

X
p X
p X
p
Var(Yj ) = λj = tr (Σ) = Var(Xj ).
j=1 j=1 j=1
Pp
The unexplained variance of the p-dimensional X = µ + Yj ej , with respect to the p1 -dimensional
Pp1 j=1
X∗ = µ + j=1 Yj ej , is
2
X
p X
p
EkX − X∗ k = E2
Yj ej = λj .
j=p1 +1 j=p1 +1

The fraction of variation explained by the first p1 PCs is then defined to be

, p
X
p1 X
FVE = λj λj ,
j=1 j=1

which reflects the model fidelity. However, p1 characterizes the model complexity. We again have a trade-off
between two goals, as we already encountered with MSEs. The optimal p∗1 is usually chosen with reference
to a plot of cumulative FVE against p1 — there are no universal rules, only rules of thumb. In view of this,
the distribution of λ’s is essential, whose properties have been intensively studied in random matrix theory.
As a result, one often see that p∗1 ∼ nα with a small α ∈ (0, 1).
PCA can also be done using the correlation matrix instead of the covariance matrix. Note that the
correlation matrix is exactly the covariance matrix of the scaled/standardized variables. Results of scaled
and unscaled PCA can differ a lot, since scaling makes variables equally important. Often (but not always),
scaling (i.e., using the correlation matrix) is preferred.
Anyway, PCA is a technique for dimensionality reduction that pursues a simplified structure of the
initial data, and belongs to unsupervised learning. Sometimes the PCs are used in supervised learning, e.g.,
regression/classification, where the explanation and selection will become much more clear.
The interpretation can be difficult. When all measurements are positively correlated, the first PC is
often some kind of average of the measurements (e.g., market general index). Then the other PCs give
important information about the remaining pattern (e.g., industrial sector index).
xxiv) cf. [Link]

26
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

The geometrical interpretation of PCs based on data is the minimization of perpendicular distances, as
shown in Figure 5, not to be confused with vertical distances in OLS, where responses and predictors are

Figure 5: PCA (left-hand side) vs. OLS (right-hand side)xxv)

unequalxxvi) . Formally, let x1 , . . . , xn ∈ Rp be centered observations, and Πq be the subspace spanned by

the first q vectors of an orthonormal basis f1 , . . . , fq , fq+1 , . . . , fp . It follows that
X
q
projΠq xi = hxi , fj i fj .
j=1

If
2 X
p
2
d2 (xi , Πq ) = xi − projΠq xi = hxi , fj i ,
j=q+1

then
X
n
span (e1 , . . . , eq ) = arg min d2 (xi , Πq ).
Πq ⊂Rp :dim(Πq )=q i=1
Pn
Proof. Denote by S = 1
n−1 i=1 xi x′i the sample covariance matrix. Direct calculations suggest that

X
n   X
p
d (xi , Πq ) = (n − 1) tr S
2
Π⊥
= (n − 1) fj′ Sfj
q
i=1 j=q+1

attains its minimum if and only if fq+1 , . . . , fp are the last p − q eigenvectors of S.

§18 Canonical Correlation & Factor Analysis (2019/12/10)

   
X ΣXX ΣXY
Given a random p-vector X and a random q-vector Y with Var = , the first pair
Y ΣY X ΣY Y
′ ′
of canonical variables consists of U1 = e1 X and V1 = f1 Y such that Corr(U1 , V1 ) is maximized among
−1/2 −1/2
e1 ∈ Rp and f1 ∈ Rq . Actually, e1 = ΣXX ẽ1 and f1 = ΣY Y f˜1 , where ẽ1 and f˜1 are the first left and
−1/2 −1/2
right singular vectors of the canonical correlation matrix C = ΣXX ΣXY ΣY Y . To see this, note that
1/2 1/2
if c = ΣXX a for a ∈ Rp and d = ΣY Y b for b ∈ Rq , then
a′ ΣXY b c′ Cd
Corr(a′ X, b′ Y ) = √ √ = .
a′ ΣXX a b′ ΣY Y b kck · kdk
Recall that the left and right singular vectors of C are eigenvectors of CC ′ and C ′ C, respectively. Se-
quentially, the k th pair of canonical variables is Uk = e′k X and Vk = fk′ Y that maximizes Corr(Uk , Vk )
among those linear combinations uncorrelated with the preceding k − 1 canonical pair(s). It turns out that
−1/2 −1/2
Uk = ẽ′k ΣXX X and Vk = f˜k′ ΣY Y Y , where ẽk and f˜k are the k th left and right singular vectors of C.
xxv) from JD Long’s blog on 16 Sep 2010. [Link]
xxvi) cf. [Link]

27
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Factor analysis is widely used in social sciences for qualitative modeling, where observable/manifest
variables X = (X1 , . . . , Xp )′ are explained in terms of a potentially smaller number of underlying latent
variables F = (F1 , . . . , Fm )′ that are often impossible to measure directly. Such unobservable variables of
interest are called common factors, because they influence the observed variables broadly. In psychology,
intelligence is a prime example, which researchers examine indirectly by measuring variables that are
believed to be its indicators, e.g., test scores in literature, language and mathematics. The essential purpose
of factor analysis is to exhibit the relationship between the latent and the observed variables.
Let µ = EX. The classic orthogonal factor model is

X
m
X −µ= L F + ε , or Xi − µi = ℓij Fj + εi , 1≤i≤p
(p×1) (p×m)(m×1) (p×1)
j=1

where ℓij is called the loading of the ith variable on the j th factor, meaning the weight that the j th factor
carries in determining the ith variable; and εi is the error (sometimes called the specific factor) associated
only with the ith variable. Note that the loading matrix L = (ℓij ) does not vary across observations. Also,
the common and specific factors are assumed to satisfy that
• EF = 0m , Var(F ) = Im to ensure that the factors are uncorrelated;
• Eε = 0p , Var(ε) = Ψ = diag(ψ1 , . . . , ψp ) so that Cov(Xi , Xi′ |F ) = 0 for i 6= i′ ; and
• Cov(F, ε) = 0m×p .

Note that the assumptions still hold if we set F ← Q′ F and L ← LQ for an orthogonal matrix Q ∈ Rm×m ,
and thereby the model is identifiable only up to a factor rotation. From the assumptions, it follows the
covariance structure Var(X) = LL′ + Ψ, or

X
m X
m
Var(Xi ) = ℓ2ij + ψi , & Cov(Xi , X ) =
i′ ℓij ℓi′ j , ∀i 6= i′ .
j=1 j=1

Pm
Here h2i = j=1 ℓ2ij is called the ith communality and ψi is often called the specific/unique variance.
It can also be seen that
Cov(X, F ) = Cov(LF, F ) + Cov(ε, F ) = L.
There is a strong connection between PCA and factor analysis (FA), but the two are not identical. Both
methods are mostly used in exploratory data analysis for dimension reduction. However, PCA aims at
explaining no covariances but the total variance, whereas FA is more concerned about accounting for the
covariances or correlations. For comparison, write the PCA as

X − µ = Q1 Y (1) + Q2 Y (2) = Q1 Y (1) + u,

(1)
where (Q1 ; Q2 ) = (e1 , . . . , em ; em+1 , . . . , ep ) = (qij ) consists of eigenvectors of Var(X), and Y P and Y (2) are
the first m and the last p−m principal component(s) of X, respectively. Note that Var(u) = j=m+1 λj ej e′j
p
Pm
is usually not diagonal. It’s difficult to interpret the ith variable Xi = j=1 qij Yj + ui . And we prefer the
Pm
entire extraction X = j=1 Yj ej + u. Indeed, PCA transforms observed variables to obtain principle com-
ponents, while FA transforms factors to obtain observed variables. It’s sensible to say how Xi ’s collectively
form/interpret Yj = e′j X, rather than how Yj ’s explain Xi .
The objective of factor analysis is to find the factor loadings L and the specific variances Ψ. In addition,
we need to select the number of factors m. The first step is to estimate Var(X) by calculating the sample
covariance matrix S, or the sample correlation matrix R, since the factor model is scale invariant in that
scaled (Xi /σi )1≤i≤p admits the same latent factors. The estimates L̂ of L and Ψ̂ of Ψ are deduced from
the covariance structure
S = L̂L̂′ + Ψ̂.
d.o.f.: p(p+1)/2 pm p

Whether the model has a solution or not is determined by the degree of freedom s = 21 p(p + 1) − (pm + p).

• If s < 0, then the model is meaningless: the number of factors m is “too large” relative to the number
of original parameters p and there are infinitely many solutions.

28
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

• If s = 0, then the solution is unique, but not necessarily proper — the specific variances Ψ are required
to be non-negative.
• In practice we usually have s > 0: an exact solution may not exist and approximate solutions are
used. This case is the most interesting for data/dimension reduction with interpretation.
There are two popular methods of estimation, the principal component (and the related principal factor)
method and the maximum likelihood method.
Principal component method leads to an iterative solution as follows:
1. Calculate eigenpairs (λi , ei ) of the sample covariance matrix S;
√ √  
2. Initialize L̂ = ( λ1 e1 , . . . , λm em ) and Ψ̂ = diag S − L̂L̂′ ;

3. Update L̂ using the eigen-decomposition of S − Ψ̂;

 
4. Update Ψ̂ = diag S − L̂L̂′ ;

5. Repeat steps 3. and 4. until convergence.

§19 Factor Analysis Cont’d & Classification (2019/12/12)

As an alternative for estimation of the factor model, maximum likelihood method assumes normality.
Recall that
n

log f (X; µ, Σ) = − log(det(2πΣ)) + tr Σ−1 Sn + (X̄ − µ)′ Σ−1 (X̄ − µ) ,
2
Pn
where Sn = n1 i=1 (Xi − X̄)(Xi − X̄)′ = n−1 ′
n S. Replacing µ with µ̂ = X̄ and substituting Σ = LL + Ψ,
the estimates L̂ and Ψ̂ are obtained by numerical minimization of


log(det(LL′ + Ψ)) + tr Sn (LL′ + Ψ)−1 .
It is desirable to make L̂ well defined by imposing the computationally convenient uniqueness condition
(check!) that L′ Ψ−1 L is diagonal. MLE method has an advantage: one can test whether the number of
factors m is adequate, i.e., Σ has the form LL′ + Ψ or not. The likelihood-ratio statistic is −2 log(Λ) = nF ,
where !
det(L̂L̂′ + Ψ̂)  
F = log + tr Sn (L̂L̂′ + Ψ̂)−1 − p.
det(Sn )
Bartlett has suggested that

d
n − 1 − 61 (2p + 5) − 23 m F −−→ χ2ν ,
H0
where the degree of freedom
   
ν = 12 p(p + 1) − pm + p − 12 m(m − 1) = 1
2 (p − m)2 − (p + m)
√
must be positive. In practice, we often start with the largest possible m < 21 (2p + 1 − 8p + 1). The
backward tests proceed as follows: if the adequacy of the m-factor model is not rejected, continue with
m − 1 factors, until rejection. Also, forward testing procedure works.

Figure 6: Orthogonal rotation Figure 7: Oblique rotation

As is shown, an orthogonal transformation of factors, referred to as factor rotation, preserves the

ability to reproduce the covariance. We can utilize the non-uniqueness property to deliberately change the
loadings in order to improve their interpretability. Figures 6 & 7 are visualsxxvii) of what happens during
xxvii) from [Link]

29
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

a rotation of two factors. Note that oblique rotations allow for correlated factors. The loadings can be
interpreted easily if a “simple clusters structure” is achieved — each variable is highly loaded on a single
factor and little on others. Given the estimated factor loadings L̂ = (ℓ̂ij ), Kaiser proposed an analytical
measure of simple structure known as the varimax criterion. The idea of the varimax rotation method is to
find the rotation that maximizes the variances of the squared loadings within each column since each factor
should have a few large and many negligible loadings. More precisely, the varimax criterion of L = (ℓij ) is
defined to be  !2 
Xm Xp Xp
VC(L) = 1 ℓe4ij −
1
ℓeij
2 ,
j=1
p i=1
p i=1
Pm
where ℓeij = ℓij /ĥi using the square root of the estimated ith communality ĥ2i = j=1 ℓ̂2ij . The varimax
optimal rotation is chosen so that VC(L̂Q) is maximized among orthogonal matrices Q ∈ Rm×m .
The aim of discriminant analysis is to classify objects into groups that are known a priori. Assume
that data come from g distinct classes/populations D1 , . . . , Dg . Given an observed random vector x, which
population should have generated x? Denote by fj (x) the densities of each population Dj , respectively.
The maximum likelihood principle is implemented by allocating x to DJ so that
J = arg max fj (x).
j∈{1,··· ,g}

Consider g = 2 univariate normal populations D1 = N (µ1 , σ12 ) and D2 = N (µ2 , σ22 ). Then
   
1 (x − µ1 )2 1 (x − µ2 )2
f1 (x) > f2 (x) ⇐⇒ exp − > exp −
σ1 2σ12 σ2 2σ 2
     2 2 2  
1 1 µ1 µ2 µ1 µ2 σ2
⇐⇒ − x 2
− 2 − x + − < 2 log ,
σ12 σ22 σ12 σ22 σ12 σ22 σ1
which simplifies to  
µ1 + µ2
(µ1 − µ2 ) x − >0
2
provided that σ1 = σ2 .

Figure 8: Maximum likelihood rule for univariate normal distributionsxxviii)

In multivariate setting, suppose Dj = Np (µj , Σ). Note that the covariance matrices are equal for all
classes. Then
f1 (x) > f2 (x) ⇐⇒ (x − µ1 )′ Σ−1 (x − µ1 ) < (x − µ2 )′ Σ−1 (x − µ2 )
 

⇐⇒ log ff21 (x)
(x)
= (µ1 − µ2 )′ Σ−1 x − 12 (µ1 + µ2 ) > 0.
p
Note that kx − µkΣ = (x − µ)′ Σ−1 (x − µ) is called the Mahalanobis distance. The classification rule
based on linear combinations belongs to the family of linear discriminant analysis (LDA) methods.

§20 Classification Cont’d (2019/12/17)



The decision boundary x ∈ Rp : (µ1 − µ2 )′ Σ−1 x − 12 (µ1 + µ2 ) = 0 is often indicated in contour plots,
as a hyperplane that goes through the point 12 (µ1 + µ2 ) and has normal vector a = Σ−1 (µ1 − µ2 ).
xxviii) generated by [Link]

30
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Figure 9: LDA decision boundary for bivariate normal distributionsxxix)

For Dj = Np (µj , Σj ) with Σ1 6= Σ2 , we have f1 (x) > f2 (x) if and only if

( ) [ ( ) ]
f1 (x) 1 1 det Σ1 1 ( ′ −1 )
log = − x′ (Σ−1
1 − Σ−1
2 )x + (µ ′ −1
1 Σ1 − µ ′ −1
2 Σ2 )x − log + µ 1 Σ 1 µ 1 − µ ′ −1
2 Σ2 µ 2 > 0.
f2 (x) 2 2 det Σ2 2

This classification rule belongs to the family of quadratic discriminant analysis (QDA) methods.
Denote the classifier by G, i.e.,
G = j ⇐⇒ X ∼ Dj .
In Bayesian view, past experience gives rise to prior probabilities

πj = p(Dj ) = P(G = j), j ∈ {1, · · · , g}

before data is collected. Then, Bayes’ theorem allows us to update our previous beliefs with new information.
Maximization of the posterior probability

πj fj (x)
p(Dj |x) = P(G = j|x) = Pg ∝ πj fj (x)
k=1 πk fk (x)

leads to the naïve Bayesian classifier (see also Lecture §3)

Ĝ(x) = arg max p(Dj |x) = arg max πj fj (x).

j∈{1,··· ,g} j∈{1,··· ,g}

Note that  
 
f1 (x) π2
π1 f1 (x) > π2 f2 (x) ⇐⇒ log > log .
f2 (x) π1
The Bayes rule of discrimination will coincide with the maximum likelihood principle if, πj ∝ 1 correspond
to the uniform distribution, i.e., the prior is non-informative.
To apply these methods, one gathers a training sample, for which the correct classes are known (labeled).
Then it suffices to estimate parameters as follows:
• µj ← x̄j , sample mean of observations that belong to class j.
• Σ ← Spooled , sample covariance matrix of pooled observations, applicable to LDA.
• Σj ← Sj , sample covariance matrix of observations that belong to class j, applicable to QDA.
• πj ← nj /n, proportion of observations in the data that belong to class j.
We may encounter a misclassification error when the observation is assigned to a certain group. Denote
by Aj ⊂ Rp the classification region (a.k.a. the action space in statistical decision theory) such that we
allocate x ∈ Aj to Dj . Let c(k|j) ∈ R+ be the cost when an observation actually from Dj falls into Ak and
is mistakenly assigned to Dk . The probability of putting X ∼ Dj into Dk can be calculated as
Z
p(k|j) = P(X ∈ Ak |Dj ) = fj (x) dx.
Ak

xxix) from [Link]

31
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Suppose there are only two populations. The expected cost of misclassification is then

ECM = c(2|1)p(2|1)π1 + c(1|2)p(1|2)π2 .

The rule minimizing the ECM among all partitions A1 t A2 = Rp is given by (check!)
   
f1 (x) c(1|2)π2
x ∈ A1 ⇐⇒ c(2|1)π1 f1 (x) > c(1|2)π2 f2 (x) ⇐⇒ log > log .
f2 (x) c(2|1)π1
This is exactly the Bayes discrimination rule as long as the misclassification costs are equal. Thus, the
Bayes rule is optimal in that it minimizes the total probability of misclassification

TPM = p(2|1)π1 + p(1|2)π2 .

Note that if two or more classification rules are available, it is important to evaluate the performance of each
rule. The smallest possible TPM is called the optimum error rate (OER). Let us derive an expression
for the OER when π1 = π2 = 12 and Dj = Np (µj , Σ). LDA yields
 
A1 = x ∈ Rp : a′ x > 12 a′ (µ1 + µ2 ) , & A2 = x ∈ Rp : a′ x ≤ 12 a′ (µ1 + µ2 ) ,

where a = Σ−1 (µ1 − µ2 ). After tedious calculations, one can see that

p(2|1) = p(1|2) = Φ(−∆/2) =⇒ OER = Φ(−∆/2),

√
where ∆ = a′ Σa = kµ1 − µ2 kΣ . Generally speaking, estimated TPM is bitterly subtle and involved.

§21 Classification Cont’d (2019/12/24)

The performance of classifiers can be empirically evaluated by calculating the estimate of TPM, i.e.,
Z Z
[ = π̂1 p̂(2|1) + π̂2 p̂(1|2) = n1
TPM fˆ1 (x) dx +
n2
fˆ2 (x) dx
n Â2 n Â1

using a sample of size n = n1 + n2 .

The probabilities of misclassification may also be estimated by the re-substitution method, which does
not depend on the form of the parent populations. We reclassify the testing sample according to a rule
obtained from the training sample. Then we have p̂(k|j) = njk /nj as an estimate of p(k|j), where njk is
the number of individuals coming actually from Dj that are classified into Dk .

Estimated
class 1 class 2
class 1 n11 n12 n1
True
class 2 n21 n22 n2

Table 4: Confusion matrix in binary classification

The estimated TPM is called error rate here, defined as

X
g X Xg X
njk nj 1X
g
[ =
TPM p̂(k|j)π̂j = =1− njj .
j=1 k̸=j j=1
nj n n j=1
k̸=j

[ all tends to be accompanied by downward bias, which

If all samples are used to train a classifier, then TPM
means that the classification rule appears overly optimistic.
An approach that corrects for this bias is based on the holdout procedure. The original observations
are split into a training sample and a testing sample recurrently, and the cross-validation (CV, see also
Lecture §16) with the 0-1 loss applies. The leave-one-out method, which is a special case of K-fold CV,
[ CV . Thereby, we usually have TPM
is frequently used and often results in an unbiased TPM [ CV > TPM[ all .
For example, see exercise 2 in homework 5.

32
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

By the way, data are often partitioned into three parts if there exist tuning parameters. The tuning
sample is used to numerically assess and determine parameters first. The training and testing sample are
then used to compare the proposed method with others.
Fisher’s discriminant analysis is a distance-based method, which is extensively applied due to its
geometrical interpretation. The idea is to find a linear rule (in other words, base the classification rule on
a projection, say, x ∈ Rp 7→ a′ x ∈ R) that separates the classes best. Denote the observations by xij ∈ Rp ,
i = 1, · · · , nj ; j = 1, · · · , g. Let yij = a′ xij for some a ∈ Rp . In the case when g = 2, the objective is
!
|ȳ·1 − ȳ·2 | 1 Xn1 Xn2
max q 2
, where sy = (yi1 − ȳ·1 ) +
2
(yi2 − ȳ·2 ) .
2
a∈Rp
s2 n1 + n2 − 2 i=1 i=1
y

Note that Fisher’s approach does not assume that the populations are normal. It does, however, implicitly
assume that the population covariance matrices are equal. Write
!
1 Xn1 Xn2
′ ′
Sx = (xi1 − x̄·1 )(xi1 − x̄·1 ) + (xi2 − x̄·2 )(xi2 − x̄·2 ) .
n1 + n2 − 2 i=1 i=1

Then q
|ȳ·1 − ȳ·2 | |a′ (x̄·1 − x̄·2 )|
q = √ ′ ≤ (x̄·1 − x̄·2 )′ Sx−1 (x̄·1 − x̄·2 ) = kx̄·1 − x̄·2 kSx
s2y a Sx a

by Cauchy-Schwarz inequality, with equality holding if and only if a ∝ â = Sx−1 (x̄·1 − x̄·2 ), which has

Figure 10: Fisher’s linear discriminantxxx)

already occurred in Gaussian LDA. An observation x ∈ Rp will be allocated to D1 if and only if

|â′ (x − x̄·1 )| < |â′ (x − x̄·2 )| ⇐⇒ â′ x − 12 (x̄·1 + x̄·2 ) > 0.

Generally, we seek a ∈ Rp to maximize the ratio of the between-group variance to the within-group variance,
for g classes with equal covariance matrices. Let (see Lecture §12, table 1)
X
g X
g X
nj
B= nj (x̄·j − x̄·· )(x̄·j − x̄·· )′ , & W = (xij − x̄·j )(xij − x̄·j )′ .
j=1 j=1 i=1

Note that Spooled = W /(n − g). Then,

X
g X
g X
nj
nj (ȳ·j − ȳ·· )2 = a′ Ba, & (yij − ȳ·j )2 = a′ W a.
j=1 j=1 i=1

xxx) from Bishop’s PRML; see also [Link]

33
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

One can show that

a′ Ba (b=W 1/2 a) b′ W −1/2 BW −1/2 b
arg max ==========
a∈Rp a′ W a b′ b
is given (check!) by the first eigenvector â1 of W −1 B, producing the first Fisher’s sample discriminant.
Now a classification rule is easy to obtain: we allocate x ∈ Rp to class J1 where

J1 = arg min |â′1 (x − x̄·j )| .

j∈{1,··· ,g}

We can utilize more linear discriminants in order to summarize the distinction between the classes in higher
dimension. Let
âk = arg max a′ Ba s.t. a′ Spooled a = 1; a′ Spooled âℓ = 0, ∀ℓ < k.
a∈Rp
−1
It turns out that âk is the k th eigenvector of Spooled B for 1 ≤ k ≤ p. The number of âk ’s corresponding to
nonzero eigenvalues is

s = rank W −1 B ≤ min(p, g − 1),
Pg
since B is spanned by those x̄·j − x̄·· such that j=1 nj (x̄·j − x̄·· ) = 0p . Also, for any k > s,

Bâk = 0p =⇒ â′k Bâk = 0 =⇒ â′k (x̄·j − x̄·· ) = 0, ∀j = 1, · · · , g.

Allocate x ∈ Rp to class
X
r
|â′k (x − x̄·j )| .
2
Jr = arg min
j∈{1,··· ,g} k=1

Note that !
X
p X
p
|â′k (x ′
âk â′k
2 2
− x̄·j )| = (x − x̄·j ) (x − x̄·j ) = kx − x̄·j kSpooled .
k=1 k=1
−1
= Spooled
Pp
If r < s discriminants are used for classification, there is a loss of squared distance of k=r+1 |â′k (x − x̄·j )| ,
2
Ps ′ 2 ′
where k=r+1 |âk (x − x̄·j )| is the part useful for classification. Note that âk (x̄·j ′ − x̄·j ) = 0 for k > s, and
Pp
therefore k=s+1 |â′k (x − x̄·j )| is constant — does not depend on j.
2

§22 Classification Cont’d & Final Review (2019/12/26)

If prior probabilities πj are considered, we may allocate x ∈ Rp to class
( r )
X
′ 2
arg min |âk (x − x̄·j )| − 2 log(πj ) ,
j∈{1,··· ,g} k=1

or equivalently, ( )
1X ′
r
2
arg max − |âk (x − x̄·j )| + log(πj ) .
j∈{1,··· ,g} 2
k=1

Assume that there are only two classes, indexed by an indicator Y that takes the value from {0, 1}. In
logistic regression, one directly models

p(x) := P(Y = 1|X = x) = E[Y |X = x].

 
The logarithm of the odds p
1−p is called the logit of the probability p, i.e., logit : p ∈ [0, 1] 7→ log p
1−p ∈ R.
xxxi)
It is a special case of a link function in a generalized linear model (GLM ) — for the Bernoulli
distribution. As in linear regression,
logit(p(x)) = α + x′ β.
xxxi) see
[Link] for a brief (p)review. More generally, when an exact distribution is difficult
to handle, we often appeal to corresponding large sample properties.

34
© RandomSeminar
 Notes for Applied Multivariate Statistical Analysis

Then
exp(α + x′ β)
p(x) = .
1 + exp(α + x′ β)
It follows that the log-likelihood function is
! n h
Y
n X
i
ℓn (α, β) = log [p(xi )]yi [1 − p(xi )]1−yi = (α + x′i β)yi − log 1 + exp(α + x′i β) .
i=1 i=1

To obtain the MLE for θ = (α, β), we often solve

 Pn  
 ∂ℓn exp(α+x′i β)

 ∂α = yi − 1+exp(α+x′i β) =0
i=1
 ∂ℓ Pn  

 n = yi −
exp(α+x′i β)
xi = 0p
∂β 1+exp(α+x′i β)
i=1

numerically; for instance, the Newton–Raphson method is iterated as

θ̂(t+1) = θ̂(t) − ℓ̈n (θ̂(t) )−1 ℓ̇n (θ̂(t) ).
Note that ℓ̈n (θ) depends on θ and xi ’s, but not on yi ’s. In contrast to LDA, logistic regression relies on
conditional likelihood and requires a large sample size, which means that logistic regression is less efficient.
LDA, relying only on the equal covariance matrices assumption, is less robust.
Nonparametric classification is based mainly onP density estimation. A generalization of the kNN
1 n xi −x
estimate (see Lecture §3) is given by p̂n (x) = nh i=1 K( h ), called the kernel density estimate,
where K is a kernel and h is the bandwidth.

 As for the final exam, important topics are

 eigen-decomposition and SVD (see Lecture 1,2)
 multivariate normal distr — conditional distr, normal matrix, Wishart distr (see Lecture 5,6,8,9)
 inference about means — one sample, two samples (see Lecture 9–11)
 Cochran’s Theorem — sum of squares decomposition (see Lecture 9)
 multivariate linear regression — basics (see Lecture 13,14)
 PCA and FA — statistical interpretation (see Lecture 16–18)
 classification — LDA applied to two populations (see Lecture 19–21)
The End K

35
© RandomSeminar