In this chapter, we turn our attention to the visualization of high-dimensional data with the aim to discover interesting patterns. We cover heatmaps, i.e., image representation of data matrices, and useful re-ordering of their rows and columns via clustering methods. To scale up visualization to very high-dimensional data, we furthermore introduce Principal Component Analysis as a dimension reduction technique.
Lower cases are used for scalars (e.g. \(x\)), bold lower cases for vectors (e.g. \(\mathbf{x}\)) and bold upper cases for matrices (e.g. \(\mathbf{X}\)). The transpose of a matrix or of a vector is denoted with a T-superscript (e.g. \(\mathbf{X}^\top\)). The Euclidean norm of vector \(\mathbf{x}\) is denoted \(||\mathbf{x}||\).
Methods of this chapter assume numeric variables. If encountered, categorical variables can be transformed to numeric variables by one-hot encoding 1.
We denote \(n\) the number of observations, \(p\) the number of variables, and \(\mathbf{X}\) the \(n \times p\) data matrix.
We use a subset of the base R mtcars dataset consisting of 10 rows (cars) and four selected variables. We store this data into the numeric matrix mat. For ease, we give full names (rather than abbreviations) to the columns and keep the row names (car names).
from plotnine.data import mtcars
mat = mtcars.iloc[1:10][
["mpg", "carb", "hp", "wt"]
].copy()
mat.index = mtcars.name[1:10]
mat.columns = [
"Miles.per.gallon",
"Carburetor",
"Horsepower",
"Weight"
]
mat.head()| Miles.per.gallon | Carburetor | Horsepower | Weight | |
|---|---|---|---|---|
| name | ||||
| Mazda RX4 Wag | 21.0 | 4 | 110 | 2.875 |
| Datsun 710 | 22.8 | 1 | 93 | 2.320 |
| Hornet 4 Drive | 21.4 | 1 | 110 | 3.215 |
| Hornet Sportabout | 18.7 | 2 | 175 | 3.440 |
| Valiant | 18.1 | 1 | 105 | 3.460 |
Beyond 5 to 10 variables, matrix scatterplots cannot be rendered with enough resolution to be useful. However, heatmaps which simply display data matrices as an image by color-coding its entries, become handy. Heatmaps allow visualization of data matrices of up to ca. 1,000 rows and columns (order of magnitude), i.e the pixel resolution of your screen (and maybe of your eyes!).
To draw heatmaps, we recommend using the library sns.clustermap (pretty heatmaps), which offers convenient functionalities in particular for clustering (See Section 5.4.2).
Here is a basic call to sns.clustermap on our data matrix mat:
import seaborn as sns
g = sns.clustermap(
mat,
row_cluster=False,
col_cluster=False,
figsize=(4,4)
)
Strikingly, the horsepower variable saturates the color scale because horsepower lives in a different scale than the other variables. Consequently, we barely see variations in the other variables.
Bringing variables to a common scale is useful for visualization but also for computational and numerical reasons. Moreover, it makes analysis independent of the units chosen. For instance the mtcars dataset provide car weights in 1,000 pounds and gas consumption in miles per gallon. We would certainly want our analysis to be the same if these variables were expressed with the metric system.
The widely used operations to bring variables to a same scale are:
centering: subtracting the mean
standard scaling or Z-score normalization: centering then dividing by the standard deviation
These operations are implemented in the base scikit-learn function StandardScaler() and are often offered as parameters of other functions. Scaling is usually done for variables (data matrix columns). However, in some application contexts, row-scaling can be considered as well.
With sns.clustermap we can also scale the data by rows or columns by setting the argument z_score accordingly. We do it here in the classical way, by column:
sns.clustermap(
mat,
row_cluster=False,
col_cluster=False,
z_score=1,
figsize=(4,4)
)
Scaling allows us to appreciate the variation for all variables. The default color scale is centered on 0. Because each column is centered, we find positive (red-ish) and negative (blue-ish) values in each column.
While all data are rendered in the Figure 5.1, it is hard to see a pattern emerging. Which cars are similar to each others? Which variables are similar to each other? Clustering is the task of grouping observations by similarities. Clustering helps finding patterns in data matrices. Clustering can also be applied to variables, by simply applying clustering methods on the transpose of the data matrix.
There are several clustering algorithms. In the following sections, we explain two widely used clustering methods: K-means clustering and hierarchical clustering.
K-Means clustering aims to partition the observations into \(K\) non-overlapping clusters. The number of clusters \(K\) is predefined. The clusters \(C_1,...C_K\) define a partition of the observations, i.e., every observation belongs to one and only one cluster. To this end, one makes use of so-called cluster centroids, denoted \(\boldsymbol{\mu_1}, ..., \boldsymbol{\mu_K}\), and associates each observation to its closest centroid (Figure 5.2).
Formally, one aims to determine the clusters \(C_1,...,C_K\) and the centroids \(\boldsymbol{\mu_1},...,\boldsymbol{\mu_K}\) in order to minimize the within-cluster sum of squares:
\[ \min_{C_1,...,C_K, \boldsymbol{\mu_1},...,\boldsymbol{\mu_k}} \sum_{k=1}^K\sum_{i \in C_k}|| \mathbf{x_i} - \boldsymbol{\mu_k}||^2 \tag{5.1}\]
where \(||\mathbf{x_i} - \boldsymbol{\mu_k}||^2 = \sum_{j=1}^p(x_{i,j} - \mu_{k,j})^2\) is the squared Euclidean distance between observation \(\mathbf{x_i}\) (the \(i\)-th row vector of the data matrix) and the centroid \(\boldsymbol{\mu_k}\).
The minimization problem (Equation 5.1) is difficult because of the combinatorial number of partitions of \(n\) observations into \(K\) clusters. However, two useful observations can be made.
First, if we assume that the positions of the centroids _k are given, then each summand \(|| \mathbf{x_i} - \boldsymbol{\mu_k}||^2\) in Equation 5.1 can be minimized by including the observation \(\mathbf{x_i}\) to the cluster of its closest centroid. The number of clusters \(K\) depends on the dataset.
Second, if we now assume that the clusters are given, then the values of the centroids that minimize \(\sum_{i \in C_k}|| \mathbf{x_i} - \boldsymbol{\mu_k}||^2\) are the observation means, i.e. \(\boldsymbol{\mu_k} = \frac{1}{|C_k|}\sum_{i \in C_k} \mathbf{x_i}\). Hence, the centroids are the cluster means, giving the name of the algorithm.
These two observations lead to an iterative algorithm that provides in practice good solutions, even though it does not guarantee to find the optimal solution.
K-Means algorithm (Figure 5.3)
Choose the \(K\) initial centroids (one for each cluster). Different methods such as sampling random observations are available for this task.
Assign each observation \(\mathbf{x_i}\) to its nearest centroid by computing the Euclidean distance between each observation to each centroid.
Update the centroids \(\boldsymbol{\mu_k}\) by taking the mean value of all of the observations assigned to each previous centroid.
Repeat steps 2 and 3 until the difference between new and former centroids is less than a previously defined threshold.
At every iteration, and at every step 2 and 3, the within-cluster sum of squares (Equation 5.1) decreases.
However, there is no guarantee to reach the optimal solution. In particular, the final clustering depends on the initialization (step 1). To not overly depend on the initialization, the K-means algorithm is typically executed with different random initializations. The clustering with lowest within-cluster sum of squares is then retained.
We have to make sure that the following assumptions are met when performing k-Means clustering (Figure 5.4, Source: scikit-learn 2):
Let us now apply K-means to the mat data matrix searching for 2 clusters. This can be easily achieved with the scikit-learn function KMeans(). While not necessary, it is a good idea to scale the variables, in order not to give the variables with larger scales too much importance (by dominating the Euclidean distances). Another way to look at it, is that the scaling reduces to some extent the problem of anisotropic clusters. In the following code, we first scale the data matrix with StandardScaler(). We also use the argument n_init of KMeans() to perform multiple random initializations.
import sklearn
scaler = sklearn.preprocessing.StandardScaler()
X = pd.DataFrame(
scaler.fit_transform(mat),
index=mat.index,
columns=mat.columns
)
k = 2
clust_km = sklearn.cluster.KMeans(n_clusters=k, n_init=20)
clust_km.fit(X)
clust_km.labels_array([1, 1, 1, 0, 0, 0, 1, 1, 0], dtype=int32)We now update our heatmap with the results of the clustering. We make use of the row_colors argument of sns.clustermap() which generates color-coded row annotations on the left side of the heatmap. We furthermore order the rows of the data matrix by cluster.
row_ann = pd.Series(
["C"+str(s) for s in clust_km.labels_],
index = mat.index
)
unique_categories = np.unique(row_ann.values)
palette = sns.color_palette("husl", len(unique_categories))
color_mapping = dict(zip(unique_categories, palette))
row_col = row_ann.map(color_mapping)
sns.clustermap(
X.iloc[np.argsort(clust_km.labels_)],
row_colors=row_col,
row_cluster=False,
col_cluster=False,
figsize=(4,4)
)
Cluster \(C_1\) appears to group the heavy, powerful and gas-consuming cars and cluster \(C_2\) the light, less powerful and more economic cars.
A major limitation of the K-means algorithm is that it relies on a predefined number of clusters. What if the interesting number of clusters is larger or smaller? Hierarchical clustering allows exploring multiple levels of clustering granularity at once by computing nested clusters. It results in a tree-based representation of the observations, called a dendrogram. Figure 5.6 shows an example of a hierarchical clustering using two variables only.
Unlike with K-means, there is no objective function associated with hierarchical clustering. Hierarchical clustering is simply defined by how it operates.
We describe bottom-up (a.k.a. agglomerative) hierarchical clustering:
Initialization: Compute all the \(n(n − 1)/2\) pairwise dissimilarities between the \(n\) observations. Treat each observation as its own cluster. A typically dissimilarity measure is the Euclidean distance. Other dissimilarities can be used (1-correlation), Manhattan distance, etc.
For \(i=n, n-1, ..., 2\):
Fuse the two clusters that are least dissimilar. The dissimilarity between these two clusters indicates the height in the dendrogram at which the fusion should be placed.
Compute the new pairwise inter-cluster dissimilarities among the \(i − 1\) remaining clusters using the linkage rule.
The linkage rules define dissimilarity between clusters. Here are four popular linkage rules:
Complete: The dissimilarity between cluster A and cluster B is the largest dissimilarity between any element of A and any element of B.
Single: The dissimilarity between cluster A and cluster B is the smallest dissimilarity between any element of A and any element of B. Single linkage can result in extended, trailing clusters in which single observations are fused one-at-a-time.
Average: The dissimilarity between cluster A and cluster B is the average dissimilarity between any element of A and any element of B.
Centroid: The dissimilarity between cluster A and cluster B is the dissimilarity between the centroids (mean vector) of A and B. Centroid linkage can result in undesirable inversions.
In Python, Hierarchical clustering can be performed using either the AgglomerativeClustering class from scikit-learn, or the linkage function from scipy. In both cases, the clustering entails two steps. First, it computes the distance between observations (rows of a DataFrame) across variables (columns of a DataFrame). Here, the Euclidean distance between rows is computed by default. Alternatives include the Manhattan or Minkowski distance. As for K-means, it is recommended to work on scaled variables to not give too much importance to variables with large variance. We therefore compute the pairwise Euclidean distance to our scaled data matrix X. Second, we use the resulting Euclidean distance matrix as a dissimilarity matrix to perform Hierarchical clustering. The default linkage method is single for linkage, and ward 3 for Agglomerative Clustering. Here we use average.
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
hc_scipy = linkage(X, method="average")
hc_sklearn = sklearn.cluster.AgglomerativeClustering(
linkage="average"
).fit_predict(X)The results of hierarchical clustering can be shown using a dendrogram (i.e., a tree representation). Here, observations that are determined to be similar by the clustering algorithm are displayed close to each other in the x-axis. The height in the dendrogram at which two clusters are merged represents the distance between those two clusters.
print(dendrogram(hc_scipy)){'icoord': [[25.0, 25.0, 35.0, 35.0], [55.0, 55.0, 65.0, 65.0], [45.0, 45.0, 60.0, 60.0], [75.0, 75.0, 85.0, 85.0], [52.5, 52.5, 80.0, 80.0], [30.0, 30.0, 66.25, 66.25], [15.0, 15.0, 48.125, 48.125], [5.0, 5.0, 31.5625, 31.5625]], 'dcoord': [[0.0, np.float64(1.686236409517311), np.float64(1.686236409517311), 0.0], [0.0, np.float64(0.8551281831143587), np.float64(0.8551281831143587), 0.0], [0.0, np.float64(1.3039349663134179), np.float64(1.3039349663134179), np.float64(0.8551281831143587)], [0.0, np.float64(1.5981431443987575), np.float64(1.5981431443987575), 0.0], [np.float64(1.3039349663134179), np.float64(2.1765166664822457), np.float64(2.1765166664822457), np.float64(1.5981431443987575)], [np.float64(1.686236409517311), np.float64(2.457718283406245), np.float64(2.457718283406245), np.float64(2.1765166664822457)], [0.0, np.float64(3.140614226663734), np.float64(3.140614226663734), np.float64(2.457718283406245)], [0.0, np.float64(4.212933000617484), np.float64(4.212933000617484), np.float64(3.140614226663734)]], 'ivl': ['5', '1', '0', '8', '2', '6', '7', '3', '4'], 'leaves': [5, 1, 0, 8, 2, 6, 7, 3, 4], 'color_list': ['C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C0', 'C0'], 'leaves_color_list': ['C0', 'C0', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1', 'C1']}
As illustrated before, the library sns.clustermap enables the easy creation of heatmaps. In the previous example, we set the parameters row_cluster and col_cluster to False to avoid the default computation of hierarchical clustering. If we want to include hierarchical clustering, we can simply set these parameters to True or let Python consider its default values.
sns.clustermap(X, row_cluster=True, col_cluster=False, figsize=(4,4))
Compared to K-means (Figure 5.5), the hierarchical clustering shows useful different degrees of granularity of the clusters. We see the most similar cars grouping together (the Mercedes 230 and the Mercedes 240D). At the high level, the Duster 360 stands out as an outlier.
After having inspected the result of a hierarchical clustering, it is often interesting to define distinct clusters by cutting the dendrogram at a certain height.
Typically, one cuts dendrogram either at a given height, or in order to obtain a certain number of clusters. The functions fcluster from scipy and AgglomerativeClustering from scikit-learn support both options. Here is an example of cutting the dendrogram to get 3 clusters.
hc_sklearn = sklearn.cluster.AgglomerativeClustering(
n_clusters=3,
linkage="average").fit(X)
hc_sklearn.labels_array([0, 2, 0, 0, 0, 1, 0, 0, 0])fcluster(hc_scipy, t=3, criterion="maxclust")array([1, 2, 1, 1, 1, 3, 1, 1, 1], dtype=int32)Both K-means and hierarchical clustering are well established and widely used. Here, we briefly state a few differences that may be considered when deciding which algorithm to apply in practice.
The time complexity of K-Means clustering is linear, while that of hierarchical clustering is quadratic. This implies that hierarchical clustering can not handle extremely large datasets as efficiently as K-Means clustering.
In K-Means clustering, we start with a random choice of centroids for each cluster. Hence, the results produced by the algorithm depend on the initialization. Therefore, the results might differ when running the algorithm multiple times. Hierarchical clustering outputs reproducible results.
Another difference is that K-Means clustering requires the number of clusters a priori. In contrast, the number of clusters we find appropriate in hierarchical clustering can be decided a posteriori by interpreting the dendrogram.
Let us first visualize the outcome of K-means and of the hierarchical clustering cut for 3 clusters thanks to the row annotation option of sns.clustermap.
row_ann = pd.DataFrame(
{
"kmeans":["C"+str(s) for s in clust_km.labels_],
"hc":["C"+str(s) for s in hc_sklearn.labels_]
},
index = X.index
)
unique_categories = np.unique(row_ann.values)
palette = sns.color_palette("husl", len(unique_categories))
color_mapping = dict(zip(unique_categories, palette))
row_col = row_ann.map(lambda x: color_mapping[x])
sns.clustermap(
X,
row_colors=row_col,
col_cluster=False,
figsize=(5,5)
)
Comparing clustering results is a challenging task. When clustering into two groups, we could use evaluation measures from classification, which we will introduce later. Moving from two partitions of the data into arbitrarily many groups requires new ideas.
We remark that a partition is, in our context, the result from a clustering algorithm and, therefore, the divided dataset into clusters. Generally, two partitions (from different clustering algorithms) are considered to be similar when many pairs of points are grouped together in both partitions.
The Rand index is a measure of the similarity between two partitions. Formally, we introduce the following definitions:
\(S = \{o_1, \dots, o_n\}\) a set of \(n\) elements (or observations)
First partition \(X = \{X_1, \dots, X_k\}\) of \(S\) into \(k\) sets
Second partition \(Y = \{Y_1, \dots, Y_l\}\) of \(S\) into \(l\) sets
\(a\) number of pairs of elements of \(S\) that are in the same set in \(X\) and in \(Y\)
\(b\) number of pairs of elements of \(S\) that are in different sets in \(X\) and in \(Y\)
\({n}\choose{2}\) total number of pairs of elements of \(S\)
Then, the Rand index can be computed as \[ R = \frac{a + b}{ {n}\choose 2 } \]
where \({n}\choose{k}\) (reads “n choose 2”), the binomial coefficient for \(k=2\), is the number of pairs of observations and is equal to:
\[ {{n}\choose{2}} = \frac{(n-1) \cdot n }{2}. \]
By definition, the Rand index has values between 0 and 1, including them. A Rand index of 1 means that all pairs that are in the same cluster in the partition \(X\) are also in the same cluster in the partition \(Y\) and all pairs that are not in the same cluster in \(X\) are also not in the same cluster in \(Y\). Hence, the two partitions are identical with a Rand index of 1. In general, the higher the Rand index, the more similar are both partitions.
We can compute the Rand index between our k-means result and the cut of the hierarchical clustering for a given number of groups. See exercise sheet.
A heatmap is a visualization method of choice for data matrices as long as rows and columns are visually resolved because it satisfies the two main data visualization principles, i.e.,: i. having a high data/ink ratio and ii. showing the data as raw as possible.
However, beyond dimensions exceeding the thousands of variables, dimension reduction techniques are needed. The idea of dimension reduction is simple: if the dimension of our data \(p\) is too large, let us consider instead a representation of lower dimension \(q\) which retains much of the information of the dataset.
For Principal Component Analysis (Pearson, 1901), this representation is the projection of the data on the subspace of dimension \(q\) that is closest to the data according to the sums of the squared Euclidean distances.
Principal Component Analysis is not only the mother of all data reduction techniques but also still widely used. PCA enjoys several noticeable statistical properties which we will now look at.
To get an intuition of PCA, let us consider reducing a 2-dimensional dataset into a single dimension. This application has no visualization purposes but could help defining a linear combination of two variables into a aggregated score. For example, one can want to summarize the weight and horsepower of the cars into a single score.
Geometrically, we search for a line lying as close as possible to the data, in the sense of least squared Euclidean distances.
We will assume all variables to be centered and admit the line passes therefore through the origin. Let us denote \(\mathbf{w}\) a direction vector of the line of length 1.
The closest point of the observation vector \(\mathbf{x_i}\) to the line is its orthogonal projection \(p_{\top}(\mathbf{x_i})\) which is equal to the scalar product of the direction vector and the observation vector, times the direction vector:
\[ \begin{align} p_{\top}(\mathbf{x_i}) &= (\mathbf{w}^\top\mathbf{x_i})\mathbf{w} \end{align} \]
Hence, we look for \(\mathbf{w}\) such that:
\[ \begin{align} \min_{\mathbf{w}} & \sum_{i=1}^n || \mathbf{x} - (\mathbf{w}^\top\mathbf{x_i})\mathbf{w} ||^2 \\ \text{subject to} & ||\mathbf{w}||=1 \end{align} \tag{5.2}\]
There is typically a unique solution to this optimization problem (up to a sign). This direction vector \(\mathbf{w}\) is called the first principal component (PC1) of the data.
We defined PC1 with a minimization problem (Equation 5.2). One can also see it as a maximization problem. To this end, consider the orthogonal triangle defined by an observation vector \(\mathbf{x_i}\), its projection \(p_{\top}(\mathbf{x_i})\), and the origin. Pythagoras’ theorem implies that:
\[ \begin{align} ||\mathbf{x_i}||^2 = ||p_{\top}(\mathbf{x_i})||^2 + || \mathbf{x} - p_{\top}(\mathbf{x_i}) ||^2 \end{align} \]
Over the entire dataset, the sum of the \(||\mathbf{x_i}||^2\) is constant independently of the choice of \(\mathbf{w}\). Therefore minimization problem (Equation 5.2) is equivalent to:
Over the entire dataset, the sum of the \(||\mathbf{x_i}||^2\) is constant independently of the choice of \(\mathbf{w}\). Therefore minimization problem Equation (Equation 5.2) is equivalent to:
\[ \begin{align} \max_{\mathbf{w}} & \sum_{i=1}^n ||p_{\top}(\mathbf{x_i})||^2 \\ \text{subject to} & ||\mathbf{w}||=1 \end{align} \]
As we have centered the data, the origin is the mean of the data. Hence the sum of squared norms of the observation vectors is \(n\) times their total variance. By linearity, the origin is also the mean of the projected data and thus:
\[ \sum_i||p_{\top}(\mathbf{x_i})||^2 = n \mathrm{Var}(p_{\top}(\mathbf{X})) \] Hence, one can equivalently consider that PC1 maximizes the variance of the projected data.
Result PC1 maximizes the variance of the projected data.
The proportion of variance captured by PC1 is defined as the ratio of the variance of the projected data over the total variance of the data. It is a proportion, hence lies between 0 and 1. The higher it is, the smaller the sum of squared distances, the closer the line is to the data. The proportion of variance hence quantifies how good our dimension reduction is.
In the general case, with \(p\) variables and \(n\) observations, one searches for the \(q\)-dimensional plane that is closest to the data in terms of sums of squared Euclidean distances. This is also the \(q\)-dimensional plane that maximizes the variance of the projected data.
An important property relates principal components to the eigendecomposition of the covariance matrix.
The covariance matrix is \(\frac{1}{n}\mathbf{X}^{\top}\mathbf{X}\). It is a symmetric positive matrix. We denote \(\mathbf{w_1},...,\mathbf{w_j},...\) its eigenvectors ordered by decreasing eigenvalues \(\lambda_1 >...> \lambda_j>...\).
Result. The PCA \(q\)-dimensional plane, i.e., the \(q\)-dimensional plane that is closest to the data in terms of sums of squared Euclidean distances, is the plane spanned by the first \(q\) eigenvectors of the covariance matrix.
Result. The variance explained by the PCA \(q\)-dimensional plane equals to the sum of the \(q\) first eigenvalues of the covariance matrix.
See Bishop (2007) for proofs.
These results have several implications:
The PCA planes are nested: the PCA 2D-plane contains PC1, the PCA 3D-plane contains the PCA 2D-plane, etc.
We call second principal component (PC2) the second eigenvector of the covariance matrix, etc.
The principal components are linearly uncorrelated. This is because the eigenvectors of a positive matrix are orthogonal to each other. If \(n>p\) (more observation than variables) the PCs form an orthonormal basis.
A 3D illustration of PCA is available at [https://setosa.io/ev/principal-component-analysis/].
PCA can be easily performed in Python by using the scikit-learn function PCA().
In most applications, scaling the data beforehand is important. Because PCA is based on minimizing squared Euclidean distances, scaling allows to not give too much importance to variables living on larger scales than the other ones. Be careful: PCA always centers data to 0, but it does NOT scale the variables to unit variance.
In the following examples, we perform PCA on our mat dataset. We run StandardScaler to scale the data before the PCA.
X = sklearn.preprocessing.StandardScaler().fit_transform(mat)
pca_res = sklearn.decomposition.PCA().fit(X) The output can be stored in pca_res, which contains information about the values of each sample in terms of the principal components (pca.transform(X)), variance (pca.explained_variance_) of each principal component, and the relationship between the initial variables and the principal components (i.e. loadings, pca.components_).
An overview of the PCA result can be obtained from the object attributes. We remark that the cumulative proportion is always equal to one for the last principal component.
sum_df = pd.DataFrame(
{
"Variance": pca_res.explained_variance_,
"Proportion of Variance": pca_res.explained_variance_ratio_,
"Cumulative Proportion": np.cumsum(pca_res.explained_variance_ratio_)
},
index=[f"PC{i+1}" for i in range(pca_res.n_components_)]
)
sum_df| Variance | Proportion of Variance | Cumulative Proportion | |
|---|---|---|---|
| PC1 | 3.003031 | 0.667340 | 0.667340 |
| PC2 | 0.819505 | 0.182112 | 0.849453 |
| PC3 | 0.574486 | 0.127664 | 0.977116 |
| PC4 | 0.102977 | 0.022884 | 1.000000 |
In this example, the first principal component explains 66.7% of the total variance and the second one 18.2%. So, just PC1 and PC2 can explain approximately 84.9% of the total variance.
Plotting the results of PCA is particularly important. The so-called scree plot is a good first step for visualizing the PCA output, since it may be used as a diagnostic tool to check whether the PCA worked well on the selected dataset or not.
(
ggplot(
sum_df,
aes("sum_df.index", "Variance", group=1)
)
+ geom_line()
+ geom_point()
+ mytheme
+ theme(figure_size=(4,3))
)
The scree plot shows the variance in each projected direction. The y-axis contains the eigenvalues, which essentially stand for the amount of variation. We can use a scree plot to select the principal components to keep. If the scree plot has an ‘elbow’ shape, it can be used to decide how many principal components to use for further analysis. For example, we may achieve dimensionality reduction by transforming the original four dimensional data (first four variables of mtcars) to a two-dimensional space by using the first two principal components.
A variant of the scree plot can be considered by plotting the proportion of the total variance for every principal component.
(
ggplot(
sum_df,
aes("sum_df.index", "Proportion of Variance", group=1)
)
+ geom_line()
+ geom_point()
+ mytheme
+ theme(figure_size=(4,3))
)
We can access the projection of the original data on the principal components by using the method .transform(). A scatterplot shows the projection of the data on the first two principal components.
scores_df = pd.DataFrame(
pca_res.transform(X),
columns = [
f"PC{i+1}" for i in range(pca_res.n_components_)
],
index=mat.index
)
p = (
ggplot(
scores_df,
aes("PC1", "PC2")
)
+ geom_point()
+ geom_text(
aes(label="scores_df.index"),
adjust_text=adjust_text_dict
)
+ mytheme
)
p
The biplot shows the projection of the data on the first two principal components. It includes both the position of each sample in terms of PC1 and PC2 and also shows how the initial variables map onto this. The sign of the variable loadings in the new dimensional space indicates their contribution to the PCs. Here Miles per Gallon correlates positively with PC1, while the other three variables contribute negatively to the first component. Of note, concluding about the similarity of high-dimensional variables or observations from their mere vicinity on a 2D representation (be PCA or a non-linear mapping such as UMAP or t-SNE) is a common mistake in data science. The PCA projection is blind to variations that are orthogonal to the projection plane. In particular, if two variables (or two observations) are highly correlated then their projection on the plane will be correlated (small angle). However, the converse is not true. The projections of two variables (or two observations) can be strongly correlated (small angle) even though the underlying variables (or observations) are not correlated in the complete space. Correlation of the projections is a necessary but not a sufficent condition for correlation.
scale_x = scores_df["PC1"].std() * 3 # Scale loadings to fit the space
scale_y = scores_df["PC2"].std() * 3
load_df = pd.DataFrame(
{
"xstart": 0.0,
"ystart": 0.0,
"label": mat.columns,
"xend": pca_res.components_.T[:,0] * scale_x,
"yend": pca_res.components_.T[:,1] * scale_y
}
)
(
p
+ geom_segment(
load_df,
aes("xstart", "ystart", xend="xend", yend="yend"),
color="tomato",
arrow=arrow(length=0.1)
)
+ geom_text(
load_df,
aes("xend", "yend", label="label"),
color="tomato",
adjust_text=adjust_text_dict
)
)
PCA is a statistical procedure that uses an orthogonal transformation to convert a set of possibly correlated variables into a set of linearly uncorrelated variables, which are denoted as principal components.
Each principal component explains a fraction of the total variation in the dataset. The first principal component has the largest possible variance. Respectively, the second principal component has the second-largest possible variance.
In this manner, PCA aims to reduce the number of variables, while preserving as much information from the original dataset as possible. High-dimensional data is often visualized by plotting the first two principal components after performing PCA.
One limitation of PCA is that it is restricted to linear transformation of the data. What if the data lies closer to a parabola rather than a straight line? There are many non-linear alternatives to PCA including Independent Component Analysis, kernel PCA, t-SNE, UMAP. Details of these techniques are beyond the scope of this lecture. However, as long as one uses these techniques as visualization and exploratory tools rather than for making any claim on the data, a profound understanding of their theory is not necessary.
We illustrate here with one example of a PCA and a UMAP representation of same single-cell gene expression matrix of mouse. The input matrix has 15,604 rows (single cells) and 1,951 columns (genes) and comes for the Tabula muris project 4.
PCA on mouse single-cell transcriptome data. Source: Laura Martens, TUM
Clustering and dimension reduction techniques belong to the family of unsupervised learning methods. Unlike supervised learning methods (e.g. regression, classification) unsupervised learning methods shall discover patterns in the data without being guided by some ground truth.
There is no “right” clustering, or “right” subspace in real-life datasets
Clustering and dimension reduction techniques are exploratory tools meant to help deriving some hypotheses
These hypotheses are then best tested on independent data
By now, you should be able to:
G. James, D. Witten, T. Hastie and R. Tibshirani. An Introduction to Statistical Learning with Applications in R. Book and R code available at: https://www.statlearning.com/
Advanced (PCA proofs): C. Bishop, Pattern Recognition and Machine Learning. https://www.microsoft.com/en-us/research/people/cmbishop/prml-book/
