The theme of this week is statistical modeling. Given the prerequisites, you have seen statistical models already, as well as the method of maximum likelihood estimation (MLE) to obtain a set of model parameters which best describe the data. If not, here is a refresher on statistical models. A straightforward approach to obtain the MLE parameters, is to write out the joint density function for the observed data, and often we then take the logarithm of the likelihood. If this log likelihood has an amenable form, we can compute derivatives with respect to the parameters, set this equal to zero and solve the equations. For example, this procedure gives the sample mean as the MLE for the \(\mu\) parameter of a Normal distribution.

In computational biology, we often would like to specify models and find the MLE for the model parameters, but the above approach does not work due to the complexity of the problem. We can lump many such complex problems together as those sharing a common property: we can formulate these problem as having a latent variable, a special model parameter, such that if we knew its value, it would be much easier to estimate the other model parameters. Typically the latent variable is represented by \(z_i\), the (other) parameters as \(\theta\), and the data itself as \(x_i\). Typically, the z takes discrete values (e.g. it may be the membership of the observations into multiple groups), x is either discrete- or continuous-valued, and \(\theta\) is typically continuous-valued.

Expectation maximization (EM) is a class of algorithms that help us tackle such problems. The algorithm had been proposed in various domains as early as the 1950s, but a paper in 1977 unified the description of the incomplete data problem (equivalent to latent variable) and the EM algorithm, providing a unified statistical context and guarantees on convergence. Note that the authors of the 1977 work do not use z, \(\theta\), and x notation, as above.

Motivation for EM

We would like to maximize the log likelihood with respect to \(\theta\), and we don’t observe \(z_i\) for each observation of data \(x_i\). We typically get stuck in that we cannot distribute the log to all the terms, which would help simplify the equations. This is described in Ramesh Sridharan’s excellent lecture notes for the mixture of Normals model. We usually have a sum of densities that cannot be separated. If z were known, we could split apart the sums based on which x belong to which latent groupings according to z. Once the sums of densities are split apart, we could easily take derivatives of the log likelihood to maximize over the \(\theta\) parameters.

The EM algorithm proceeds by taking the expectation of the log likelihood with respect to the conditional distribution of Z given X and \(\theta^t\), our best guess of the parameters \(\theta\) at step t in the algorithm. This quantity will fill in for our actual objective, which was the log likelihood marginalized over assignments for Z. Because of Jensen’s Inequality, we have that the expected value of the log of a random variable is less than or equal to the log of the expected value. So we work with a quantity that is a lower bound for our actual objective.

The part about Jensen’s Inequality can be remembered by sketching a diagram of the logarithm function, and the density of a random variable X and log(X), here recreated from Sridharan’s notes. The expected values are drawn with arrows. Intuitively, the logarithm stretches out the left tail of the distribution, and squishes the right tail.

By taking the expectation of the log likelihood with respect to Z given X and \(\theta^t\), we obtain better than simply a lower bound: we obtain a function which is bounded above by the log likelihood and equal to it at \(\theta^t\) (the current estimate). One can furthermore show that by optimizing this replacement log likelihood and obtaining \(\theta^{t+1}\), we have \(\ell(\theta^t) \le \ell(\theta^{t+1})\), that is, monotonic increases in the original log likelihood at each iteration (see Andrew Ng’s course notes for a proof). A diagram of the process of iteratively maximizing functions which lower bound the log likelihood can be found here, which is a supplemantary PDF associated with this review on the EM algorithm by Chuong Do and Serafim Batzoglou.

The actual algorithm begins with a random initialization of parameters \(\theta^0\). We then calculate conditional probabilities for Z given X and \(\theta^0\), and use these in calculating the expectation of the log likelihood. These “guesses” for z help to break apart the log likelihood and facilitate maximization. Finally, we obtain a new estimate of \(\theta\) and repeat until convergence. Note that EM is not guaranteed to find global maxima, but will find local maxima. It is therefore sometimes necessary to try multiple random starts.

Ok, that was a lot of text, but most of all I want to show you how it’s useful. Let’s dive into a real example:

Mixture of Normals example

Suppose we have univariate observed data \(x_i\) which we happen to know are being drawn from two Normal components. We will index the two components by c, and the observations by i. We consider a latent variable Z, which is a binomial random variable, such that \(z_i\) tells us which component a given \(x_i\) comes from, where 0 corresponds to the first component and 1 corresponds to the second component. If \(z_i = 0\), this means that observation i comes from the 1st Normal component. Define \(\theta = (\pi_1, \mu_1, \sigma^2_1, \pi_2, \mu_2, \sigma^2_2)\), where the \(\pi_c\) represent the proportion of observations from component c. Then the binomial distribution for Z has a success probability \(\pi_2 = (1 - \pi_1)\).

Before we start with the EM steps, it’s useful to recall that if we observed the data \(x_i\) and the latent variable \(z_i\) we could work on the log likelihood:

\[ \ell^*(\theta) = \sum_i \sum_c \mathbb{1}(z_i = c) \log \left( \pi_c f(x_i; \mu_c, \sigma^2_c) \right) \quad \textrm{* : if we observed} \, z_i \]

However, we don’t observe the \(z_i\) and so if we attempted maximum likelihood estimation, we would be stuck with a sum over components c preventing us from distributing the log and breaking up terms for differentiation with respect to elements of \(\theta\):

\[ \ell(\theta) = \sum_i \log \sum_c \pi_c f(x_i; \mu_c, \sigma^2_c) \]

Let’s see what we get when we write out the expectation of the log likelihood with respect to the conditional distribution of Z given X and \(\theta^t\). For the moment, we focus on \(\mu_c\).

\[ \textrm{E}_{Z|X,\theta^t} \left[ \ell(\theta) \right] \]

\[ = \textrm{E}_{Z|X,\theta^t} \left[ \log \prod_i \prod_c \left( \pi_c f(x_i; \mu_c, \sigma^2_c) \right)^{\mathbb{1}(z_i = c)} \right] \]

\[ = \textrm{E}_{Z|X,\theta^t} \left[ \sum_i \sum_c \mathbb{1}(z_i = c) (\log \pi_c + \log f(x_i; \mu_c, \sigma^2_c)) \right] \]

\[ \sum_i \sum_c \textrm{E}_{Z|X,\theta^t} \mathbb{1}(z_i = c) \left( \log \pi_c + \log \frac{1}{\sqrt{2 \pi \sigma^2_c}} - \frac{1}{2 \sigma^2_c} (x_i - \mu_c)^2 \right) \]

Now differentiate this double sum with respect to one of the parameters, say \(\mu_1\), and set it equal to zero. Remember we have \(z_i = 0\) if observation i is from the 1st Normal component.

\[ \sum_i \textrm{E}_{Z|X,\theta^t} \mathbb{1}(z_i = 0) (x_i - \mu_1) = 0 \]

Define \(\textrm{E}_{Z|X,\theta^t} \mathbb{1}(z_i = 0) \equiv w_i^1\), which we can calculate using our probability model, \(x\), and \(\theta^t\). Again, \(w_i^1\) is the conditional probability that observation i is from component 1, given the data and the current best estimate of parameters.

We end up with:

\[ \sum_i w_i^1 (x_i - \mu_1) = 0 \]

\[ \hat{\mu}^{t+1}_1 = \frac{\sum_i w_i^1 x_i}{\sum_i w_i^1} \]

We obtain a weighted sum of the observations that belong to this component for our estimate of the mean. If we repeat the same for \(\sigma^2_1\), we will obtain a similar weighted sum of the squares \((x_i - \mu_1)^2\).

Running EM on simulated data

Let’s see how this algorithm performs first on simulated data:

n <- 4000
true.mus <- c(0, 5)
true.sigma2 <- c(1, 4)
x <- c(rnorm(n/2, true.mus[1], sqrt(true.sigma2[1])),
       rnorm(n/2, true.mus[2], sqrt(true.sigma2[2])))
hist(x, breaks=50, col="grey")

We start with arbitrary assignments for our parameters:

pi <- c(.5, .5)
mus <- c(-2, 7)
sigma2 <- c(.5, .5)

Given this parameter set, what do we obtain for the condition probabilities for Z?

w <- cbind(pi[1] * dnorm(x, mus[1], sqrt(sigma2[1])),
           pi[2] * dnorm(x, mus[2], sqrt(sigma2[2])))
w <- w / rowSums(w)
head(w)
##      [,1]         [,2]
## [1,]    1 2.367093e-22
## [2,]    1 4.314743e-14
## [3,]    1 1.835610e-13
## [4,]    1 1.915909e-26
## [5,]    1 3.008355e-12
## [6,]    1 1.046684e-15

We can quickly see what the conditional probabilities look like by plotting one column of w:

plot(w[,1])

We can also see how this would give us new parameter estimates:

N <- colSums(w)
N / n # new estimate of pi
## [1] 0.54404 0.45596
1/N * c(sum(w[,1] * x), sum(w[,2] * x)) # new estimate of mu
## [1] 0.09370729 5.38102194
1/N * c(sum(w[,1] * (x - mus[1])^2), # new estimate of sigma2
        sum(w[,2] * (x - mus[2])^2))
## [1] 5.523739 5.535781

Now we build these two steps into a loop. Notice how quickly the estimates go from our initial, quite bad estimates to much better.

niter <- 10
# lay out a plot where we fill in iterative estimates:
library(rafalib)
nullplot(true.mus[1]-3,true.mus[2]+3,0,niter+1,
         ylab="iterations",xlab="estimates of mu, sigma2")
abline(v=true.mus,lty=2)
points(true.mus, rep(0,2), pch=1:2, col="blue")
arrows(true.mus - sqrt(true.sigma2), rep(0,2),
       true.mus + sqrt(true.sigma2), rep(0,2),
       length=.1, code=3, col="blue")
points(mus, rep(1,2), pch=1:2)
arrows(mus - sqrt(sigma2), rep(1,2),
       mus + sqrt(sigma2), rep(1,2),
       length=.1, code=3)
# the EM loop:
for (i in 1:niter) {
  w <- cbind(pi[1] * dnorm(x, mus[1], sqrt(sigma2[1])),
             pi[2] * dnorm(x, mus[2], sqrt(sigma2[2])))
  w <- w / rowSums(w)
  N <- colSums(w)
  pi <- N / n
  mus <- 1/N * c(sum(w[,1] * x), sum(w[,2] * x))
  sigma2 <- 1/N * c(sum(w[,1] * (x - mus[1])^2),
                    sum(w[,2] * (x - mus[2])^2))
  points(mus, rep(i+1,2), pch=1:2)
  arrows(mus - sqrt(sigma2), rep(i+1,2),
         mus + sqrt(sigma2), rep(i+1,2),
         length=.1, code=3)
}

Promoter CpG count

Now let’s try running our mixture of Normals EM algorithm on some human genome sequence data. We will look at the sequence surrounding the transcription start site (TSS) of the genes in the UCSC known gene database. These regions are called promoters, although the exact definition of the size of the promoter may vary for different analyses. Here we will use the definition from the following article:

A genome-wide analysis of CpG dinucleotides in the human genome distinguishes two distinct classes of promoters.

The authors created windows of size 3 kb surrounding the TSS of genes, and then counted the number of CpG’s, which are C followed by G. Note that the count of CpG does not depend on which strand is used, as the sequence is its own reverse complement. The authors normalized the count of CpG’s by the expected count given the amount of C’s and G’s.

From the paper’s abstract:

A striking feature of the human genome is the dearth of CpG dinucleotides (CpGs) interrupted occasionally by CpG islands (CGIs), regions with relatively high content of the dinucleotide. CGIs are generally associated with promoters; genes, whose promoters are especially rich in CpG sequences, tend to be expressed in most tissues.

Further description in the paper:

The abundance of CpG dinucleotides in human DNA is much lower than expected based on the GC content, which results from the inherent mutability of methylated cytosine. Whereas the product of cytosine deamination, uracil, is readily recognized as aberrant and is repaired, the deamination product of methylated cytosine is thymine, leading to transition mutations in the next round of replication. Consequently, methylated CpGs in the germ line are likely to be lost over time. Ostensibly, CGIs are retained because their CpGs are hypomethylated in the germ line…

We will go over the meaning of this section of the paper in class, but the main idea is that the genes which are methylated in the germ line lose CpG’s, whereas those which are not methylated keep their CpG’s. Very loosely, one can think of genes which are methylated in germ line as having specialized roles in differentiated cell-types (e.g. tissue-specific genes). The paper states that:

Evidence from other studies suggests that CGIs are more frequently associated with “house-keeping” genes than with tissue-specific genes.

Let’s take a look at the distribution of normalized CpG count in the human genome, using UCSC known genes:

library(BSgenome.Hsapiens.UCSC.hg38)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
suppressMessages({
  g <- genes(TxDb.Hsapiens.UCSC.hg38.knownGene)
})
pro <- promoters(g, upstream=1500, downstream=1500)
seq <- getSeq(Hsapiens, pro)
gc <- as.vector( letterFrequency(seq, letters="GC") / 3000 )
cpg <- as.vector( vcountPDict(PDict("CG"), seq) / 3000 )
norm.cpg <- cpg / (gc/2)^2

The following histogram recreates Fig 2A from the paper above.

hist(norm.cpg, breaks=50, col="grey")

We can use the same EM algorithm to fit two densities to the observed bimodal distribution:

pi <- c(.5, .5)
mus <- c(0, 1)
sigma2 <- c(.1, .1)
x <- norm.cpg
n <- length(x)
niter <- 30
for (i in 1:niter) {
  w <- cbind(pi[1] * dnorm(x, mus[1], sqrt(sigma2[1])),
             pi[2] * dnorm(x, mus[2], sqrt(sigma2[2])))
  w <- w / rowSums(w)
  N <- colSums(w)
  pi <- N / n
  mus <- 1/N * c(sum(w[,1] * x), sum(w[,2] * x))
  sigma2 <- 1/N * c(sum(w[,1] * (x - mus[1])^2),
                    sum(w[,2] * (x - mus[2])^2))
}

The paper reported means of 0.23 and 0.61, with standard deviations 0.07 and 0.14 for the two components. We can see we get very similar parameter estimates. The paper reported the lower CpG component with \(\pi = 0.28\), though this depends heavily on the set of genes used for creating promoter regions.

round(mus, 3)
## [1] 0.219 0.605
round(sqrt(sigma2), 3)
## [1] 0.072 0.143
round(pi, 3)
## [1] 0.351 0.649

We can also plot the estimated mixture components onto the histogram to assess our fit qualitatively:

hist(x, breaks=50, col="grey", freq=FALSE)
s <- seq(0,1,length=200)
lines(s, pi[1]*dnorm(s, mus[1], sqrt(sigma2[1])), col="dodgerblue", lwd=4)
lines(s, pi[2]*dnorm(s, mus[2], sqrt(sigma2[2])), col="green4", lwd=4)