Quarto - Data 311: Machine Learning

Introduction

So far we’ve seen a few methods for clustering: $k$ -means, hierarchical.
Recall that clustering refers to a very broad set of techniques for finding subgroups, or clusters, in a data set.
This lecture will talk about Gaussian Mixture Models (GMM):

Recall Unsupervised Learning

Recall that clustering is an unsupervised task; we do not have any labelled data (i.e. we don’t have a $Y$ )
It is often easier to obtain unlabled date than labeled data (which can require human intervention)
The unsupervised task is inherently more challenging than the supervised task.

Clustering methods

There are various clustering techniques, each with its own strengths, assumptions, and use cases.
For heirachical clustering and K-means, we relied on Euclidean distance.
In contrast to these distance-based algorithns, GMM is a probabilistic model-based clustering technique that assumes that the data is generated from a mixture of multiple Gaussian distributions.

Old Faithful Data

The faithful dataset in R contains the measurements on time to eruption, and length of eruption time of the Old Faithful geyser in Yellowstone National Park.
The data frame consists of 272 observations on 2 variables:
- eruptions Eruption time in mins (numeric)
- waiting Waiting time to next eruption in mins (numeric)

N.B. to see the help file for any dataset (or function in R) use ?name_of_thing

head(faithful)

ABCDEFGHIJ0123456789

	eruptions <dbl>	waiting <dbl>
1	3.600	79
2	1.800	54
3	3.333	74
4	2.283	62
5	4.533	85
6	2.883	55

jump to the iris data set.

plot(faithful, 
     ylab = "Waiting time to next eruption (in mins)", 
     xlab = "Eruption duration (mins)")

Determining Clusters

Rather than creating the clusters ourselves, we look for automated ways for computers to find these natural groupings for us.
While this is conveneient in lower dimensional space, it is absolutely crucial in higher dimensions in which we can no longer visualize our data.
For instance we could run something like mclust¹ to obtain a optimal* partition of the data.

faithful clustered

library(mclust)
mod2 <- Mclust(faithful)
plot(mod2, what = "classification")

The points are coloured by group allocation.

Interestingly, mclust determines that there are 3 groups in the data.

But how did mclust determine that 3 clusters was optimal?

Resources

[MP] Mixture Model-Based Classification by Paul D. McNicholas

[FS] Finite Mixture and Markov Switching Models by Fruhwirth-Schnatter, S. (2006)

[MG] Finite Mixture Models by McLachlan, Geoffrey J; Peel, David/ New York: John Wiley & Sons

[LB] Mixture models: theory, geometry, and applications by Lindsay, Bruce G

Motivating Example

As in the previous units, by way of example we will look at some labelled data and then perform cluster (pretending that we don’t know the labels)
Common “benchmark” data sets in the clustering world include: iris, crabs, wine
This exercise is quite common when trying to assess a clustering algorithm (remember in the genuine cluster problem we wouldn’t calculate a confusion matrix and accuracy measure for example)

Iris

Description

Here i am using “clusters” as groups of data that clump together, but really these are just know classifcation corresponding to the species are Iris setosa, versicolor, and virginica.

Data frame

iris

ABCDEFGHIJ0123456789

Sepal.Length <dbl>	Sepal.Width <dbl>	Petal.Length <dbl>	Petal.Width <dbl>	Species <fct>
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa
4.6	3.4	1.4	0.3	setosa
5.0	3.4	1.5	0.2	setosa
4.4	2.9	1.4	0.2	setosa
4.9	3.1	1.5	0.1	setosa

jump to the geyser data set.

Plot

Code

library(ggplot2)
library(GGally)

ggpairs(iris,                 # Data frame
        columns = 1:4,        # Columns
        aes(color = Species,  # Color by group (cat. variable)
            alpha = 0.5))     # Transparency

We can see that these species form distinct clusters (esp in the Sepal.Width vs Petal.Width scatterplot)

Fitted GMM

Here we have fit a multivariate normal distribution to each sub-group in the data.

That is we find the maximum likelihood estimates (MLEs) ${\hat{μ}}_{g}$ , ${\hat{Σ}}_{g}$ for $g = 1, 2, 3$

Mixture Models

Mixture models are just a convex combination of statistical distributions. In general, for a random vector $X$ , a mixture distribution is defined by $f (x) = \sum_{g = 1}^{G} π_{g} p_{g} (x ∣ θ_{g})$ where $G$ is the number of groups, $0 \leq π_{g} \leq 1$ and $\sum_{g = 1}^{G} π_{g} = 1$ and $p_{g} (x ∣ θ_{g})$ are statistical distributions with parameters $θ_{g}$ .

Terminology

The $π_{g}$ are called mixing proportions and we will write $π = (π_{1}, π_{2}, \dots, π_{G})$
The $p_{g} (x)$ are called component densities
Since the $p_{g} (x)$ are densities, $f (x)$ describes a density called a $g$ -component finite mixture density.

Statistical Distributions

In theory, $p_{g} (x ∣ θ_{g})$ could be a different distribution type for each $g$ .
In practice, however, $p_{g} (x ∣ θ_{g})$ is assumed to be the same across groups, but with different $θ_{g}$ .
Most commonly, we assume a a Gaussian mixture model (GMM) where $p_{g} (x ∣ θ_{g})$ is assumed to be multivariate Gaussian (aka multivariate Normal), denoted $ϕ (x | μ_{g}, Σ_{g})$

Gaussian Mixture model {#gmm}

Thus $f (x)$ — sometimes written $f (x | Θ)$ or $f (x; Θ)$ — is a Gaussian mixture model would have the form: $\begin{aligned} f (x | Θ) & = \sum_{g = 1}^{G} π_{g} ϕ (x ∣ μ_{g}, Σ_{g}) \\ = \sum_{g = 1}^{G} π_{g} \frac{| Σ |^{- 1 / 2}}{(2 π)^{k / 2}} \exp (- \frac{1}{2} (x - μ_{g})^{T} Σ^{- 1} (x - μ_{g})) \end{aligned}$ where $ϕ$ denote the normal distribution, $μ_{g}$ is the mean vector for group $g$ and $Σ_{g}$ is the covariance matrix for group $g$ , and $Θ = {μ_{1}, . . ., μ_{G}, Σ_{1}, . . ., Σ_{G}, π_{1}, . . ., π_{G}}$

Distribution Alternatives

Other robust alternatives include:

Mixtures of Student- $t$ distributions (heavier tails): (Peel & McLachlan, 2000), (Andrews et. al, 2011)
Mixtures of skewed distributions
- skew-normal (Lin et al, 2007)
- skewed- $t$ (Vrbik and McNicholas, 2012)

faithful finite mixture 1D

hist(faithful$waiting, freq = FALSE)
lines(density(faithful$waiting))

This can be viewed as a mixture of 2 univariate normals

Two Univariate Normals I

curve(dnorm(x, mean = 2, sd = 1), col = 2, xlim = c(-2,9),
      lwd = 2, ylab="Density")
curve(dnorm(x, mean = 5, sd = 1), add = TRUE, col = 4, lwd = 2)

Mixture of 2 univariate Normals I

curve(0.5*dnorm(x, mean = 2, sd = 1) + 0.5*dnorm(x, mean = 5, sd = 1), 
      xlim = c(-2,9),
      lwd = 2, ylab="Density")

Mixture of 2 univariate Normals II

Code

par(mfrow=c(1,2))
mu1=4; mu2 = 2
sd1 = sd2 = 1
meanvec = c(mu1, mu2)
sdvec = c(sd1, sd2)
lind = which.min(meanvec)
rind = which.max(meanvec)
xrange = c(meanvec[lind]-3.5*sdvec[lind],
           meanvec[rind]+3.5*sdvec[rind])

# Two-curve plot
curve(dnorm(x, mean = mu1, sd = sd1), col = 2, xlim = xrange,
      lwd = 2, ylab="Density")
curve(dnorm(x, mean = mu2, sd = sd2), add = TRUE, col = 4, lwd = 2)

# One mixture curve plot
curve(0.5*dnorm(x, mean = mu1, sd = sd1) + 0.5*dnorm(x, mean = mu2, sd = sd2), 
      xlim = xrange,
      lwd = 2, ylab="Density")

Mixture of 2 univariate Normals II

Code

par(mfrow=c(1,2))
mu1=3; mu2 = 2
sd1 = sd2 = 1
meanvec = c(mu1, mu2)
sdvec = c(sd1, sd2)
lind = which.min(meanvec)
rind = which.max(meanvec)
xrange = c(meanvec[lind]-3.5*sdvec[lind],
           meanvec[rind]+3.5*sdvec[rind])

# Two-curve plot
curve(dnorm(x, mean = mu1, sd = sd1), col = 2, xlim = xrange,
      lwd = 2, ylab="Density")
curve(dnorm(x, mean = mu2, sd = sd2), add = TRUE, col = 4, lwd = 2)

# One mixture curve plot
curve(0.5*dnorm(x, mean = mu1, sd = sd1) + 0.5*dnorm(x, mean = mu2, sd = sd2), 
      xlim = xrange,
      lwd = 2, ylab="Density")

Mixture of 2 univariate Normals III

Code

par(mfrow=c(1,2))
mu1=5; mu2 = 2
sd1 = 1;  sd2 = 2
meanvec = c(mu1, mu2)
sdvec = c(sd1, sd2)
lind = which.min(meanvec)
rind = which.max(meanvec)
xrange = c(meanvec[lind]-4*sdvec[lind],
           meanvec[rind]+4*sdvec[rind])

# Two-curve plot
curve(dnorm(x, mean = mu1, sd = sd1), col = 2, xlim = xrange,
      lwd = 2, ylab="Density")
curve(dnorm(x, mean = mu2, sd = sd2), add = TRUE, col = 4, lwd = 2)

# One mixture curve plot
curve(0.5*dnorm(x, mean = mu1, sd = sd1) + 0.5*dnorm(x, mean = mu2, sd = sd2), 
      xlim = xrange,
      lwd = 2, ylab="Density")

`faithful` GMM 1D fit

Code

mod1 = Mclust(faithful$waiting)
theta = summary(mod1, parameters = TRUE)

hist(faithful$waiting, freq = FALSE, ylim = c(0,0.045), main = "MCLUST fit", xlab = "waiting time (in minutes)")
#lines(density(faithful$waiting)) 
curve(theta$pro[1]*dnorm(x, theta$mean[1], sqrt(theta$variance[1])), col = 2, add = TRUE, lwd=4) 
curve(theta$pro[2]*dnorm(x, theta$mean[2], sqrt(theta$variance[2])), col = 4, add = TRUE, lwd=4) 
curve(theta$pro[1]*dnorm(x, theta$mean[1], sqrt(theta$variance[1])) + theta$pro[2]*dnorm(x, theta$mean[2], sqrt(theta$variance[2])), col = 1, add = TRUE, lty=2, lwd = 2)

Mixture Models for clustering

Finite mixture models are ideally suited for clustering
Clusters correspond to the $1, . . ., g, . . . G$ sub-populations/groups/components in a finite mixtures model.
To derive the MLEs of the Gaussian mixture model, we first need a likelihood function …

The Likelihood

The likelihood function for this mixture is $\begin{aligned} L (Θ) & = f (x | Θ) = Π_{i} [\sum_{g = 1}^{G} π_{g} ϕ (x ∣ μ_{g}, Σ_{g})] \end{aligned}$ and the so-called log-likelihood is: $\begin{aligned} ℓ (Θ) & = \log (L (Θ)) = \sum_{i = 1}^{N} \log [\sum_{g = 1}^{G} π_{g} ϕ (x ∣ μ_{g}, Σ_{g})] \end{aligned}$

EM algorithm

Most commonly, we use the expectation-maximization (EM) algorithm, which helps us find parameter estimates with missing data (group memberships are missing).
More generally, the EM is an iterative algorithm for finding maximum likelihood parameters in statistical models, where the model depends on unobserved latent variables.
We assume that each observed data point has a corresponding unobserved latent variable specifying the mixture component to which each data point belongs.

Latent Variable

We introduce a new latent¹ variable that represents the group memberships: $Z_{i g} = {\begin{cases} 1 & if x_{i} belongs to group g \\ 0 & otherwise \end{cases}$

We call ${X, Z}$ our “complete” data.

Now that we have that, let’s look at a non-technical summary of the EM algorithm…

Latent Variable Representation

We can now represent our mixture model $\begin{aligned} f (x | Θ) & = \sum_{g = 1}^{G} π_{g} ϕ (x ∣ μ_{g}, Σ_{g}) \\ = \sum_{g = 1}^{G} p (Z_{i g} = 1) p (x ∣ Z_{i g} = 1) \end{aligned}$ This is sometimes called the latent variable representation of the mixture model.

Complete Likelihood

The “complete” likelihood takes the form: $L_{c} (Θ) = \prod_{i = 1}^{n} \prod_{g = 1}^{G} π_{g}^{{\hat{z}}_{i g}} ϕ (x_{i} | μ_{g}, σ_{g})^{{\hat{z}}_{i g}}$ so the complete log-likelihood, $\log (L_{c} (Θ))$ , takes the form: $\begin{aligned} ℓ_{c} (Θ) & = \sum_{i = 1}^{n} \sum_{g = 1}^{G} {\hat{z}}_{i g} [\log (π_{g}) + \log (ϕ (x_{i} | μ_{g}, σ_{g}))] \end{aligned}$

We’ll work with univariate normal distribution for simplicity.

EM Pseudocode

Start the algorithm with random values for ${\hat{z}}_{i g}$ .
Mstep: Assuming those ${\hat{z}}_{i g}$ are correct, find the maximum likelihood estimates (MLE) of $μ_{g}$ and $Σ_{g}$
Estep: Assuming those parameters are correct, find the expected value of group memberships ${\hat{z}}_{i g} = \frac{π_{g} ϕ (x ∣ μ_{g}, Σ_{g})}{\sum_{g = 1}^{G} π_{g} ϕ (x ∣ μ_{g}, Σ_{g})}$
Repeat 2. and 3. until “changes” are minimal.

Fuzzy z’s

The resulting ${\hat{z}}_{i g}$ are not going to be 0’s and 1’s when the algorithm converges.
For each observation $x_{i}$ we end up with vector of probabilities ( ${\hat{z}}_{i 1}$ ,…, ${\hat{z}}_{i G}$ ) which we denote by ${\hat{z}}_{i}$
So, for example, observation $x_{i}$ might have a corresponding estimate of ${\hat{z}}_{i} = (0.25, 0.60, 0.15)$ –> in this case $x_{i}$ would get assigned to group 2
But! We usually still want a “best” guess as to which group the observation belongs to. Continuing with example above what do you think might happen …

MCLUST example

To see this in action, lets take a closer look at the MCLUST output when it was fit to our old faithful geyser data.

To recap:

MCLUST imposes the GMM structure.
We introduce indicator $z_{i g}$ such that $z_{i g} = 1$ is observation $i$ belongs to group $g$ and $z_{i g} = 0$ otherwise.
The model is fitted using the EM alogrithm

Classification vector

Our estimates of $z_{i g}$ , denoted ${\hat{z}}_{i g}$ determine our clusters: $MAP {\hat{z}}_{i g} = {\begin{cases} 1 & if max_{g} {{\hat{z}}_{i g}} \\ 0 & otherwise \end{cases}$

# fuzzy z matrix
round(head(mod1$z) , 5)

        [,1]    [,2]
[1,] 0.00010 0.99990
[2,] 0.99991 0.00009
[3,] 0.00415 0.99585
[4,] 0.96757 0.03243
[5,] 0.00000 1.00000
[6,] 0.99981 0.00019

# "hard" zs = MAP(zhat)
head(cbind(ifelse(mod1$classification==1,1,0), ifelse(mod1$classification==2,1,0)))

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

# classification vector
head(mod1$classification)

[1] 2 1 2 1 2 1

Soft Clustering

GMM is a type of soft clustering.
This means that a data object can exist in more than one cluster with a certain probability or degree of membership.
This is unlike hard clustering methods which define a non-overlapping partition of the data set (eg. kmeans)
The matrix of ${\hat{z}}_{i g}$ is often referred to has a “fuzzy” z matrix.
The rows of this matrix should necessarily sum to 0.

Problem 1

EM algorithm is a deterministic, monotonic optimization algorithm; however, it is susceptible to local maxima!
Hence if we begin the algorithm at different starting points, we risk getting different clustering results.

Image Source: CLP-1

Smart/multiple starts

One solution to problem 1 is is start your algorithm and several different initialization
Another potential solution it to start in a smart location (to be determined by some other clustering method, eg. kmeans)

Problem 2

The total number of “free” parameters in a GMM is: $(G - 1) + G p + G p (p + 1) / 2$

$G p (p + 1) / 2$ of these are from the group covariance matrices $Σ_{1}, . . ., Σ_{G}$
For the EM, we have to invert these $G$ covariance matrices for each iteration of the algorithm.
This can get quite computationally expensive…

Common Covariance

One potential solution to problem 2 is to assumes that the covariance are equal across groups
The constraint of $Σ_{g} = Σ$ for all $g$ this reduces to $p (p + 1) / 2$ free parameters in the common covariance matrix.
This however is a very strong assumption about the underlying data, so perhaps we could impose a softer constraint …

Towards MCLUST

MCLUST exploits an eigenvalue decomposition of the group covariance matrices to give a wide range of structures that uses between 1 and $G p (p + 1) / 2$ free parameters.
The eigenvalue decomposition has the form $Σ_{g} = λ_{g} D_{g} A_{g} D_{g}^{'}$ where $λ_{g}$ is a constant, $D_{g}$ is a matrix consisting of the eigenvectors of $Σ_{g}$ and $A_{g}$ is a diagonal matrix with entries proppotional the the eigenvalues of $Σ_{g}$ .

Constraints & MCLUST

We may constrain the components of the eigenvalue decomposition of $Σ_{g}$ across groups of the mixture model
Furthermore, the matrices $D_{g}$ and $A_{g}$ may be set to the identity matrix $I$
Fraley & Ratery (2000) describe the mclust software, which is available as an R package that efficiently performs this constrained model-based clustering.

MCLUST models

Notes: E = equal; S=spherical; V=variable; AA=axis-aligned
ID	Volume	Shape	Orientation	Covariance Decomposition	Number of free Covariance parameters
EII	E	S		$λ I$	1
VII	V	S		$λ_{g} I$	$G$
EEI	E	E	AA	$λ A$	$p$
VEI	V	E	AA	$λ_{g} A$	$p + G - 1$
EVI	E	V	AA	$λ A_{g}$	$p G - G + 1$
VVI	V	V	AA	$λ_{g} A_{g}$	$p G$
EEE	E	E	E	$λ D A D^{'}$	$p (p + 1) / 2$
EEV	E	E	V	$λ D_{g} A D_{g}^{'}$	$G p (p + 1) / 2 - (G - 1) p$
VEV	V	E	V	$λ_{g} D_{g} A D_{g}^{'}$	$G p (p + 1) / 2 - (G - 1) (p - 1)$
VVV	V	V	V	$λ_{g} D_{g} A_{g} D_{g}^{'}$	$G p (p + 1) / 2$

MCLUST Cluster Shapes

Cluster shapes that correspond to the covariance structures on the previous slide

faithful finite mixture 2D

mod2 <- Mclust(faithful)
summary(mod2, parameters = TRUE)

---------------------------------------------------- 
Gaussian finite mixture model fitted by EM algorithm 
---------------------------------------------------- 

Mclust EEE (ellipsoidal, equal volume, shape and orientation) model with 3
components: 

 log-likelihood   n df       BIC       ICL
      -1126.326 272 11 -2314.316 -2357.824

Clustering table:
  1   2   3 
 40  97 135 

Mixing probabilities:
        1         2         3 
0.1656784 0.3563696 0.4779520 

Means:
               [,1]      [,2]      [,3]
eruptions  3.793066  2.037596  4.463245
waiting   77.521051 54.491158 80.833439

Variances:
[,,1]
           eruptions    waiting
eruptions 0.07825448  0.4801979
waiting   0.48019785 33.7671464
[,,2]
           eruptions    waiting
eruptions 0.07825448  0.4801979
waiting   0.48019785 33.7671464
[,,3]
           eruptions    waiting
eruptions 0.07825448  0.4801979
waiting   0.48019785 33.7671464

The “best” model is 3 components

This assumes the EEE covariance structure

But how was this determined?

BIC plot

plot(mod2$BIC)

By default, MCLUST will fit a g-component mixture model for all 14 models models for $g$ = 1, 2, … 9.

For each of these 126 models, a BIC value is calculated.

The BIC

The Bayesian Information Criterion or BIC (Schwarz, 1978) is one of the most widely known and used tools in statistical model selection. The BIC can be written: $BIC = 2 ℓ (x, \hat{θ}) - p_{f} \log n$ where $ℓ$ is the log-likelihood, $\hat{θ}$ is the MLE of $θ$ , $p_{f}$ is the # of free parameters and $n$ is the # of observations.
It can be shown that the BIC gives consistent estimates of the number of components in a mixtures model under certain regulator conditions (Keribin, 1988 , 2000)

BIC plot (zoomed)

plot(mod2$BIC, ylim=c(-2400, -2300))

The BIC plot reveals that the model which obtains the largest BIC value is the 3-component mixture model with the ‘EEE’ covariance structure.

Pros

Why Mixture Model-based Clustering?

Mathematically and statistically sound using standard theory
Parameter estimation is built-in, along with number of clusters via model selection
Provides a probabilistic answer to the question “does observation $i$ belong to group $G$ ?”
Adaptable, eg. data with outliers/skewness can be handled by non-Gaussian densities

Cons

Cons (a.k.a. open avenues for research)

Computationally intensive (takes time to estimate and invert covariance matrices; not feasible for small $n$ large $p$ without some decomposition)
Model-selection method is still viewed as an open problem by many practitioners — BIC grudgingly accepted
Traditional parameter estimation via EM algorithm prone to local maxima, overconfidence in clustering probabilities, and several other issues