Chapter 9 - Markov Chain Monte Carlo

This chapter 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

9-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. Visualize the posteriors for both models. Does the different prior have any detectible influence on the posterior distribution of sigma? Why or why not?

library(tidybayes)

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)

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 ...
## Exponential prior for sigma
m8.3 <- 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(1)
  ) , 
  data=dat_slim , chains=4 , cores=4)

precis(m8.3 , depth=2)
##             mean          sd        5.5%      94.5%    n_eff     Rhat4
## a[1]   0.8864430 0.016083583  0.86059050  0.9123971 3115.539 0.9992618
## a[2]   1.0504920 0.009847275  1.03489837  1.0660036 3472.358 0.9989867
## b[1]   0.1309124 0.075929156  0.01134495  0.2509957 2546.858 1.0006731
## b[2]  -0.1433005 0.054010464 -0.22717900 -0.0581127 3152.074 0.9987046
## sigma  0.1116048 0.006034337  0.10285635  0.1218343 1997.679 1.0004723
pairs(m8.3)

## Uniform prior for sigma
m9.1 <- 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 ~ dunif(0,1)
  ) , 
  data=dat_slim , chains=4 , cores=4)

precis(m9.1 , depth=2)
##             mean          sd         5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8863792 0.015974827  0.860705898  0.91160493 2679.552 0.9994461
## a[2]   1.0505665 0.009802903  1.034791162  1.06588777 3415.629 0.9993598
## b[1]   0.1316338 0.079179675  0.009303561  0.26154217 1529.286 1.0024445
## b[2]  -0.1420801 0.056643085 -0.233510553 -0.05355969 2638.422 1.0001900
## sigma  0.1115311 0.006167395  0.102210672  0.12182211 3106.715 0.9990486
pairs(m9.1)

# In comparison, parameter estimations are very similar. Except for sigma, which has a higher n_eff in the uniform-prior model, the individual model outputs are also quite close, with all parameter estimates in the model with the exponential-prior model on sigma being significantly higher. As a result, we discover that while different priors impact n_eff, posterior distributions are unaffected.

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

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

# In the posterior, there are almost no distinctions. There seemed to have been so much data that the previous was simply overloaded. The uniform prior, on the other hand, produces a broader and less peaky posterior, indicating more uncertainty.

9-2. Modify the terrain ruggedness model again. This time, change the prior for b[cid] to dexp(0.3). What does this do to the posterior distribution? Can you explain it?

m9.2 <- 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] ~ dexp(0.3),
    sigma ~ dexp(1)
  ) , 
  data=dat_slim , chains=4 , cores=4)

precis(m9.2 , depth=2)
##             mean          sd        5.5%      94.5%    n_eff     Rhat4
## a[1]  0.88770052 0.016764998 0.860899533 0.91371665 1546.522 1.0007586
## a[2]  1.04857040 0.010556231 1.032261075 1.06530357 2102.163 0.9996791
## b[1]  0.15089380 0.074190808 0.039328209 0.27808485 1054.174 1.0037941
## b[2]  0.01857028 0.017409480 0.001143584 0.05156502 1968.452 1.0018607
## sigma 0.11424714 0.006032012 0.104861690 0.12437732 1393.575 1.0021685
pairs(m9.2)

# We can see that there isn't much of a difference; nevertheless, the most significant difference is in the b parameter, which appears to have fewer samples and no longer accepts negative values.

9-3. Re-estimate one of the Stan models from the chapter, but at different numbers 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)

m9.3 <- 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(1)
  ),
  data = dat_slim,
  start = start,
  chains = 4, 
  cores = 4,
  iter = 100
)

# define warmup values to run through
warm_list <- c(5, 10, 100, 500, 1000) 

# first make matrix to hold n_eff results
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(m9.3, 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 # columns show parameters, rows show n_eff
##            a[1]       a[2]       b[1]       b[2]      sigma
## 5      41.04688   74.69378   13.07987   27.54887   13.94382
## 10   4138.04982 4337.11006 1745.13790 2665.85524 2557.45368
## 100  4515.09445 4449.07975 1292.06603 2364.02285 2293.04010
## 500  4709.10436 6038.61322 4733.62364 4521.05374 5430.88874
## 1000 5046.16223 5911.42070 5130.62437 5384.67425 4892.02122
# As can be seen, n_eff does not vary significantly after 10 warmup samples. In this situation, a warmup of 10 would suffice.

9-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.027 seconds (Warm-up)
## Chain 1:                0.067 seconds (Sampling)
## Chain 1:                0.094 seconds (Total)
## Chain 1:
precis(mp , depth=2)
##         mean         sd      5.5%    94.5%     n_eff    Rhat4
## a 0.04711378  0.9439197 -1.325154 1.673254  89.68108 1.009323
## b 0.25856370 10.6380647 -6.185702 6.802573 323.91735 1.001232
pairs(mp)

traceplot(mp)
## [1] 1000
## [1] 1
## [1] 1000
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")

# The normal distribution, as can be seen, has been reconstructed well but the cauchy distribution hasn't.

Compare the samples for the parameters a and b. 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?

From the plot we can see plot a should be a normal distribution as the prior is around 0 and spread between 2 and -2, while plot b is Cauchy distribution because it contains some extreme value that’s up to over 30 and -40.

9-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.

# load data and copy
data(WaffleDivorce)
d <- WaffleDivorce

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

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

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_trim,
  chains = 4, 
  cores = 4,
  log_lik = TRUE # this is needed to get the terms for calculating PSIS or WAIC instead log_log
)

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_trim,
  chains = 4, 
  cores = 4,
  log_lik = TRUE # this is needed to get the terms for calculating PSIS or WAIC instead log_log
)

m5.3 <- ulam(
  alist(
    D ~ dnorm(mu, sigma),
    mu <- a + bA * A + bM * M,
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  data = d_trim,
  chains = 4, 
  cores = 4,
  log_lik = TRUE # this is needed to get the terms for calculating PSIS or WAIC instead log_log
)

set.seed(77)
compare( m5.1 , m5.2 , m5.3 , func=WAIC )
##          WAIC        SE    dWAIC      dSE    pWAIC      weight
## m5.1 125.6568 12.691826  0.00000       NA 3.587915 0.742878898
## m5.3 127.7844 12.858339  2.12763 0.742790 4.833452 0.256394688
## m5.2 139.5171  9.876543 13.86034 9.293044 3.058939 0.000726414
# Model m5.1 with only median age at marriage as a predictor performs best, but is not really distinguishable from model m5.3. However, the model with marriage rate only, m5.2 clearly performs worse than both.