Chapter 9 - Markov Chain Monte Carlo

This week has been an informal introduction to Markov chain Monte Carlo (MCMC) estimation. The goal has been to introduce the purpose and approach MCMC algorithms. The major algorithms introduced were the Metropolis, Gibbs sampling, and Hamiltonian Monte Carlo algorithms. Each has its advantages and disadvantages. The ulam function in the rethinking package was introduced. It uses the Stan (mc-stan.org) Hamiltonian Monte Carlo engine to fit models as they are defined in this book. General advice about diagnosing poor MCMC fits was introduced by the use of a couple of pathological examples.

Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them.

Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points.

Questions

8-1. Re-estimate the terrain ruggedness model from the chapter, but now using a uniform prior for the standard deviation, sigma. The uniform prior should be dunif(0,1). Visualize the priors. Use ulam to estimate the posterior distribution of sigma. Visualize the posterior distribution of sigma for both models. Do not use a ‘pairs’ plot. Does the different prior have any detectable influence on the posterior distribution of sigma? Why or why not?

library(usethis)
library(devtools)
library(bayesplot)
library(ggplot2)
library(tibble)
library(rethinking)
library(dplyr)
library(dagitty)
library(tidyr)
library(tidybayes)
library(ggdist)

data(rugged)
d <- rugged
d$log_gdp <- log(d$rgdppc_2000)
dd <- d[ complete.cases(d$rgdppc_2000) , ]
dd$log_gdp_std <- dd$log_gdp / mean(dd$log_gdp)
dd$rugged_std <- dd$rugged / max(dd$rugged)
dd$cid <- ifelse( dd$cont_africa==1 , 1 , 2 )

m9.1 <- quap(
alist(
log_gdp_std ~ dnorm( mu , sigma ) ,
mu <- a[cid] + b[cid]*( rugged_std - 0.215 ) ,
a[cid] ~ dnorm( 1 , 0.1 ) ,
b[cid] ~ dnorm( 0 , 0.3 ) ,
sigma ~ dexp( 1 )
) ,
data=dd )
precis( m9.1 , depth=2 )
##             mean          sd        5.5%       94.5%
## a[1]   0.8865635 0.015675440  0.86151111  0.91161587
## a[2]   1.0505696 0.009936443  1.03468928  1.06644999
## b[1]   0.1325008 0.074203288  0.01390961  0.25109198
## b[2]  -0.1425756 0.054748523 -0.23007436 -0.05507693
## sigma  0.1094923 0.005935052  0.10000693  0.11897765
m9.1_unif <- quap(
  alist(
    log_gdp_std ~ dnorm( mu , sigma ) ,
    mu <- a[cid] + b[cid]* (rugged_std-0.215) ,
    a[cid] ~ dnorm(1,0.1),
    b[cid] ~ dnorm(0,0.3),
    sigma ~ dunif(0,1)
  ) , 
  data=dd)
precis(m9.1_unif , depth=2)
##             mean          sd        5.5%       94.5%
## a[1]   0.8865652 0.015680674  0.86150443  0.91162592
## a[2]   1.0505688 0.009939814  1.03468306  1.06645455
## b[1]   0.1325024 0.074227139  0.01387305  0.25113166
## b[2]  -0.1425730 0.054766659 -0.23010067 -0.05504527
## sigma  0.1095298 0.005940139  0.10003632  0.11902330
pairs(m9.1)

pairs(m9.1_unif)

Plot_d <- data.frame(
  Posteriors = c(
    extract.samples(m9.1, n = 1e4)$sigma,
    extract.samples(m9.1_unif, n = 1e4)$sigma
  ),
  Name = rep(c("Exp", "Uni"), each = 1e4),
  Model = rep(c("m9.1", "m9.1_unif"), each = 1e4)
)

ggplot(Plot_d, aes(y = Model, x = Posteriors)) +
  stat_halfeye() +
  labs(x = "Parameter Estimate", y = "Model") +
  theme_bw()

dat_slim <- list(
log_gdp_std = dd$log_gdp_std,
rugged_std = dd$rugged_std,
cid = as.integer( dd$cid )
)
str(dat_slim)
## List of 3
##  $ log_gdp_std: num [1:170] 0.88 0.965 1.166 1.104 0.915 ...
##  $ rugged_std : num [1:170] 0.138 0.553 0.124 0.125 0.433 ...
##  $ cid        : int [1:170] 1 2 2 2 2 2 2 2 2 1 ...
# Parameter estimates are comparable. Except for sigma, which has a greater n eff in the uniform-prior model, individual model outputs are similar, with all parameter values in the exponential-prior model much higher on sigma. Thus, while various priors alter n eff, posterior distributions are unchanged.

8-2. Modify the terrain ruggedness model again. This time, change the prior for b[cid] to dexp(0.3). Plot the joint posterior. Do not use a ‘pairs’ plot. What does this do to the posterior distribution? Can you explain it?

m9.2_exp <- quap(
  alist(
    log_gdp_std ~ dnorm( mu , sigma ) ,
    mu <- a[cid] + b[cid]* (rugged_std-0.215) ,
    a[cid] ~ dnorm(1,0.1),
    b[cid] ~ dnorm(0,0.3),
    sigma ~ dexp(0.3)
  ) , 
  data=dd)
precis(m9.2_exp , depth=2)
##             mean          sd        5.5%      94.5%
## a[1]   0.8865663 0.015678988  0.86150827  0.9116244
## a[2]   1.0505727 0.009938723  1.03468869  1.0664567
## b[1]   0.1325179 0.074219428  0.01390088  0.2511348
## b[2]  -0.1425710 0.054760809 -0.23008940 -0.0550527
## sigma  0.1095177 0.005938497  0.10002684  0.1190086
#Visualize joint posterior
library(data.table)

b <- extract.samples(m9.2_exp)$b[,1]
sigma <- extract.samples(m9.2_exp)$sigma

b_df <- data.table(x=b)
sigma_df <- data.table(x=sigma)

b_posterior_plot <- ggplot(b_df, aes(x=x)) +
  geom_density(fill="gray", alpha=0.5) +
  ggtitle("Posterior for b") +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x="Posterior", y="")
  
b_posterior_plot

# Not much change is observed. b contains fewer values and takes less negative value. This is the most important difference observed.

8-3. Re-estimate one of the Stan models from the chapter, but at different numbers (at least 5) of warmup iterations. Be sure to use the same number of sampling iterations in each case. Compare the n_eff values. How much warmup is enough?

# use fixed start values for comparability of runs
start <- list(a = c(1, 1), b = c(0, 0), sigma = 1)

m3 <- ulam(
  alist(
    log_gdp_std ~ dnorm(mu, sigma),
    mu <- a[cid] + b[cid] * (rugged_std - 0.215),
    a[cid] ~ dnorm(1, 0.1),
    b[cid] ~ dnorm(0, 0.3),
    sigma ~ dexp(0.3)
  ),
  data = dat_slim,
  start = start,
  chains = 4, 
  cores = 4,
  iter = 100
)

# Warm-up values to run
warm_list <- c(5, 10, 100, 500, 1000) 

# Make matrix to hold n_eff outputs
n_eff <- matrix(NA, nrow = length(warm_list), ncol = 5)

for (i in 1:length(warm_list)) { # loop over warm_list and collect n_eff
  w <- warm_list[i]
  m_temp <- ulam(m3, chains = 4, cores = 4, iter = 1000 + w, warmup = w, start = start)
  n_eff[i, ] <- precis(m_temp, 2)$n_eff
}

colnames(n_eff) <- rownames(precis(m_temp, 2))
rownames(n_eff) <- warm_list
n_eff
##           a[1]       a[2]       b[1]        b[2]      sigma
## 5     198.2246   34.80307   34.57076    8.120709   21.85870
## 10   1747.3564 1414.01891   78.53167  556.632325   80.17199
## 100  4397.7121 4392.84010 1863.69611 2467.675694 2830.43656
## 500  5351.3744 5820.47650 5216.88238 5091.916825 4624.05272
## 1000 5282.3486 5687.19049 5210.38312 5480.474467 4708.38331
# As observed, after 10 warm-up samples, n eff does not change considerably. In this case, a ten-minute warm-up would adequate.

8-4. Run the model below and then inspect the posterior distribution and explain what it is accomplishing.

mp <- ulam(
 alist(
   a ~ dnorm(0,1),
   b ~ dcauchy(0,1)
 ), data=list(y=1) , chains=1 )
## 
## SAMPLING FOR MODEL 'bcf56ee89f6cf2a4224a4139ff01c7d4' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.037 seconds (Warm-up)
## Chain 1:                0.018 seconds (Sampling)
## Chain 1:                0.055 seconds (Total)
## Chain 1:
precis(mp , depth=2)
##           mean        sd      5.5%    94.5%    n_eff     Rhat4
## a -0.005958679 0.9491638 -1.609185 1.479621 141.3132 1.0107814
## b  0.247553816 5.9078048 -5.849092 6.193462 129.2424 0.9995635
pairs(mp)

traceplot(mp)

post <- extract.samples(mp)
par(mfrow = c(1, 2))

dens(post$a)
curve(dnorm(x, 0, 1), from = -4, to = 4, add = T, lty = 2)
legend("topright", lty = c(1, 2), legend = c("Sample", "Exact density"), bty = "n")
mtext("Normal")

dens(post$b, col = "red", xlim = c(-10, 10))
curve(dcauchy(x, 0, 1),
  from = -10, to = 10, add = T, lty = 2,
  col = "red"
)
mtext("Cauchy")

# As seen, the normal distribution has been properly recreated, but the cauchy distribution has not.

Compare the samples for the parameters a and b. Plot the trace plots. Can you explain the different trace plots? If you are unfamiliar with the Cauchy distribution, you should look it up. The key feature to attend to is that it has no expected value. Can you connect this fact to the trace plot?

8-5. Recall the divorce rate example from Chapter 5. Repeat that analysis, using ulam this time, fitting models m5.1, m5.2, and m5.3. Use compare to compare the models on the basis of WAIC or PSIS. To use WAIC or PSIS with ulam, you need add the argument log_log=TRUE. Explain the model comparison results.

data(WaffleDivorce)
d <- WaffleDivorce


# Standardize variables
d$A <- scale(d$MedianAgeMarriage)
d$D <- scale(d$Divorce)
d$M <- scale(d$Marriage)

d_slim <- list(M=d$M, D=d$D, A=d$A)

# Model m5.1
m5.1 <- ulam(
    alist(
        D ~ dnorm(mu , sigma) ,
        mu <- a + bA * A ,
        a ~ dnorm(0 , 0.2) ,
        bA ~ dnorm(0 , 0.5) ,
        sigma ~ dexp(1)
    ), data = d_slim, chains=4, cores=4, log_lik=TRUE)

# Model m5.2
m5.2 <- ulam(
    alist(
        D ~ dnorm(mu, sigma) ,
        mu <- a + bM * M ,
        a ~ dnorm(0 , 0.2) ,
        bM ~ dnorm(0 , 0.5) ,
        sigma ~ dexp(1)
    ), data = d_slim, chains=4, cores=4, log_lik=TRUE)

# Model m5.3
m5.3 <- ulam(
    alist(
        D ~ dnorm(mu, sigma) ,
        mu <- a + bM*M + bA*A ,
        a ~ dnorm(0 , 0.2) ,
        bM ~ dnorm(0 , 0.5) ,
        bA ~ dnorm(0 , 0.5) ,
        sigma ~ dexp(1)
    ), data = d_slim, chains=4, cores=4, log_lik=TRUE)

# WAIC
set.seed(42)
compare(m5.1 , m5.2 , m5.3 , func=WAIC)
##          WAIC        SE     dWAIC       dSE    pWAIC       weight
## m5.1 125.8548 12.793450  0.000000        NA 3.713248 0.6606721744
## m5.3 127.1922 12.372827  1.337384 0.8763441 4.472777 0.3385141673
## m5.2 139.2537  9.816579 13.398945 9.3580020 2.928360 0.0008136583
# Model m5.1 performs best when just median age at marriage is used as a predictor, however it is difficult to distinguish from model m5.3. However, the model with simply the marriage rate, m5.2, obviously outperforms both.