• k-means clustering is a method of cluster analysis which aims to partition \(n\) observations into \(k\) clusters in which each observation belongs to the cluster with the nearest mean.
• It is similar to the expectation-maximization algorithm for mixtures of Gaussians in that they both attempt to find the centers of natural clusters in the data.

• Initialize: choose \(k\) points as cluster means

• Repeat until convergence:
  • Assignment: place each point \(X_i\) in the cluster with the closest mean.
  • Update: recalculate the mean for each cluster
• K-means always converges.
  • The assignment and update steps always either reduce the objective function or leave it unchanged.

• Advantage: gives sharp partitions of the data
• Disadvantage: need to specify the number of clusters (K).
• Goal: find a set of \(k\) clusters that minimizes the distances of each point in the cluster to the cluster mean:

\[centroid_j = \hat{Y_j} = \frac{1}{N_{Y_j}}\sum_{i \in Y_j}{X_i}\]

\[argmin_C \sum_{i=1}^k \sum_{j \in C(i)} {|X_j - \hat{Y_i}|^2}\]

• Final results depend on starting values
• How do we chose \(K\)? There are methods but not much theory saying what is best.
• Where are the pretty pictures?

• Self organizing map (SOM) is a learning method which produces low dimension data (e.g. \(2D\)) from high dimension data (\(nD\)) through the use of self-organizing neural networks
• E.g. an apple is different from a banana in more then two ways but they can be differentiated based on their size and color only.

• Dimensionality Reduction (DR) essentially aims to find low dimensional representations of data while preserving their key properties
• DR methods can be modeled after physical models with attracting and repelling forces (Force Directed Methods), projections onto low dimensional planes (PCA, ICA), divergence of statistical distributions (SNE family), the reconstruction of local spaces or points by their neighbors (LLE)

Guido Kraemer, Markus Reichstein, and Miguel D. Mahecha, “DimRed and CoRanking Unifying Dimensionality Reduction in R,” The R Journal, 2018, https://journal.r-project.org/archive/2018/RJ-2018-039/index.html.

• Linearly decompose the dataset into components that have a desired property.

• There are largely two kinds of projection methods: principal component analysis (PCA) and independent component analysis (ICA).

• PCA produces a low-dimensisonal representation of a dataset.
• Each successive principal component is selected to be orthonormal to the previous ones, and to capture the maximum information that is not already present in the previous components.
• Components are linear combinations of the original data

• PCA also serves as a tool for data visualization

• Start with many measurements (high dimensional).
• Want to reduce to few features (lower-dimensional space).
• One way is to extract features based on capturing groups of variance.
• Another could be to preferentially select some of the current features most representative of the data.

• Eigenvectors of covariance matrix.
  
• Find orthogonal groups of variance.
  
• Order from most to least variance.
  
• Components of variation.
  
• Linear combinations explaining the variance.

http://setosa.io/ev/principal-component-analysis/

• Performs a rotation of the data that maximizes the variance in the new axes
• Projects high dimensional data into a low dimensional sub-space (visualized in 2-3 dims)
• Often captures much of the total data variation in a few dimensions (< 5)
• Exact solutions require a fully determined system (matrix with full rank), i.e. a “square” matrix with independent rows

• The first principal component of a set of features \(X_1, X_2,...,X_p\) is the normalized linear combination of the features

\[Z_1 = \phi_{11}X_1 + \phi_{21}X_2 + ... + \phi_{p1}X_p\]

that has the largest variance. Note “normalized” - \(\sum_{j=1}^p\phi_{j1}^2=1\)

• The elements \(\phi_{11}, \phi_{21},... \phi_{p1}\) are the loadings of the first principal component. Together, them make up the principal component loading vector \(\phi_1 = (\phi_{11}, \phi_{21},... \phi_{p1})^T\)
• The loadings are constrained so that their sum of squares is equal to one, since othewise setting these elements to be arbitrary large in absolute value could result in an arbitrary large variance.

• Input: a \(n \times p\) data set \(X\). Since we are only interested in variance, we assume that each of the variables in \(X\) has been centered to have mean zero (that is, the column means of \(X\) are zero).
• We then look for the linear combination of the sample feature falues of the form

\[z_{i1} = \phi_{11}x_{i1} + \phi_{21}x_{i2} + ... + \phi_{p1}x_{ip}\]

for \(i=1,...,n\) that has largest sample variance under the constraint that \(\sum_{j=1}^p\phi_{j1}^2=1\)

• Since each of the \(x_{ij}\) has mean zero, so does \(z_{i1}\).
• Hence the sample variance of the \(z_{i1}\) can be written as \(\frac{1}{n}\sum_{i=1}^n z_{i1}^2\)

• Plugging in the sample variance equation the first principal component loading vector solves the optimization problem

\[maximize_{\phi_{11},...\phi_{p1}}\frac{1}{n}\sum_{i=1}^n(\sum_{j=1}^p\phi_{j1}x_{ij})^2\]

subject to \(\sum_{j=1}^p\phi_{j1}^2=1\)

• The problem can be solved via a singular value decomposition of the matrix \(X\)
• \(Z_1\) is the first principal component with values \(z_{11}, ... z_{n1}\)

• The loading vector \(\phi_1\) with elements \(\phi_{11}, \phi_{21},... \phi_{p1}\) defines a direction in feature space along which the data vary the most

• If we project the \(n\) data points \(x_1,..., x_n\) onto this direction, the projected values (the new coordinates) are the principal component scores \(z_{11}, ... z_{n1}\) themselves

• The second principal component is the linear combination of \(X_1, ...X_p\) that has maximal variance among all linear combinations that are uncorrelated with \(Z_1\)
• The second principal component scores \(z_{12}, ... z_{n2}\) take the form

\[z_{i1} = \phi_{12}x_{i1} + \phi_{22}x_{i2} + ... + \phi_{p2}x_{ip}\]

where \(\phi_2\) is the second principal component loading vector, with elements \(\phi_{12}, \phi_{22},... \phi_{p2}\)

• Constraining \(Z_2\) to be uncorrelated with \(Z_1\) is equivalent to constraining the direction \(\phi_2\) to be orthogonal to the direction \(\phi_1\). And so on for the other components

• The principal component directions \(\phi_1, \phi_2, \phi_3,...\) are the ordered sequence of right singular vectors of the matrix \(X\)
• The variances of the components are the \(\frac{1}{n}\) times the squares of the singular values

• There are at most \(min(n-1, p)\) principal components

• Given a gene-by-sample matrix \(X\) we decompose (centered and scaled) \(X\) as \(USV^T\)
• We don’t usually care about total expression level and the dynamic range which may be dependent on technical factors
• \(U\), \(V\) are orthonormal
• \(S\) diagonal-elements are eigenvalues = variance explained

• Columns of \(V\) are
  • Principle components
  • Eigengenes/metagenes that span the space of the gene transcriptional responses
• Columns of \(U\) are
  • The “loadings”, or the correlation between the column and the component
  • Eigenarrays/metaarrays - span the space of the gene transcriptional responses
• Truncating \(U\), \(V\), \(D\) to the first \(k\) dimensions gives the best \(k\)-rank approximation of \(X\)

• Kaiser’s criterion: eigenvalue \(\ge\) 1

• Cattell’s scree test: The number of components to retain corresponds to the point where the graph of eigenvalues of successive components vs. the component number begins to level off

• PCA assumes multivariate normally distributed data - gene expression data are super-Gaussian
• ICA models observations as a linear combination of latent feature variables, or components, which are chosen to be as statistically independent as possible. Think about a number of sound sources mixed together, ICA tries to “un-mix” them
• ICA is a linear rotation of the data just as PCA but instead of recovering the maximum variance, it recovers statistically independent components
• For microarray data, observations consist of microarray gene expression measurements, and independent components are interpreted to be transcriptional modules that often correspond to specific biological processes

http://www.sciencedirect.com/science/article/pii/S1532046410001000

• \(X\) - an \(m \times n\) matrix of \(n\) genes and \(m\) experiments
• ICA models this expression matrix as a linear combination of intependent biological processes by decomposing \(X\) as: \[X = AS\]
• \(S\) is a \(k \times n\) source matrix
• \(A\) is a \(m \times k\) mixing matrix
• \(k\) is a user supplied parameter \(\leq min(m,n)\)

Same preprocessing as for PCA - filter, center, scale

• \(S\) is a \(k \times n\) source matrix

• The components, or rows of \(S\), are independent in the sense that the gene weights in each component reflect samplings of independent random variables.
• In the context of gene expression, this suggests that the sets of genes comprising the groups strongly contributing to each component have independent compositions.

• Columns of \(A\) are the distribution of the component’s expression in arrays (rows of \(S\))

fastICA R package, https://cran.r-project.org/web/packages/fastICA/index.html

http://www.sciencedirect.com/science/article/pii/S1532046410001000

• The source matrix \(S\) is used to biologically interpret the components by studying their contributing genes

• The matrix \(A\) is used to associate the component with sample features by studying the distribution of the samples on the components according to their characteristics (e.g clinical or molecular variables).

• MineICA - Analysis of an ICA decomposition obtained on genomics data https://bioconductor.org/packages/release/bioc/html/MineICA.html

• \(a_{.j}=\sum_{k=1}^Kw_{.k}h_{kj},\ or\ A=WH\) (\(A\), \(W\), \(H\) are non-negative)
• \(A\) is an \(M \times N\) matrix, input gene expression by samples matrix
• \(W\) is an \(M \times K\) matrix, contains \(K\) basis vectors
• \(H\) is a \(K \times N\) matrix, contains the coefficient vectors
• \(H\) defines a meta-gene space: similar to eigengenes
• Classification can be done in the meta-gene space

• A key feature of NMF is the ability to identify nonsubtractive patterns that together explain the data as a linear combination of its basis vectors.
• \(A=WH\) (\(A\), \(W\), \(H\) are non-negative)
• \(A\) is an \(M \times N\) matrix, input gene expression by samples matrix
• \(W\) is an \(M \times K\) matrix, contains \(K\) basis vectors
• \(H\) is a \(K \times N\) matrix, contains the coefficient vectors
• The \(K\) basis vectors in \(W\) can be regarded as the ‘building blocks’ of the data
• The \(K\) coefficient vectors describe how strongly each ‘building block’ is present in the data

NMF is more than 30-year old!

• previous variants referred to as:
  • nonnegative rank fatorization (Jeter and Pye, 1981; Chen, 1984);
    
  • positive matrix factorization (Paatero and Tapper, 1994);
• popularized by Lee and Seung (1999) for “learning the parts of objects”.

Since then, widely used in various research areas for diverse applications

NMF can handle overlapping clusters and provides soft cluster membership indications.

• Many computational methods
  • Cost function \(|A-WH|\)
  • The squared Euclidean error function: \(F(W,H)=\parallel A - WH \parallel_F^2\)
  • Kullback–Leibler divergence \(D_{KL}(A \parallel \hat{A}) = \sum_{ij}(A_{ij}log\frac{A_{ij}}{\hat{A_{ij}}})\), where \(\hat{A}=WH\)
• Optimization procedure
  • Most use stochastic initialization, and the results don’t always converge to the same answer

Lee, Daniel D., and H. Sebastian Seung. “Algorithms for Non-Negative Matrix Factorization.” In Advances in Neural Information Processing Systems 13, edited by T. K. Leen, T. G. Dietterich, and V. Tresp, 556–62. MIT Press, 2001. http://papers.nips.cc/paper/1861-algorithms-for-non-negative-matrix-factorization.pdf.

• \(A=WH\) : Toy Biological interpretation
• Assume \(k=2\)
• We have 2 transcription factors that activate gene signatures \(W1\) and \(W2\)
• \(H\) represents the activity of each factor in each sample
• TF effects are additive

• NMF operates in the original non-negative measurement space
• Highly expressed genes matter more
• Positivity constraint is advantageous: positive correlation among genes is more likely to be biologically meaningful
• NMF may more accurately capture the data generating process

MDS attempts to

• Identify abstract variables which have generated the inter-object similarity measures
• Reduce the dimension of the data in a non-linear fashion
• Reproduce non-linear higher-dimensional structures on a lower-dimensional display

\[stress=\sqrt{\frac{\sum{(d_{ij} - \hat{d_{ij}})^2}}{\sum{d_{ij}^2}}}\]

• Goodness-of-fit - Measures degree of correspondence between distances among points on the MDS map and the matrix input.
• Start with distances \(d_{ij}\)
• Fit decreasing numbers \(\hat{d_{ij}}\)
• Subtract, square, sum
• Take a square root
  
• Divide by a scaling factor

• Obtain and order the \(M\) pairs of similarities
• Try a configuration in \(q\) dimensions
  • Determine inter-item distances and reference numbers
  • Minimize Kruskal’s stress
• Move the points around to obtain an improved configuration
• Repeat until minimum stress is obtained

• PCA tries to preserve the covariance of the original data

• MDS tries to preserve the metric (ordering relations) of the original space

• Stochastic Neighbor Embedding (SNE) is a technique that minimizes the Kullback-Leibler divergence of scaled similarities of the points \(i\) and \(j\) in high dimensional space, \(p_{ij}\), and low dimensional space, \(q_{ij}\)

\[KL(P||Q) = \sum_{i \ne j}p_{ij}log\frac{p_{ij}}{q_{ij}}\]

• SNE uses a Gaussian kernel (\(exp(-\frac{||x_i - x_j||^2}{2\sigma^2})\)) to compute similarities in high and low dimensional space. \(\sigma\) is a length scale parameter accounting for the width of the kernel

• Improves on SNE by using a t-Distribution as a kernel in low dimensional space
• Because of the heavy-tailed t-distribution, t-SNE maintains local neighborhoods of the data better and penalizes wrong embeddings of dissimilar points. This property makes it especially suitable to represent clustered data and complex structures in few dimensions
• t-SNE has one parameter, perplexity, to tune which determines the neighborhood size of the kernels used.

• Often t-SNE is performed on PCA components
• Liberal number of components.
• Removes mild signal (assumption of noise).
• Faster, on less data but, hopefully the same signal.

• Awesome Blog on t-SNE parameterization: http://distill.pub/2016/misread-tsne
  
• Publication: https://lvdmaaten.github.io/publications/papers/JMLR_2008.pdf
  
• Another YouTube Video: https://www.youtube.com/watch?v=RJVL80Gg3lA
  
• Code: https://lvdmaaten.github.io/tsne/
  
• Interactive Tensor flow: http://projector.tensorflow.org/