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
mdl1 <- 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(mdl1 , depth=2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8863433 0.016505211  0.85964279  0.91304510 2827.466 0.9996252
## a[2]   1.0504657 0.010164289  1.03422780  1.06685718 2921.220 0.9997632
## b[1]   0.1282445 0.074118471  0.00929291  0.24607185 2483.153 0.9990073
## b[2]  -0.1431307 0.056371658 -0.23407131 -0.05277653 2677.243 0.9985521
## sigma  0.1115293 0.006155348  0.10216348  0.12179748 2573.969 1.0003972
## Uniform prior for sigma
mdl2 <- 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(mdl2 , depth=2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8866276 0.015544313  0.86168601  0.91190974 3106.828 0.9989677
## a[2]   1.0504408 0.009855062  1.03443547  1.06648080 3077.463 0.9992718
## b[1]   0.1317907 0.074200543  0.01457731  0.24787010 2593.137 0.9987949
## b[2]  -0.1416034 0.055180013 -0.22923032 -0.05289532 2315.896 0.9998813
## sigma  0.1118192 0.006341196  0.10175412  0.12212446 2346.930 0.9998674
coeftab(mdl1,mdl2)
##       mdl1    mdl2   
## a[1]     0.89    0.89
## a[2]     1.05    1.05
## b[1]     0.13    0.13
## b[2]    -0.14   -0.14
## sigma    0.11    0.11
## nobs      170     170
#when we compare both models the parameter estimations seems exactly same.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(mdl1, n = 1e4)$sigma,
    extract.samples(mdl2, n = 1e4)$sigma
  ),
  Name = rep(c("Exp", "Uni"), each = 1e4),
  Model = rep(c("mdl1", "mdl2"), each = 1e4)
)

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

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?

mdl3 <- 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(mdl3 , depth=2)
##             mean          sd       5.5%      94.5%     n_eff     Rhat4
## a[1]  0.88749290 0.016330170 0.86205050 0.91452857 2056.0421 0.9993134
## a[2]  1.04813933 0.010162424 1.03252913 1.06437859 2383.2816 0.9992909
## b[1]  0.14444123 0.073483788 0.02933570 0.26760119  860.5506 1.0014148
## b[2]  0.01887009 0.017350235 0.00133116 0.05218188 2028.4094 0.9997363
## sigma 0.11398220 0.006389134 0.10460598 0.12455943 1854.5017 1.0005579
Plot_df <- data.frame(
  Posteriors = extract.samples(mdl3, n = 1e4)$sigma,
  Name = rep( "Uni", each = 1e4),
  Model = rep( "mdl3", each = 1e4)
)

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

#Looking at the above plot Modifying the uniform prior model and changing the variable did not show any impact on posteriors.

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 = 100 + 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     29.55489  83.28189   3.810239   3.604834   8.88348
## 10   405.11386 416.46215 161.987316 235.608644 296.52715
## 100  423.07611 464.17887 135.315365 279.819389 294.01076
## 500  511.93142 540.84314 560.142824 620.897916 613.68157
## 1000 476.92659 589.88718 643.604320 710.415039 450.12628

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.016 seconds (Warm-up)
## Chain 1:                0.016 seconds (Sampling)
## Chain 1:                0.032 seconds (Total)
## Chain 1:
precis(mp)
##         mean        sd      5.5%    94.5%    n_eff    Rhat4
## a 0.03261893 0.9812569 -1.499700 1.594843 221.3314 1.012867
## b 0.02196166 3.9521974 -3.947496 3.478616 159.2457 1.004512
rstan::traceplot(mp, n_cols = 1, col = "Purple")
## [1] 1000
## [1] 1
## [1] 1000

# The prior is about 0 and the spread is between 2 and -2, thus plot 'a' should be a normal distribution. Plot 'b' depicts a Cauchy distribution with several extreme values ranging from 30 to -40.
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 = "Dark blue", xlim = c(-10, 10))
curve(dcauchy(x, 0, 1),
  from = -10, to = 10, add = T, lty = 2,
  col = "red"
)
mtext("Cauchy")

# Looking at the plot we can see that the normal distribution 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?

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 
)

set.seed(77)
compare( m5.1 , m5.2 , m5.3 , func=WAIC )
##          WAIC       SE     dWAIC       dSE    pWAIC      weight
## m5.1 125.3259 12.45409  0.000000        NA 3.407594 0.717424766
## m5.3 127.1948 12.49316  1.868885 0.7586439 4.536678 0.281807908
## m5.2 139.0069  9.66392 13.681023 9.0426102 2.791354 0.000767326
compare(m5.1, m5.2, m5.3, func = PSIS)
##          PSIS        SE    dPSIS       dSE    pPSIS       weight
## m5.1 125.7880 12.923993  0.00000        NA 3.638628 0.6804313433
## m5.3 127.3050 12.651498  1.51709 0.9000941 4.591814 0.3186782669
## m5.2 139.0656  9.781889 13.27765 9.4031260 2.820699 0.0008903898
#WAIC and PSIS have almost similar coefficients. In model 3 marriage predictor is not statistically significant which means models 1 and 2 are also making same predictions