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?

There is very little impact on the posterior distribution of Sigma, since it is simply an error term that does not have profound impact on how the model itself behaves. We can see from the coefficient plot that the two sigmas are almost identical. Another reason is the sigma is greater than 0, so dexp(1) and dunif(0, 1) don't make a big difference.
# Plotting priors
par(mfrow=c(1,2))
plot(density(rexp(n=10000, 1)), main="Sigma ~ dexp(1)")
plot(density(runif(n=10000, 1)), main="Sigma ~ dunif(0, 1)")

# Data cleanup
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)
)

# Model
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 ~ dexp(1)
  ), data=dat_slim, chains=1, log_lik=TRUE
)
## 
## SAMPLING FOR MODEL 'bb40ad6ecc73307845e2f6ff30698be8' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.15 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.041835 seconds (Warm-up)
## Chain 1:                0.028685 seconds (Sampling)
## Chain 1:                0.07052 seconds (Total)
## Chain 1:
precis(m9.1, depth=2)
##             mean          sd        5.5%       94.5%    n_eff     Rhat4
## a[1]   0.8876486 0.016543088  0.86227290  0.91321549 738.2927 0.9980012
## a[2]   1.0504719 0.010025340  1.03419001  1.06701222 768.0456 0.9994057
## b[1]   0.1360938 0.079293007  0.01685955  0.25599414 732.7361 0.9986346
## b[2]  -0.1402565 0.057994899 -0.22910680 -0.04642876 507.2201 0.9983379
## sigma  0.1114906 0.005992905  0.10278865  0.12088684 805.4459 0.9991742
m9.1.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=1, log_lik=TRUE
)
## 
## SAMPLING FOR MODEL 'b9b9ce68d37d9c63639a5ccefcede5fd' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.12 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.046126 seconds (Warm-up)
## Chain 1:                0.029324 seconds (Sampling)
## Chain 1:                0.07545 seconds (Total)
## Chain 1:
precis(m9.1.1, depth=2)
##             mean          sd         5.5%       94.5%     n_eff     Rhat4
## a[1]   0.8859268 0.017060840  0.859090263  0.91330195  699.3062 1.0064763
## a[2]   1.0502630 0.010708656  1.033219834  1.06676599 1054.9981 0.9999964
## b[1]   0.1296734 0.078104443  0.002731057  0.25105738  648.8509 0.9997949
## b[2]  -0.1403445 0.055611428 -0.234643794 -0.05444422  587.4740 0.9984915
## sigma  0.1115069 0.006443606  0.101501720  0.12179869  745.8308 0.9999924
compare(m9.1, m9.1.1)
##             WAIC       SE     dWAIC       dSE    pWAIC    weight
## m9.1   -258.7630 14.62862 0.0000000        NA 5.219067 0.5436971
## m9.1.1 -258.4125 14.71856 0.3504708 0.3859944 5.331814 0.4563029
plot(coeftab(m9.1, m9.1.1), par=c("a[1]", "a[2]", "b[1]", "b[2]", "sigma"))

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?

Changing the prior of b[cid] does have a significant impact on the model. As we can see, even though b[1] is very similar to the original model, b[2] is vastly different. Because in the old prior, the mean and variance of b[2] were both less than 0, in our new prior we eliminated the option to sample from negative numbers, and b[2] is bound to have a positive effect on the outcome, which is not a good prior at all.
# Plotting priors
par(mfrow=c(1,2))
plot(density(rnorm(n=10000, mean=0, sd=0.3)), main="b[cid] ~ dnorm(0, 0.3)")
plot(density(rexp(n=10000, 0.3)), main="b[cid] ~ dexp(0.3)")

m9.1.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=1, log_lik=TRUE
)
## 
## SAMPLING FOR MODEL '83ec4a9a932102ba02ab973023dbd90c' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 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.094682 seconds (Warm-up)
## Chain 1:                0.036604 seconds (Sampling)
## Chain 1:                0.131286 seconds (Total)
## Chain 1:
precis(m9.1.2, depth=2)
##             mean          sd        5.5%      94.5%    n_eff     Rhat4
## a[1]  0.88807790 0.014966078 0.863840375 0.91200177 456.1243 1.0052173
## a[2]  1.04761724 0.009749892 1.032880514 1.06456068 393.2458 1.0016368
## b[1]  0.15132262 0.070996340 0.046642057 0.26808545 268.3307 1.0034257
## b[2]  0.01786205 0.017698756 0.001243129 0.05044324 570.4102 0.9989167
## sigma 0.11440288 0.006500800 0.104083491 0.12522974 306.2590 0.9985752
compare(m9.1, m9.1.2)
##             WAIC       SE   dWAIC      dSE    pWAIC    weight
## m9.1   -258.7630 14.62862 0.00000       NA 5.219067 0.9244312
## m9.1.2 -253.7547 14.36687 5.00827 6.982703 3.257702 0.0755688
plot(coeftab(m9.1, m9.1.2), par=c("a[1]", "a[2]", "b[1]", "b[2]", "sigma"))

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?

We can see that with sampling iteration being 500, with warmup = 5, the model performed significantly worse. However, 10 warmup iterations performs very similarly to 900 warmup iterations.
m9.1.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=1, log_lik=TRUE, iter=505, warmup=5
)
## 
## SAMPLING FOR MODEL 'bb40ad6ecc73307845e2f6ff30698be8' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: WARNING: No variance estimation is
## Chain 1:          performed for num_warmup < 20
## Chain 1: 
## Chain 1: Iteration:   1 / 505 [  0%]  (Warmup)
## Chain 1: Iteration:   6 / 505 [  1%]  (Sampling)
## Chain 1: Iteration:  55 / 505 [ 10%]  (Sampling)
## Chain 1: Iteration: 105 / 505 [ 20%]  (Sampling)
## Chain 1: Iteration: 155 / 505 [ 30%]  (Sampling)
## Chain 1: Iteration: 205 / 505 [ 40%]  (Sampling)
## Chain 1: Iteration: 255 / 505 [ 50%]  (Sampling)
## Chain 1: Iteration: 305 / 505 [ 60%]  (Sampling)
## Chain 1: Iteration: 355 / 505 [ 70%]  (Sampling)
## Chain 1: Iteration: 405 / 505 [ 80%]  (Sampling)
## Chain 1: Iteration: 455 / 505 [ 90%]  (Sampling)
## Chain 1: Iteration: 505 / 505 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.000295 seconds (Warm-up)
## Chain 1:                0.011575 seconds (Sampling)
## Chain 1:                0.01187 seconds (Total)
## Chain 1:
m9.1.4 = 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=1, log_lik=TRUE, iter=510, warmup=10
)
## 
## SAMPLING FOR MODEL 'bb40ad6ecc73307845e2f6ff30698be8' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: WARNING: No variance estimation is
## Chain 1:          performed for num_warmup < 20
## Chain 1: 
## Chain 1: Iteration:   1 / 510 [  0%]  (Warmup)
## Chain 1: Iteration:  11 / 510 [  2%]  (Sampling)
## Chain 1: Iteration:  61 / 510 [ 11%]  (Sampling)
## Chain 1: Iteration: 112 / 510 [ 21%]  (Sampling)
## Chain 1: Iteration: 163 / 510 [ 31%]  (Sampling)
## Chain 1: Iteration: 214 / 510 [ 41%]  (Sampling)
## Chain 1: Iteration: 265 / 510 [ 51%]  (Sampling)
## Chain 1: Iteration: 316 / 510 [ 61%]  (Sampling)
## Chain 1: Iteration: 367 / 510 [ 71%]  (Sampling)
## Chain 1: Iteration: 418 / 510 [ 81%]  (Sampling)
## Chain 1: Iteration: 469 / 510 [ 91%]  (Sampling)
## Chain 1: Iteration: 510 / 510 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.00107 seconds (Warm-up)
## Chain 1:                0.046821 seconds (Sampling)
## Chain 1:                0.047891 seconds (Total)
## Chain 1:
m9.1.5 = 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=1, log_lik=TRUE, iter=1400, warmup=900
)
## 
## SAMPLING FOR MODEL 'bb40ad6ecc73307845e2f6ff30698be8' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 1400 [  0%]  (Warmup)
## Chain 1: Iteration:  140 / 1400 [ 10%]  (Warmup)
## Chain 1: Iteration:  280 / 1400 [ 20%]  (Warmup)
## Chain 1: Iteration:  420 / 1400 [ 30%]  (Warmup)
## Chain 1: Iteration:  560 / 1400 [ 40%]  (Warmup)
## Chain 1: Iteration:  700 / 1400 [ 50%]  (Warmup)
## Chain 1: Iteration:  840 / 1400 [ 60%]  (Warmup)
## Chain 1: Iteration:  901 / 1400 [ 64%]  (Sampling)
## Chain 1: Iteration: 1040 / 1400 [ 74%]  (Sampling)
## Chain 1: Iteration: 1180 / 1400 [ 84%]  (Sampling)
## Chain 1: Iteration: 1320 / 1400 [ 94%]  (Sampling)
## Chain 1: Iteration: 1400 / 1400 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.072771 seconds (Warm-up)
## Chain 1:                0.02963 seconds (Sampling)
## Chain 1:                0.102401 seconds (Total)
## Chain 1:
compare(m9.1.3, m9.1.4, m9.1.5)
##             WAIC       SE      dWAIC       dSE     pWAIC       weight
## m9.1.4 -259.2255 14.62543  0.0000000        NA  5.003797 6.093015e-01
## m9.1.5 -258.3368 14.72707  0.8887538 0.3598463  5.441463 3.906985e-01
## m9.1.3 -184.4591 12.37524 74.7664716 9.2929589 10.021834 3.544045e-17
plot(coeftab(m9.1.3, m9.1.4, m9.1.5), par=c("a[1]", "a[2]", "b[1]", "b[2]", "sigma"))

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

The posterior simply sampled data from the prior, we can see that from the similarity between the prior plot and posterior plot.
# Plotting priors
par(mfrow=c(1,2))
plot(density(rnorm(n=10000, mean=0, sd=1)), main="a ~ dnorm(0,1)")
plot(density(rcauchy(n=10000, 0, 1)), main="b ~ dcauchy(0,1)")

mp <- ulam(
 alist(
   a ~ dnorm(0,1),
   b ~ dcauchy(0,1)
 ), data=list(y=1) , chains=1
)
## 
## SAMPLING FOR MODEL '3bd3f4d287e9cccab124308e5415245c' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.05 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.013808 seconds (Warm-up)
## Chain 1:                0.009076 seconds (Sampling)
## Chain 1:                0.022884 seconds (Total)
## Chain 1:
plot(mp, depth=2)

traceplot(mp)
## [1] 1000
## [1] 1
## [1] 1000

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?

The traceplot for a follows a nice normal distribution, most of the samples are centered around 0. However, the traceplot for b is all over the place, even though it is still centered at around 0. Cauchy distribution is an example of pathological distribution, which has undefined expected value and variance. It aligns with what we see on the posterior 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_lik=TRUE. Explain the model comparison results.

The model with only Age paramter performed the best, Age and Marriage performed slightly worse due to the causal DAG, the model with Marriage only performed the worst. This is expected, and algined with our quap model results.
data("WaffleDivorce")
d = data.frame(
  D=standardize(WaffleDivorce$Divorce),
  M=standardize(WaffleDivorce$Marriage),
  A=standardize(WaffleDivorce$MedianAgeMarriage)
)

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, chains=1, log_lik=TRUE
)
## 
## SAMPLING FOR MODEL 'bc675254cc2159bb33743bbf24468ab3' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 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.011009 seconds (Warm-up)
## Chain 1:                0.010519 seconds (Sampling)
## Chain 1:                0.021528 seconds (Total)
## Chain 1:
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, chains=1, log_lik=TRUE
)
## 
## SAMPLING FOR MODEL '8b690dbf8986a1e80daeae2b8293f78f' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 8e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.08 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.012762 seconds (Warm-up)
## Chain 1:                0.01104 seconds (Sampling)
## Chain 1:                0.023802 seconds (Total)
## Chain 1:
m5.3 = ulam(
  alist(
    D ~ dnorm(mu, sigma),
    mu <- a + bM * M + bA * A,
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ), data=d, chains=1, log_lik=TRUE
)
## 
## SAMPLING FOR MODEL '76181c4a8cfb1100846be932fc80818b' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 1.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 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.024323 seconds (Warm-up)
## Chain 1:                0.019072 seconds (Sampling)
## Chain 1:                0.043395 seconds (Total)
## Chain 1:
compare(m5.1, m5.2, m5.3)
##          WAIC        SE     dWAIC      dSE    pWAIC       weight
## m5.1 125.2331 12.725166  0.000000       NA 3.421932 0.7526481739
## m5.3 127.4634 12.550383  2.230315 0.884348 4.688593 0.2467659640
## m5.2 139.5497  9.889213 14.316537 9.290484 3.064402 0.0005858621
plot(coeftab(m5.1, m5.2, m5.3), par=c("a", "bA", "bM", "sigma"))