Session 3: Introduction to MCMC in R (Computing Practical)

Example 1: sampling from an exponential distribution using MCMC

Plot the exponential distribution:

x <- seq(from = 0, to = 10, length=20)

Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  EOF within quoted string

y <- exp(-x)
plot(x,y)
lines(x,y) #lines adds lines to an existing plot

Now: set up a target function we want to sample from In this case we are going to use the density of the exponential distribution

target <- function(x){
  if (x < 0) {
    return(0)
  } else {
    return(exp(-x))
  }
}

Try computing a couple of values

target(3)

[1] 0.04978707

target(-1)

[1] 0

Next, we will program a Metropolis–Hastings scheme to sample from a distribution proportional to the target

x <- rep(0, 1000)
x[1] <- 3     # this is just a starting value, which I’ve set arbitrarily to 3
for (i in 2:1000) {
  currentx <- x[i - 1]
  proposedx <- currentx + rnorm(1, mean = 0, sd = 1)
  A <- target(proposedx) / target(currentx)
  if (runif(1) < A) {
    x[i] <- proposedx
  } else {
    x[i] <- currentx
  }
}

Note that x is a realisation of a Markov Chain. We can make a few plots of x:

plot(x)

hist(x)

We can wrap this up in a function to make things a bit neater, and make it easy to try changing starting values and proposal distributions

easyMCMC <- function(niter, startval, proposalsd) {
  x <- rep(0, niter)
  x[1] <- startval
  for (i in 2:niter) {
    currentx <- x[i - 1]
    proposedx <- rnorm(1, mean = currentx, sd = proposalsd)
    A <- target(proposedx) / target(currentx)
    if (runif(1) < A) {
      x[i] <- proposedx
    } else {
      x[i] <- currentx
    }
  }
  return(x)
}

Now we’ll run the MCMC scheme 3 times, and look to see how similar the results are

z1 <- easyMCMC(1000,3,1)
z2 <- easyMCMC(1000,3,1)
z3 <- easyMCMC(1000,3,1)
plot(z1, type = "l")
lines(z2, col = 2)

lines(z3, col = 3)

par(mfcol = c(3, 1)) # rather odd command tells R to put 3 graphs on a single page
hist(z1)
hist(z2)

hist(z3)

Exercises: use the function easyMCMC to explore the following:

1a) how do different starting values affect the MCMC scheme?

sv <- seq(1, 50, 1)
z1 <- list()
for (i in 1:length(sv)) {
    z1[[i]] <- easyMCMC(1000, sv[i], 1)
}
par(mfrow = c(1, 2))
plot(sv, as.numeric(lapply(z1, mean)), xlab = "starting value", ylab = "posterior mean")
plot(sv, as.numeric(lapply(z1, var)), xlab = "starting value", ylab = "posterior variance")

par(mfrow = c(1, 3))
hist(z1[[1]])
hist(z1[[25]])

hist(z1[[50]])

par(mfrow = c(1, 3))
plot(z1[[1]], type = "l")
plot(z1[[20]], type = "l")

plot(z1[[50]], type = "l")

1b) what is the effect of having a bigger/smaller proposal standard deviation?

prop_sd <- seq(1,100,1)
z2 <- list()
for (i in 1:length(prop_sd)) {
    z2[[i]] <- easyMCMC(1000, 3, prop_sd[i])
}
par(mfrow = c(1, 2))
plot(prop_sd, as.numeric(lapply(z2,mean)), xlab = "proposal_sd", ylab = "posterior mean")
plot(prop_sd, as.numeric(lapply(z2,var)), xlab = "proposal_sd", ylab = "posterior variance")

par(mfrow = c(1, 3))
hist(z2[[1]])
hist(z2[[25]])

hist(z2[[50]])

par(mfrow = c(1, 3))
plot(z2[[1]], type = "l")
plot(z2[[20]], type = "l")

plot(z2[[50]], type = "l")

A Bimodal Target

Try changing the target function to the following

target <- function(x){
return((x > 0 & x < 1) + (x > 2 & x < 3))
}

What does this target look like? BIMODAL

x <- seq(-1,6,0.01)
plot(x, target(x), type = "l")

There is positive probability only for values between 0 and 1 and between 2 and 3.

z1 <- easyMCMC(1000, 2.5, 1)
z2 <- easyMCMC(1000, 0.7, 1)
z3 <- easyMCMC(1000, 0.5, 1)
par(mfrow = c(1, 3))
hist(z1)
hist(z2)

hist(z3)

And look at the points visited

par(mfrow = c(1, 3))
plot(z1)
plot(z2)

plot(z3)

What happens if the proposal sd is too small here? (try eg 1 and 0.1) values within the two groups work

z1 <- easyMCMC(1000, 2.5, 0.1)
z2 <- easyMCMC(1000, 0.5, 0.1)
par(mfrow = c(1, 3))
hist(z1, xlim = c(0, 3))
hist(z2, xlim = c(0, 3))

plot(z1, type = "l", ylim = c(0, 3))
points(z2, type = "l", col = "dark grey")

Try starting it on a value that is not in (0,1) nor (2,3). It fails. Why?

z1 <- easyMCMC(1000, 1.5, 0.1)
hist(z1)
plot(1:1000, z1)

Example 2: Estimating an allele frequency

A standard assumption when modelling genotypes of bi-allelic loci (eg loci with alleles \(A\) and \(a\)) is that the population is “randomly mating”. From this assumption it follows that the population will be in “Hardy Weinberg Equilibrium” (HWE), which means that if \(p\) is the frequency of the allele \(A\) then the genotypes \(AA\), \(Aa\) and \(aa\) will have frequencies \(p^2\), \(2p(1-p)\), and \((1-p)^2\).

A simple prior for \(p\) is to assume it is uniform on \((0,1)\). Suppose that we sample \(n\) individuals, and observe \(n_{AA}\) with genotype \(AA\), \(n_{Aa}\) with genotype \(Aa\) and \(n_{aa}\) with genotype \(aa\).

The following R code gives a short MCMC routine to sample from the posterior distribution of \(p\). Try to go through the code to see how it works.

Note that this is very similar to the example at the end of Session 2, but it is interesting in that the likelihood is given in terms of \(n_{AA}\), \(n_{Aa}\), and \(n_{aa}\) rather than simply in terms of the number of \(A\) and \(a\) alleles found in the sample.

prior <- function(p){
  if((p < 0) || (p > 1)){ # || here means "or"
    return(0)
  } else{ 
    return(1)
  } 
}
likelihood <- function(p, nAA, nAa, naa) {
    return(p ^ (2 * nAA) * (2 * p * (1-p)) ^ nAa * (1 - p) ^ (2 * naa))
}
psampler <- function(nAA, nAa, naa, niter, pstartval, pproposalsd) {
  p <- rep(0, niter)
  p[1] <- pstartval
  for (i in 2:niter) {
    currentp <- p[i - 1]
    newp <- currentp + rnorm(1, 0, pproposalsd)
    A <- prior(newp) * likelihood(newp, nAA, nAa, naa) / (prior(currentp) * likelihood(currentp, nAA, nAa, naa))
    if (runif(1) < A){
      p[i] <- newp
    } else {
      p[i] <- currentp
    }
  }
  return(p)
}

Now run this sample for \(nAA = 50\), \(nAa = 21\), \(naa=29\).

z <- psampler(50,21,29,10000,0.5,0.01)

Now some R code to compare the sample from the posterior with the theoretical posterior (which in this case is available analytically; from lecture 1, since we observed 121 \(A\)’s, and 79 \(a\)’s, out of 200, the posterior for \(p\) is \(\mathrm{Beta}(121+1,79+1)\). See notes from Session 1.

par(mfcol = c(1, 2))
x <- seq(0,1,length = 1000)
hist(z, prob = T, xlim = c(0.4, 0.8), col = "dark grey", ylim = c(0, 20), main = "No Burn in")
lines(x,dbeta(x,122, 80))  # overlays beta density on histogram

# You might also like to discard the first 5000 z’s as "burnin"
# here’s one way in R to select only the last 5000 zs:
hist(z[5001:10000], prob = T, col = "light grey",  xlim = c(0.4, 0.8), ylim = c(0, 20), main = "With Burn in")
lines(x, dbeta(x, 122, 80))  # overlays beta density on histogram

Have fun exploring!

Some things to try: how do the starting point and proposal sd affect things?

Example 3: Estimating an allele frequency and inbreeding coefficient

(If time allows!)

A slightly more complex alternative than HWE is to assume that there is a tendency for people to mate with others who are slightly more closely-related than “random” (as might happen in a geographically-structured population, for example). This will result in an excess of homozygotes compared with HWE. A simple way to capture this is to introduce an extra parameter, the “inbreeding coefficient” \(f\), and assume that the genotypes \(AA\), \(Aa\) and \(aa\) have frequencies \(fp + (1-f)p^2\), \((1-f) 2p(1-p)\), and \(f(1-p) + (1-f)(1-p)^2\).

In most cases it would be natural to treat \(f\) as a feature of the population, and therefore assume f is constant across loci. For simplicity we will consider just a single locus. Note that both \(f\) and \(p\) are constrained to lie between 0 and 1 (inclusive). A simple prior for each of these two parameters is to assume that they are independent, uniform on [0,1]. Suppose that we sample \(n\) individuals, and observe \(n_{AA}\) with genotype \(AA\), \(n_{Aa}\) with genotype \(Aa\) and \(n_{aa}\) with genotype \(aa\).

Exercise: write a short MCMC routine to sample from the joint distribution of \(f\) and \(p\).

Here is what it looks like.

prior <- function(p) {
  if ((p < 0) || (p > 1)) { # || here means "or"
    return(0)
  } else { 
    return(1)
  } 
}
likelihood <- function(p, f, nAA, nAa, naa) {
    return(((f * p) + ((1 - f) * (p ^ 2)))^nAA * 
             ((1 - f) * 2 * p * (1 - p)) ^ nAa * 
             (f * (1 - p) + (1 - f) * ((1 - p) * (1 - p))) ^ naa
           )
}
fpsampler <- function(nAA, nAa, naa, niter, 
                      fstartval, pstartval,
                      fproposalsd, pproposalsd) {
    f <- rep(0,niter)
    p <- rep(0,niter)
    f[1] <- fstartval
    p[1] <- pstartval
    for (i in 2:niter) {
        currentp <- p[i - 1]
        currentf <- f[i - 1]
        newp <- currentp + rnorm(1, 0, pproposalsd)
        A <- prior(newp) * likelihood(newp,currentf,nAA,nAa,naa) /
          (prior(currentp) * likelihood(currentp,currentf,nAA,nAa,naa))
        
        if (runif(1) < A) {
          p[i] <- newp
        } else {
          p[i] <- currentp
        }
        newf <- currentf + rnorm(1, 0, fproposalsd)
        
        A <- prior(newf) * likelihood(p[i],newf,nAA,nAa,naa) / 
          (prior(currentf) * likelihood(p[i],currentf,nAA,nAa,naa))
        
        
        if (runif(1) < A) {
          f[i] <- newf
        } else {
          f[i] <- currentf
        }
    }
    
    return(list(f = f,p = p)) # return a "list" with two elements named f and p
}

Then run it for a sample with \(n_{AA} = 45\), \(n_{Aa} = 10\), and \(n_{aa} = 45\), which clearly is a case where there are far fewer heterozygotes than one would expect under Hardy-Weinberg equilibrium.

z <- fpsampler(45, 10, 45, 10000, 0.8, 0.5, 0.01, 0.01)

Let’s plot those results:

x <- seq(0, 1, length = 1000)
par(mfrow = c(1,2))
hist(z$f[5000:10000],prob = T, xlim = c(0, 1), col = "dark grey", ylim = c(0, 10), 
     main = "f sample", xlab = "posterior estimate")
hist(z$p[5000:10000], prob = T, xlim = c(0, 1), col = "dark grey", ylim = c(0, 10), 
     main = "p sample", xlab = "posterior estimate")

And look at some trace plots:

par(mfrow = c(1, 2))
plot(z$f, type = "l")
plot(z$p, type = "l")

If we want to summarize those samples to a single point for \(p\) and for \(f\) a good choice is the posterior mean, which is just the mean from the sample from the Markov chain.

# mean for p
p_mean <- mean(z$p)
# mean for f
f_mean <- mean(z$f)

And it is always a good idea to simulate new data from the an estimated model, or look at the expected value of data from an estimated model, to see if it looks anything like the data that went into it.

We can do that by finding the expected number of the three genotypes given the posterior mean estimates of \(f\) and \(p\):

N <- 100
c(nAA = ((p_mean ^ 2) + (f_mean * (1 - p_mean) * p_mean)) * N,
  nAa = ((1 - f_mean) * 2 * p_mean * (1 - p_mean)) * N,
  naa = ((f_mean * (1 - p_mean) * p_mean) + ((1 - p_mean)^2)) * N)

     nAA      nAa      naa 
44.23173 10.97835 44.78992 

Yep…that looks like our observed data…

Session Info

sessionInfo()

R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] compiler_3.4.3  backports_1.1.2 magrittr_1.5    rprojroot_1.3-2
 [5] htmltools_0.3.6 tools_3.4.3     base64enc_0.1-3 yaml_2.1.16    
 [9] Rcpp_0.12.16    stringi_1.1.7   rmarkdown_1.8   knitr_1.18     
[13] jsonlite_1.5    stringr_1.3.0   digest_0.6.15   evaluate_0.10.1