First Problem Set Multivariate Statistical Analysis

Edoardo Bandoni, Ecem Akın

April 7, 2024

Exersise 1
The air pollution data (file [Link]) consists of 7 measurements recorded at n = 41 cities in
the United States. Variables are SO2: Sulphur dioxide content in micrograms per cubic meter,
[Link]: Average annual temperature in F o (negative values), Manuf : Number of manufac-
turing enterprises employing 20 or more workers, Pop: Population size (1970 census) in thousands,
Wind: Average annual wind speed in miles per hour, Precip: Average annual precipitation in
inches, Days: Average number of days with precipitation per year.

u s a i r = r e a d . t a b l e ( ” data / u s a i r . t x t ” , h e a d e r=TRUE)
head ( u s a i r )

SO2 [Link] Manuf Pop Wind Precip Days

Phoenix 10 -70.3 213 582 6.0 7.05 36
Little Rock 13 -61.0 91 132 8.2 48.52 100
San Francisco 12 -56.7 453 716 8.7 20.66 67
Denver 17 -51.9 454 515 9.0 12.95 86
Hartford 56 -49.1 412 158 9.0 43.37 127
Wilmington 36 -54.0 80 80 9.0 40.25 114

We are going to ignore the SO2 variable and concentrate on the remaining 6, two of which relate
to human ecology (Manuf, Pop) and four to climate ([Link], Wind, Precip, Days).

airquality = airquality [ ,2:7]

1. Correlation
We compute the correlation matrix by using cor(): when applied to a data frame like usair with
6 columns, it returns a 6-dimensional square matrix with sample correlations computed on each
element of the [Link]. With 6 variables it is difficult to see what the pairs of variables with
the largest correlation are, although we do expect variable which relate to the same matter to be
highly correlated.

R = cor ( a i r q u a l i t y )
R[ ! l o w e r . t r i (R) ] = NA
round (R, 3 )

[Link] Manuf Pop Wind Precip Days

[Link] NA NA NA NA NA NA
Manuf 0.190 NA NA NA NA NA
Pop 0.063 0.955 NA NA NA NA

1
Wind 0.350 0.238 0.213 NA NA NA
Precip -0.386 -0.032 -0.026 -0.013 NA NA
Days 0.430 0.132 0.042 0.164 0.496 NA

s o r t ( abs (R) , d e c r e a s i n g=T)

[1] 0.95526935 0.49609671 0.43024212 0.38625342 0.34973963 0.23794683 0.21264375 0.19004216 0.1641
[10] 0.13182930 0.06267813 0.04208319 0.03241688 0.02611873 0.01299438

o r d e r ( abs (R) , d e c r e a s i n g=T, na . l a s t = NA)

[1] 9 30 6 5 4 10 16 2 24 12 3 18 11 17 23

p = dim ( a i r q u a l i t y ) [ 2 ]
s t a c k .m = matrix ( 1 : p ˆ 2 , n c o l=p ) ; s t a c k .m

[,1] [,2] [,3] [,4] [,5] [,6]

[1,] 1 7 13 19 25 31
[2,] 2 8 14 20 26 32
[3,] 3 9 15 21 27 33
[4,] 4 10 16 22 28 34
[5,] 5 11 17 23 29 35
[6,] 6 12 18 24 30 36

o r d e r . c o r = o r d e r ( abs (R) , d e c r e a s i n g=T, na . l a s t = NA)

which ( s t a c k .m == o r d e r . c o r [ 1 ] , a r r . i n d=T)

row col
[1,] 3 2

my . fun = f u n c t i o n ( x ) which ( s t a c k .m == x , a r r . i n d=T)

my . fun = V e c t o r i z e (my . fun )
my . fun ( o r d e r . c o r [ 1 : 3 ] )

[,1] [,2] [,3]

[1,] 3 6 6
[2,] 2 5 1

[Link]() returns a logical matrix with true values below the diagonal, while ![Link]() has
true above and on the diagonal. Since we do not want to order the correlations with the ones
in the diagonal and the repetitions above and below the diagonal, we put NA on the diagonal
and the upper triangular part. Firsty, we obtain the correlations in decreasing order via sort(),
then the permutation which arranges the correlations in descending order via order(). It is a
nuisance that we also get the NA in this ordering, so we tell R to avoid put NA last by means
of [Link]=NA. We also have to cope with the fact that the order we get refers to the entries of
the correlation matrix arranged as a long vector obtained by stacking up the columns. We create
the p × p matrix stack.m filled with 1, 2, . . . , p2 column wise, and use it to convert the order into
pairs of indexes, row and col. which(stack.m==[Link][1]) would return the position of the
largest correlation as a single integer, so we add the argument [Link]=T so to get the row and
col indexes. The last three lines of the code are meant to make this last operation in a vectorized
way. See below for a more polished printing of the ranking of the first 6 highest correlations.

2
Order . c o r = t (my . fun ( o r d e r . c o r ) )
out = matrix ( names ( a i r q u a l i t y ) [ Order . c o r ] , n c o l =2)
colnames ( out ) = c ( ”row” , ” c o l ” )
out = data . frame ( out , c o r=round ( s o r t ( abs (R) , d e c r e a s i n g=T) , 3 ) )
head ( out )

row col cor

1 Pop Manuf 0.955
2 Days Precip 0.496
3 Days [Link] 0.430
4 Precip [Link] 0.386
5 Wind [Link] 0.350
6 Wind Manuf 0.238

The results confirm what we guessed at the beginning, the correlation between Pop and Manuf
is almost 1, furthermore the variables regarding climate are rather correlated among them. 6th
highest observation is an exception, which is between a variable regarding ecology and one regarding
climate.

2. Univariate outliers with boxplot

For univariate outliers, we calculate the standardized values and examine them for large or small
values compared to the quantiles of standard normal distribution. Standardized values are com-
puted via the scale() function, which returns the sample means and the standard deviations, so
to extract the matrix of standardized values we use the subsetting [ , ].

head ( round ( s c a l e ( a i r q u a l i t y ) [ , ] , 2 ) )

[Link] Manuf Pop Wind Precip Days

Phoenix -2.01 -0.44 -0.05 -2.41 -2.52 -2.94
Little Rock -0.72 -0.66 -0.82 -0.87 1.00 -0.52
San Francisco -0.13 -0.02 0.19 -0.52 -1.37 -1.77
Denver 0.53 -0.02 -0.16 -0.31 -2.02 -1.05
Hartford 0.92 -0.09 -0.78 -0.31 0.56 0.49
Wilmington 0.24 -0.68 -0.91 -0.31 0.30 0.00

By means of the logical evaluation abs(scale(airquality)[,])>qnorm(1-alpha) we identify out-

lying standardized values with respect to the α-quantile of the standard normal distribution. Ob-
servation index and variable index are then obtained by using which( , [Link]=T). The usual
choice of the quantile α is 0.5/n, where n is the number of observation in the dataset: in such way
we handle the problem of having a limited number of observations.

which ( abs ( s c a l e ( a i r q u a l i t y ) [ , ] )>qnorm (1 −0.5 /n ) , a r r . i n d=T)

row col
Miami 9 1
Chicago 11 2
Chicago 11 3
Philadelphia 29 3
Phoenix 1 4
Wichita 14 4

3
Phoenix 1 5
Alburquerque 23 5
Phoenix 1 6

Alternatively, we produce boxplots of each variable and identify the index of observations
flagged as outliers, if any, that are the points outside of the whiskers of the boxplot. Since the
variables are measured on a different scale, we make a boxplot for every variable in different figures.
In order to use ggplot2 effectively, we need your data in a ”long” format where each row
represents a single observation. This format is often preferred because it allows ggplot2 to easily
map variables to aesthetics in the plot layers. The gather() function from the tidyr package in R
is commonly used to reshape data from wide to long format, making it compatible with ggplot2’s
grammar of graphics.

a i r q u a l i t y 2 = g a t h e r ( a i r q u a l i t y , key = ” V a r i a b l e ” , v a l u e = ” Value ” )

Variable Value
1 [Link] -70.3
2 [Link] -61.0
3 [Link] -56.7
4 [Link] -51.9
5 [Link] -49.1
6 [Link] -54.0

We are now able to plot the 6 boxplots each referring to a single variable. Since the commands
will be essentially the same for each variable, first we create an index variable j and use it to select
the variable.

4
j = 0
u = 41 ∗ j + 1
b = 41 ∗ ( j + 1 )
g g p l o t ( a i r q u a l i t y 2 [ u : b , ] , a e s ( x=V a r i a b l e , y=Value ) ) +
geom b o x p l o t ( o u t l i e r . c o l o u r=” r e d ” , o u t l i e r . shape =16 ,
o u t l i e r . s i z e =2) +
labs (x = ” Variables ” ,
y = ” Measuraments ” ) +
theme g r e y ( )

u and b identify the rows we are taking into account in order to make the boxplot of the variable
we are plotting and airquality2[u:b,] are the observations of a single variable. With the command
geom boxplot() we create the boxplot figure and we colour in red the outliers with specific size
and shape in order to make them more visible. Finally, labs() allows us to give titles.
The following function identify extreme outliers allow us to detect, for every variable, what
the two extreme outliers are.

i d e n t i f y extreme o u t l i e r s = f u n c t i o n ( x , max o u t l i e r s = 2 ) {
q = q u a n t i l e ( x , probs = c ( 0 . 2 5 , 0 . 7 5 ) )
iqr = q [2] − q [1]
l o w e r bound = q [ 1 ] − 1 . 5 ∗ i q r
upper bound = q [ 2 ] + 1 . 5 ∗ i q r
extreme o u t l i e r s = which ( x < l o w e r bound | x > upper bound )
w h i l e ( l e n g t h ( extreme o u t l i e r s ) > max o u t l i e r s ) {
iqr = iqr ∗ 1.1
l o w e r bound = q [ 1 ] − 1 . 5 ∗ i q r
upper bound = q [ 2 ] + 1 . 5 ∗ i q r
extreme o u t l i e r s = which ( x < l o w e r bound | x > upper bound )
}
r e t u r n ( extreme o u t l i e r s )
}

We saves in extreme outliers the univariate outliers identified for x; at this point, until these are
more than 2, we increase the iqr by 10% and recalculate the outliers. Finally, Once the loop exits
(when the number of outliers is within the specified limit), the function returns the indices of the
identified extreme outliers. In summary, this function iteratively adjusts the interquartile range
(IQR) to limit the number of outliers to at most max outliers. It starts with the conventional
1.5 * IQR rule but increases the IQR if necessary until the desired number of outliers is achieved.

i n d e x Neg . Temp = i d e n t i f y extreme o u t l i e r s ( a i r q u a l i t y [ , 1 ] )

i n d e x Manuf = i d e n t i f y extreme o u t l i e r s ( a i r q u a l i t y [ , 2 ] )
i n d e x Pop = i d e n t i f y extreme o u t l i e r s ( a i r q u a l i t y [ , 3 ] )
i n d e x Wind = i d e n t i f y extreme o u t l i e r s ( a i r q u a l i t y [ , 4 ] )
i n d e x P r e c i p = i d e n t i f y extreme o u t l i e r s ( a i r q u a l i t y [ , 5 ] )
i n d e x Days = i d e n t i f y extreme o u t l i e r s ( a i r q u a l i t y [ , 6 ] )
i n d e x Neg . Temp ; i n d e x Manuf ; i n d e x Pop ; indexWind ; i n d e x P r e c i p ; i n d e x
Days

[1] 9
[1] 11 29
[1] 11 29

5
integer(0)
[1] 1 23
[1] 1 23

The result shows us that the observations 9, 11, 29, 1, 23 are univariate outliers. Let’s colore-code
them and see what cities are actually univariate outliers.

c o l . ind = rep ( ” black ” ,n)

c o l . ind [ 9 ] = ” red ”
c o l . ind [ 1 1 ] = ” blue ”
c o l . ind [ 2 9 ] = ” green ”
c o l . i n d [ 1 ] = ”brown”
c o l . i n d [ 2 3 ] = ” magenta ”
o u t l i e r s = c (1 , 9 , 11 , 23 , 29)
rownames ( a i r q u a l i t y [ o u t l i e r s , ] )

[1] "Phoenix" "Miami" "Chicago" "Alburquerque" "Philadelphia"

3. Marginal normality with QQ plot

We examine marginal normality by producing normal qq-plot for each variable. Since the com-
mands will be essentially the same for each variable, we firstly create an index variable j, and use
it to select the variable as well as to specify the title of the plot. ”Using ggplot(), a new ggplot
object is initialized with the dataframe df, containing the sorted values of x. The sample aesthetic
is mapped to the y variable. Next, stat qq() is employed to incorporate points onto the plot,
representing the quantiles of the sorted data. The color aesthetic is utilized to color the points
based on previously identified outliers. Lastly, stat qq line() adds a dashed Q-Q line to the plot,
providing further insight into the relationship between theoretical and observed quantiles.
We conclude for normality when the points lie very nearly along a straight line, which seems to
be the case of Days, Precip, Wind and [Link]. On the other hand, Manuf and Pop seem
to be farther to the QQline, hence it is reasonable that these variables follow a similar pattern
since they where highly correlated.

6
j = 1
x = airquality [ , j ]
sorted order = order (x)
d f = data . frame ( y = s o r t ( x ) )
g g p l o t ( df , a e s ( sample = y ) ) +
s t a t qq ( c o l o r = c o l . i n d [ s o r t e d o r d e r ] ) +
s t a t qq l i n e ( l i n e t y p e = ” dashed ” ) +
l a b s ( x = ” T h e o r e t i c a l Q u a n t i l e s ” , y = ” Sample Q u a n t i l e s ” , t i t l e = ”QQ
P l o t Neg . Temp” ) +
theme g r e y ( )

4. Outliers with scatterplot

In order to detect outliers from the scatter plots, we make a scatter plot matrix through pairs()
replacing the plotting character with the observation index.

n = nrow ( a i r q u a l i t y )
p a i r s ( a i r q u a l i t y , p a n e l=f u n c t i o n ( x , y ) t e x t ( x , y , l a b e l s =1:n ) , l o w e r . p a n e l
= NULL)

It seems that observations 1, 29, 23, 9 stand out from the other ones in a few plots. We create a
vector of colors of the outliers identified at point 2, so to plot these observations in red, blue and

7
green,, brown, orange and magenta. We redo the scatter plot matrix using full points (pch=16)
of smaller font (cex=0.7) color-coding the observations identified at point 2.
p a i r s ( a i r q u a l i t y , pch =16 , cex =0.7 , c o l=c o l . ind , l o w e r . p a n e l=NULL)

Formally, bivariate outliers are not necessarily univariate outliers like observations 1, 29, 11, 23
and 29, they are more in general points that sit far off the main bulk of the data, i.e. far off
the ellipsoidal representation of the scatter plots. It is so useful to reproduce the scatter plots
with ellipses color coding the observations. Since the commands will be essentially the same for
each couple of variables, we firstly create two indexes variable j and k, and use them to select the
variable as well as to specify the title of the plot.
i n s t a l l . p a c k a g e s ( ” sp ” )
l i b r a r y ( sp )
j = 1; k = 2
x = airquality [ , j ]
y = airquality [ ,k]
data = data . frame ( x = x ,
y = y)
e l l i p s e data = e l l i p s e : : e l l i p s e ( cov ( data ) , c e n t e r = colMeans ( data ) ,
level = 0.95)
e l l i p s e p o l y = S p a t i a l P o l y g o n s ( l i s t ( Polygons ( l i s t ( Polygon ( e l l i p s e data )
) , ID = ” e l l i p s e ” ) ) )
o u t s i d e p o i n t s = ! p o i n t . i n . polygon ( data $x , data $y , e l l i p s e
poly@polygons [ [ 1 ] ] @Polygons [ [ 1 ] ] @coords [ , 1 ] , e l l i p s e poly@polygons
[ [ 1 ] ] @Polygons [ [ 1 ] ] @coords [ , 2 ] )
new o u t l i e r s = union ( o u t l i e r s , which ( o u t s i d e p o i n t s ) )
l a b e l d f = data . frame (
x = x [ new o u t l i e r s ] ,
y = y [ new o u t l i e r s ] ,
l a b e l = row . names ( a i r q u a l i t y ) [ new o u t l i e r s ]
)
g g p l o t ( data , a e s ( x = x , y = y ) ) +

8
 geom p o i n t ( c o l o u r = c o l . ind , shape = 1 6 ) +
geom path ( data = a s . data . frame ( e l l i p s e data ) ) +
geom t e x t ( data = l a b e l df , a e s ( x = x , y = y , l a b e l = l a b e l ) , s i z e =
2.5) +
l a b s ( x = names ( a i r q u a l i t y ) [ j ] , y = names ( a i r q u a l i t y ) [ k ] ) +
theme g r e y ( )

We use ellipse data() to define the ellipse with proper center and angle, given by the mean and
the covariance matrix, and the level parameter equal to 0.95. SpatialPolygons() converts the
ellipse data into an object, which can be used for geometric operations, therefore we determine
which points from the data fall outside of the ellipse using the [Link]() function. At
this point we update the outliers making the union between the previous one the new one. Finally,
label df is a data frame containing positions and names of the outliers and geom text is used to
plot the names of the cities. Observations 11 (blue), 1 (brown) appear outside of many ellipses, we
could expect them to be multivariate outlier as soon as bivariate oultiers. Conversely, observation
29 (green) consistently falls within all the ellipses, suggesting it is merely a univariate outlier.
Consequently, it is removed from the outliers list and colored black.
o u t l i e r s = o u t l i e r s [ o u t l i e r s != 2 9 ]
c o l . ind [ 2 9 ] = ” black ”
9 (red) and 23 (magenta) seem to be closer to the ellipses, but they remains outside some of them.
Looking at new outliers we can see that there appear some new, such as Huston, Charleston and

9
Buffalo, but there is an observation present in many figures which is Wichita. For such reason, we
color code only Wichita, which is the observation 14.

outliers = c( outliers ,14)

c o l . ind [ 1 4 ] = ” gold ”

5. Multivariate normality with Q-Q plot squared Mahalanobis distances

For checking multivariate normality by using all variables at once, we first compute the squared
Mahalanobis distances. Then, we produce a qq-plot against the chi-squared distribution with p = 6
degrees of freedoms.

d = mahalanobis ( a i r q u a l i t y , c e n t e r = colMeans ( a i r q u a l i t y ) , cov = cov (

airquality ) )
sorted order = order (d)
d f = data . frame ( y = s o r t ( d ) )
g g p l o t ( df , a e s ( sample = y ) ) +
s t a t qq ( c o l o r = c o l . i n d [ s o r t e d o r d e r ] ) +
s t a t qq l i n e ( l i n e t y p e = ” dashed ” ) +
l a b s ( x = ” T h e o r e t i c a l Q u a n t i l e s ” , y = ” Sample Q u a n t i l e s ” , t i t l e = ”
Chisq Q−Q P l o t ” ) +
theme g r e y ( )

d = mahalanobis(airquality, center = colMeans(airquality), cov = cov(airquality))

computes the Mahalanobis distance for each observation, where the center parameter specifies the
center of the data (calculated using colMeans()), and the cov parameter specifies the covariance
matrix of the data (calculated using cov()). We create a vector of indexes sorted order that
saves the positions of the ordered Mahalanobis distances that allow us to color the points in a right
way. Finally, we add the diagonal line via stat qq line(). The data do not deviate significantly
from a normal distribution, this suggests that the distribution of the data match sufficiently with
the theoretical normal distribution, apart from observations 1 (brown) and 11 (blue) that could
be multivariate outliers.

10
6. Multivariate outliers
We look for multivariate outliers by inspecting the squared Mahalanobis distances d2 , comparing
them with quantiles of the chi-squared distribution with 6 degrees of freedom. We can make the
comparison via the T 2 chart: we plot d2i against the observation index i and we draw horizontal
lines at eight χ2p,α for α small enough. It should be sufficiently small compared to the sample size
n. In this case, we use different values of α, among them the aformentioned α = 0.5/n, in this way
we have a wider perspective of the possible outliers.
d f = data . frame ( x = 1 : l e n g t h ( d ) , y = d )
l a b e l d f = data . frame (
x = d f $x [ o u t l i e r s ] ,
y = d f $y [ o u t l i e r s ] ,
l a b e l = row . names ( a i r q u a l i t y ) [ o u t l i e r s ] )
g g p l o t ( data = df , a e s ( x = x , y = y ) ) +
geom p o i n t ( c o l o u r = c o l . ind , shape = 1 6 ) +
geom h l i n e ( y i n t e r c e p t = q c h i s q (1 −0.05 , d f = p ) , l i n e t y p e = ” dashed ” )
+
geom h l i n e ( y i n t e r c e p t = q c h i s q (1−1/ l e n g t h ( d ) , d f = p ) , l i n e t y p e = ”
dashed ” ) +
geom h l i n e ( y i n t e r c e p t = q c h i s q (1 −0.5 / l e n g t h ( d ) , d f = p ) , l i n e t y p e = ”
dashed ” ) +
geom t e x t ( data = l a b e l df , a e s ( x = x , y = y , l a b e l = l a b e l ) , s i z e =
2.5) +
l a b s ( x = ” Index ” , y = ” Values ” ) +
theme g r e y ( )

We use ggplot2 to create a plot with the Mahalanobis distances on the y-axis and the indices of the
observations on the x-axis. geom point() adds points to the plot representing the Mahalanobis
distances, where the colour parameter is specifies the univariate outliers identified in point 2.
geom hline() adds dashed horizontal lines to the plot at three specific quantiles of the chi-square
distribution respectively, from the gratest to the smallest, α = 0.95, α = 1 − 1/lenth(d) and

11
α = 1−0.5/length(d), where length(d) is the number of observations. Observations 1 and 11 can
be considered multivariate outliers, because they appear over the dashed lines, while observations
9, 23 and 28 cannot be considered multivariate outliers.

12
Exercise 2
Let X = (X1 , . . . , Xp ) be distributed according to

X =µ+VZ +ϵ

where µ ∈ R, Z ∼ Nq (0, σz2 Iq ), q < p, V a p × q matrix with q orthogonal unit vectors as columns
and ϵ ∼ Nq (0, σ 2 Ip ). Here Iq and Ip are identity matrices of dimension q × q and p × p, respectively.
Let also Z and ϵ be independent random vectors and σz2 = (1 + δ)σ 2 for δ > −1.

1.
Firstly, let us find the distribution of X, proceeding by steps:
   
V Z ∼ Np V 0, V σz2 I q V T = Np 0, σz2 V V T

   
V Z + ϵ ∼ Np 0, σz2 V V T + σ 2 I p = Np 0, σ 2 ((1 + δ)V V T + I p )
 
X = µ + V Z + ϵ ∼ Np µ, σ 2 ((1 + δ)V V T + I p ) = Np (µ, Σ)

Observe that, since the columns of V are orthogonal, they are also linear independent, furthermore
the rank of V is equal to q < p. At this point, observe that V V T is an orthogonal projection of a
subspace of dimension q in a space of dimension p. Indeed,

(V V T )(V V T ) = V (V T V )V T = V V T

where V T V = Iq , because the columns of V are orthogonal, furthermore V V T is linear and finally
(V V T )T = V V T . Recall that orthogonal projections have eigenvalues equal to 0 or with module
equal to 1. Indeed, let P be an orthogonal projection and v an eigenvector of P with eigenvalue
λ, then
P 2 v = P λv = λP v = λ2 v

and at the same time

P 2 v = P v = λv

and hence λ2 = λ, that is λ = 0 or λ = 1. Since V V T has rank q, it has exactly q eigenvalues

equal to 1 and p − q equal to 0. Let us now compute the rotation matrix, denote with vj the j-th
column of V , then
(V V T )vj = V (V T vj ) = V ej = vj

where ej is the j-th vector of the standard base and V T vj = ej because it is the scalar product
between the columns of of V that by hypothesis are orthogonal. If we now complete the set of
vectors {v1 , . . . , vq } to an orthonormal basis of Rp , e.g. B = {v1 , . . . , vq , u1 , . . . , up−q } and we
construct the matrix U , whose columns are the vectors of B, we obtain the spectral decomposition
of V V T as " #
T T T Iq 0q,p−q
U VV U =U U
0p−q,q 0p−q,p−q
where 0i,j is a i × j matrix with zero entries.
Therefore, we obtain the principal component decomposition by
   
U T ΣU = U T σ 2 (1 + δ)V V T + I p U = σ 2 (1 + δ)U T V V T U + U T U =

13
 " #
2

T T
 σ 2 (2 + δ)I q 0
=σ (1 + δ)U V V U + I p = ,
0 σ 2 I p−q

Finally, we are able to find the values of δ such that the first q (population) principal components
of X account for more than 80% of total variation. The percentage of variation given by the first
q principal components is given by

qσ 2 (2 + δ)
≥ c = 0.8.
qσ 2 (2 + δ) + (p − q)σ 2

Let’s solve the inequality with respect to δ, obtaining:

q(δ + 2) ≥ cq(δ + 2) + cp − cq

2cq + cp − cq − 2q q(c − 2) + pc 4p − 6q
δ≥ = = .
q − cq q(1 − c) q

2.
Let δ = 2, σ 2 = 1/3 and
 
−1 2
1
V = 2 −1

3
2 2

We are looking for the distribution of (X1 , X2 )|X3 = −1. Recall that, if
" # " # " #!
X1 µ1 Σ11 Σ12
X= ∼ Np , ,
X2 µ2 Σ21 Σ22

where X is p-dimensional, X1 q-dimensional and X2 p − q-dimensional, then the distribution of

X1 |X2 = x2 is a q-variate normal distribution with mean µ1 + Σ12 Σ−1
22 (x2 − µ2 ) and covariance
matrix Σ22 − Σ12 Σ−1
22 Σ21 . By point 1 we have an expression that gives us the distribution of X,
hence we can make use of R to find the distribution of (X1 , X2 )|X3 = −1.

V = matrix ( c ( −1 ,2 ,2 ,2 , −1 ,2) , 3 , 2 ) / 3
## mi and Sigma a r e t h e mean and t h e c o v a r i a n c e o f X
Sigma = V%∗%t (V) + d i a g ( 3 ) / 3
mi = matrix ( c ( 1 , 1 , 1 ) , n c o l = 1 )
S11 = Sigma [ 1 : 2 , 1 : 2 ]
S12 = Sigma [ 1 : 2 , 3 ]
S21 = Sigma [ 3 , 1 : 2 ]
S22 = Sigma [ 3 , 3 ]
## m and S a r e t h e mean and t h e c o v a r i a n c e o f (X 1 ,X 2 ) c o n d i t i o n e d t o
X 3 = −1
m = mi [ 1 : 2 ] + S12%∗%s o l v e ( S22 )%∗%(−1 − mi [ 3 ] )
S = S11 − S12%∗%s o l v e ( S22 )%∗%S21
m; S

[,1]
[1,] 0.6363636
[2,] 0.6363636
[,1] [,2]
[1,] 0.8484848 -0.4848485
[2,] -0.4848485 0.8484848

14
e l l i p s e data = e l l i p s e : : e l l i p s e ( S , c e n t r e = m, l e v e l = 0 . 9 5 )
ggplot () +
geom p o i n t ( a e s ( x = m[ 1 ] , y = m[ 2 ] ) , c o l o r = ” r e d ” ) +
geom path ( data = a s . data . frame ( e l l i p s e data ) , a e s ( x = x , y = y ) ,
color = ” black ” ) +
l a b s ( x = e x p r e s s i o n (X [ 1 ] ) , y = e x p r e s s i o n (X [ 2 ] ) ) +
theme g r e y ( )

ellipse data = ellipse::ellipse(S, centre = m, level = 0.95) calculates the parameters for
an ellipse based on a covariance matrix S, with the center specified by vector m, and at a confi-
dence level of 0.95. geom path(data = [Link](ellipse data), aes(x = x, y = y),
color = ”black”) adds a path (i.e., the ellipse) to the plot. The data argument specifies the
data frame containing the ellipse data (ellipse data), which is converted to a data frame using
[Link]().

15
Exercise 3
The pen digit data set (file [Link]) was created by collecting 250 samples from 44 writers.
These writers were asked to write 250 digits in random order inside boxes of 500 by 500 tablet
pixel resolution. The raw data on each of n = 10992 handwritten digits consisted of a sequence,
(xt , yt ), t = 1, 2,. . . , T, of tablet coordinates of the pen at fixed time intervals of 100 milliseconds,
where xt and yt were integers in the range 0-500. These data were then normalized to make
the representations invariant to translation and scale distortions. The new coordinates were such
that the coordinate that had the maximum range varied between 0 and 100. Usually xt stays
in this range, because most integers are taller than they are wide. Finally, from the normalized
trajectory of each handwritten digit, 8 regularly spaced measurements, xt ,yt ), were chosen by
spatial resampling, which gave a total of p = 16 variables. The data includes a class attribute,
column digit, coded 0, 1,. . . , 9, about the actual digit.
p e n d i g i t s = r e a d . t a b l e ( ” data / p e n d i g i t s . t x t ” , s e p=” , ” , head=F)
names ( p e n d i g i t s ) = c ( p a s t e 0 ( r e p ( c ( ”x” , ”y” ) , 8 ) , r e p ( 1 : 8 , each =2) ) , ” d i g i t ” )
head ( p e n d i g i t s )

x1 y1 x2 y2 x3 y3 x4 y4 x5 y5 x6 y6 x7 y7 x8 y8 digit
1 47 100 27 81 57 37 26 0 0 23 56 53 100 90 40 98 8
2 0 89 27 100 42 75 29 45 15 15 37 0 69 2 100 6 2
3 0 57 31 68 72 90 100 100 76 75 50 51 28 25 16 0 1
4 0 100 7 92 5 68 19 45 86 34 100 45 74 23 67 0 4
5 0 67 49 83 100 100 81 80 60 60 40 40 33 20 47 0 1
6 100 100 88 99 49 74 17 47 0 16 37 0 73 16 20 20 6

1. Principal Component Analysis

Principal Component Analysis (PCA) is a widely used technique for dimensionality reduction and
data visualization. It aims to transform a set of possibly correlated variables into a new set of
uncorrelated variables called principal components. These components are ordered by the amount
of variance they explain in the original data.
p e n d i g i t s . pca = prcomp ( x = p e n d i g i t s [ , 1 : 1 6 ] , s c a l e = TRUE)
prcomp() is used to perform PCA on the data matrix, where the scale = TRUE argument
standardizes the variables.
summary ( p e n d i g i t s . pca )

Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6
Standard deviation 2.1718 1.7970 1.6052 1.10894 1.03107 0.89294
Proportion of Variance 0.2948 0.2018 0.1610 0.07686 0.06644 0.04983
Cumulative Proportion 0.2948 0.4966 0.6577 0.73452 0.80096 0.85079
PC7 PC8 PC9 PC10 PC11
Standard deviation 0.77888 0.74044 0.64096 0.54611 0.45882
Proportion of Variance 0.03792 0.03427 0.02568 0.01864 0.01316
Cumulative Proportion 0.88871 0.92297 0.94865 0.96729 0.98045
PC12 PC13 PC14 PC15 PC16
Standard deviation 0.33511 0.28369 0.24084 0.18511 0.16673
Proportion of Variance 0.00702 0.00503 0.00363 0.00214 0.00174
Cumulative Proportion 0.98747 0.99250 0.99612 0.99826 1.00000

16
summary([Link]) provides a summary of the PCA results, including the standard devia-
tions (square roots of the eigenvalues) of the principal components and their proportion of variance
with respect to the entire sample.

e i g e n v a l u e s = p e n d i g i t s . pca $ sdev
p l o t data = data . frame (
e i g e n v a l u e number = 1 : l e n g t h ( e i g e n v a l u e s ) ,
eigenvalue size = eigenvalues
)
g g p l o t ( p l o t data , a e s ( x = e i g e n v a l u e number , y = e i g e n v a l u e s i z e ) ) +
geom p o i n t ( ) +
geom l i n e ( ) +
geom v l i n e ( x i n t e r c e p t = c ( 4 , 7 ) , l i n e t y p e = ” d o t t e d ” ) +
l a b s ( x = ” E i g e n v a l u e Number” , y = ” E i g e n v a l u e S i z e ” , t i t l e = ” Trend
o f P r i n c i p a l Components ” ) +
theme g r e y ( )

[Link]$sdev extracts the standard deviations of the principal components obtained from
the PCA analysis performed with the prcomp function. The plot depicted above illustrates the

progression of eigenvalues, enabling us to discern notable inflection points, commonly referred to

as ”elbows,” where the graph experiences a sharp decline. To emphasize these significant points,
vertical dotted lines have been superimposed at points 4 and 7. However, it’s equally crucial to
grasp the extent of variation attributable to these principal components in relation to the overall
variation. Hence, let us look at the percentage of variation given by, respectively, the first 3 and 6
principal components.

sum ( ( p e n d i g i t s . pca $ sdev [ 1 : 3 ] ) ˆ 2 ) /sum ( ( p e n d i g i t s . pca $ sdev ) ˆ 2 ) ;

sum ( ( p e n d i g i t s . pca $ sdev [ 1 : 6 ] ) ˆ 2 ) /sum ( ( p e n d i g i t s . pca $ sdev ) ˆ 2 )

[1] 0.6576565
[1] 0.8507928

Despite the initial eight principal components encapsulating approximately 85% of the variance
within the dataset, their inclusion fails to achieve a sufficiently substantial reduction in dimension-
ality. Consequently, in light of this inadequacy, we deliberate over our options and ultimately opt
to retain only the first three principal components. While this selection may entail sacrificing some
explanatory power, capturing only 66% of the total variance, we deem it a pragmatic compromise,
balancing the need for dimensionality reduction with the preservation of meaningful information.

17
2. Normality Principal Components
In order to have a clear overview to address multivariate normality of the first three components,
we invesitigate univariate and bivariate normality first. The code remains essentially the same as
the one done for the third point of the first exercise just substituting airquality with pc pendigits.
pc p e n d i g i t s = p e n d i g i t s . pca $ [ , 1 : 3 ]
j = 1
x = pc p e n d i g i t s [ , j ]
sorted order = order (x)
d f = data . frame ( y = s o r t ( x ) )
g g p l o t ( df , a e s ( sample = y ) ) +
s t a t qq ( ) +
s t a t qq l i n e ( l i n e t y p e = ” dashed ” ) +
l a b s ( x = ” T h e o r e t i c a l Q u a n t i l e s ” , y = ” Sample Q u a n t i l e s ” , t i t l e = ”QQ
P l o t F i r s t P r i n c i p a l Component” ) +
theme g r e y ( )

The qqplots of PC1, PC2 and PC3 raise some questions about the univariate normality, as the
tails appear to be far from the red qqline.
For inspecting bivariate normality we plot the three components against each other and we add
the 0.5 level ellipses. If the observations were normally distributed we would expect nearly 50% of
the points to be enclosed by the related [Link] the commands will be essentially the same
for each couple of variables, we firstly create two indexes variable j and k, and use them to select
the variable as well as to specify the title of the plot
j = 1
k = 2
x = pc p e n d i g i t s [ , j ]
y = pc p e n d i g i t s [ , k ]
d f = data . frame ( x , y )
cov mat = cov ( d f )
e l l i p s e data = e l l i p s e : : e l l i p s e ( cov mat , c e n t r e = c ( 0 , 0 ) , l e v e l = 0 . 5 )
e l l i p s e p o l y = S p a t i a l P o l y g o n s ( l i s t ( Polygons ( l i s t ( Polygon ( e l l i p s e data )
) , ID = ” e l l i p s e ” ) ) )
o u t s i d e p o i n t s = p o i n t . i n . polygon ( d f $x , d f $y , e l l i p s e poly@polygons
[ [ 1 ] ] @Polygons [ [ 1 ] ] @coords [ , 1 ] , e l l i p s e poly@polygons [ [ 1 ] ] @Polygons
[ [ 1 ] ] @coords [ , 2 ] )
sum ( o u t s i d e p o i n t s ) / l e n g t h ( x )
g g p l o t ( df , a e s ( x = x , y = y ) ) +
geom p o i n t ( ) +
geom path ( data = a s . data . frame ( e l l i p s e data ) , a e s ( x = x , y = y ) ,
c o l o r = ” blue ” ) +
l a b s ( x = ” F i r s t P r i n c i p a l Component” , y = ” Second P r i n c i p a l Component
” , t i t l e = ” S c a t t e r p l o t with E l l i p s e Overlay ” ) +

18
 theme g r e y ( )

[1] 0.3862809

[1] 0.4049309

[1] 0.3773654

19
The line sum(...) calculates the percentage of points inside an ellipse. It uses the Mahalanobis
distance concept. Specifically, it computes the squared Mahalanobis distance of each point (x, y)
relative to the mean of the data (which is zero, because the principal component has been normal-
ized) and checks if it falls within the chi-squared distribution’s critical value for a 50% confidence
level. We calculate the probability of points falling inside the ellipses by dividing the number of
points meeting our criterion by the total number of points. The 0.5 levels reveal some problems
as they barely reach a proportion of 40%. Therefore the normality of the components is still
debatable.
We revisit the analysis previously conducted in point 5 of the first exercise, focusing on the
first three principal components. The code remains fundamentally unchanged, with the only
modification being the replacement of airquality with pc pendigits.

d = mahalanobis ( pc p e n d i g i t s , c e n t e r = colMeans ( pc p e n d i g i t s ) , cov =

cov ( pc p e n d i g i t s ) )
d f = data . frame ( y = s o r t ( d ) )
g g p l o t ( df , a e s ( sample = y ) ) +
s t a t qq ( s i z e = 0 . 7 ) +
s t a t qq l i n e ( c o l o r = ” r e d ” , l i n e t y p e = ” dashed ” ) +
l a b s ( x = ” T h e o r e t i c a l Q u a n t i l e s ” , y = ” Sample Q u a n t i l e s ” , t i t l e = ”
Chisq Q−Q p l o t o f f i r s t t h r e e P r i n c i p a l Components ” ) +
theme g r e y ( )

20
The resulting plot reveals a deviation between sample quantiles and theoretical quantiles. Although
discrepancies are observed primarily for larger values of theoretical quantiles, indicating a departure
from strict adherence to the theoretical distribution.

3. Scatter Plots for Principal Components

We revisit the analysis previously conducted in point 4 of the first exercise, focusing on the first
three principal components. The code remains fundamentally unchanged, with only a few modifica-
tions: that are the replacement of airquality with pc pendigits and the coloring of the observations
with respect to the digit class.

lookup = c ( ” d a r k g r e e n ” , ”brown” , ” l i g h t b l u e ” , ” magenta ” , ” p u r p l e ” ,

” b l u e ” , ” r e d ” , ” l i g h t g r e e n ” , ” o r a n g e ” , ” cyan ” )
c o l . i n d = lookup [ p e n d i g i t s $ d i g i t + 1 ]
p a i r s ( pc p e n d i g i t s , pch =16 , cex =0.7 , c o l=c o l . ind , l o w e r . p a n e l=NULL)

The vector lookup comprises ten distinct color names, while [Link] is a vector whose length
matches the number of observations in the dataset, representing the colors corresponding to the
digit classes. As the values in the digit variable range from integers 0 to 9, and the indices of
lookup span from integers 1 to 10, we augment the index by 1 in lookup[pendigits$digit + 1]
to establish the association between the colors and the digit classes.

21
4. Outliers
To identify multivariate outliers, we intend to utilize the Mahalanobis distance, a metric sensitive to
the distribution’s shape and correlation structure. However, its reliable application for multivariate
outliers hinges on the assumption of normality. To validate this assumption, we recall by point 2
that even if the theoretical quantiles might not exactly match with the sample ones, this discrepancy
can be attributed to chance. Consequently, we conclude that the distribution conforms to a normal
pattern and hence we can make use of Mahalanobi’s distance.
We revisit the analysis previously conducted in point 6 of the first exercise, focusing on the first
three principal components. The code remains fundamentally unchanged, with the only modifica-
tion being the replacement of airquality with pc pendigits.

p = dim ( pc p e n d i g i t s ) [ 2 ]
d = mahalanobis ( pc p e n d i g i t s , c e n t e r=colMeans ( pc p e n d i g i t s ) , cov=cov ( pc
pendigits ) )
d f = data . frame ( x = 1 : l e n g t h ( d ) , y = d )
g g p l o t ( data = df , a e s ( x = x , y = y ) ) +
geom p o i n t ( c o l o u r = c o l . ind , shape = 1 6 ) +
geom h l i n e ( y i n t e r c e p t = q c h i s q (1 −0.02 , d f = p ) , l i n e t y p e = ” dashed ” )
+
geom h l i n e ( y i n t e r c e p t = q c h i s q (1−1/ l e n g t h ( d ) , d f = p ) , l i n e t y p e = ”
dashed ” ) +
geom h l i n e ( y i n t e r c e p t = q c h i s q (1 −0.5 / l e n g t h ( d ) , d f = p ) , l i n e t y p e = ”
dashed ” ) +
l a b s ( x = ” Index ” , y = ” Values ” ) +
theme g r e y ( )

22
In the code, we adjusted the quantile threshold from 0.05 to 0.02, focusing on outliers of greater
significance. Additionally, the large number of observations causes the value of 1/n to approach 1,
resulting in quantiles with continuity correction to be very high. As a result, these quantiles may
not provide particularly meaningful insights.
So, in summary, the following observations can be considered multivariate outliers of the first
three principal components of the dataset pendigits according to the level of 0.98.

which ( d>q c h i s q (1 −0.02 , d f = p ) )

[1] 35 63 105 248 355 357 471 488 628 634 679 694
[13] 705 1011 1116 1197 1299 1446 1460 1483 1554 1608 1635 1696
[25] 1716 1885 2113 2136 2273 2291 2373 2452 2577 2637 2640 2683
[37] 2687 2757 2937 2948 2998 3216 3271 3498 3799 3906 4041 4056
[49] 4114 4170 4210 4251 4287 4337 4346 4371 4732 5074 5331 5334
[61] 5454 5648 5908 5948 6075 6085 6151 6169 6273 6291 6302 6494
[73] 6528 6615 6833 6878 6882 6953 7257 7350 7464 7489 7555 7645
[85] 7921 8066 8291 8351 8360 8740 8902 9262 9445 9536 9781 9889
[97] 10340 10443 10881

23