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)
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.8865630 0.015675487  0.86151051  0.9116154
## a[2]   1.0505700 0.009936473  1.03468955  1.0664504
## b[1]   0.1325066 0.074203500  0.01391504  0.2510981
## b[2]  -0.1425770 0.054748686 -0.23007595 -0.0550780
## sigma  0.1094926 0.005935098  0.10000720  0.1189781
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.8865489 0.015675977  0.86149570  0.91160218
## a[2]   1.0505596 0.009936812  1.03467863  1.06644052
## b[1]   0.1324083 0.074205935  0.01381289  0.25100372
## b[2]  -0.1429736 0.054750079 -0.23047484 -0.05547244
## sigma  0.1094963 0.005935590  0.10001011  0.11898255
pairs(m9.1)

pairs(m9.1_unif)

They are nearly identical. This is because we have enough data to get a good estimate of sigma, regardless of the prior. However, if sigma were greater than 1, the uniform prior would not give the same estimate, as that prior gives 0 probability to values outside of the [0, 1] interval.

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.8866034 0.015678801  0.8615457  0.91166119
## a[2]   1.0505712 0.009938560  1.0346875  1.06645498
## b[1]   0.1327230 0.074218091  0.0141082  0.25133789
## b[2]  -0.1424292 0.054760075 -0.2299464 -0.05491204
## sigma  0.1095159 0.005938255  0.1000254  0.11900636
pairs(m9.2_exp)

No posterior distribution difference. The new prior constrains the parameter to be positive, so we see all of the posterior density piled up right next to 0.

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?

formula <- 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.8866034 0.015678801  0.8615457  0.91166119
## a[2]   1.0505712 0.009938560  1.0346875  1.06645498
## b[1]   0.1327230 0.074218091  0.0141082  0.25133789
## b[2]  -0.1424292 0.054760075 -0.2299464 -0.05491204
## sigma  0.1095159 0.005938255  0.1000254  0.11900636

We can see that for this particular model, both the Bulk and Tail level off fairly quickly, meaning that not much warmup is needed in this particular case.So you can see that for this model and this data, very little warmup is needed.Basically anything more than 10 is the same.This isn’t always going to be true,however.For the more complicated non-linear models and especially multi level models,more warmup is often helpful.

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)
# 
# precis(mp, depth=2)
# 
# 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")

This model is just sampling prior distributions.We see that a is sampling as we would expect, given the prior, bouncing around between -3 and 3. The b parameter looks a little odd. However, this is normal for the Cauchy distribution. This distribution has an undefined mean and variance, with very heavy tails. Because of this, it will occasionally jump to extreme values.

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

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_stan <- 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
# )
# 
# m5.2_stan <- 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
# )
# 
# m5.3_stan <- 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 
# )
# 
# compare(m5.1_stan, m5.2_stan, m5.3_stan, func = WAIC)
# 
# precis(m5.3_stan)

##Model m5.1 uses only median age at marriage to predict the divorce rate, m5.2 uses only marriage rate, and m5.3 uses both predictors.The first model with only age at marriage is preferred. However, the score for m5.3 is nearly the same. This indicates that the two models are not really being differentiated very well.