# Practice: Chp. 14#

source("iplot.R")
suppressPackageStartupMessages(library(rethinking))


14E1. Add to the following model varying slopes on the predictor x.

\begin{split} \begin{align} y_i & \sim Normal(\mu_i, \sigma) \\ \mu_i & = \alpha_{GROUP[i]} + \beta x_i \\ \alpha_{GROUP} & \sim Normal(\alpha, \sigma_{\alpha}) \\ \alpha & \sim Normal(0, 10) \\ \beta & \sim Normal(0, 1) \\ \sigma & \sim Exponential(1) \\ \sigma_{\alpha} & \sim Exponential(1) \end{align} \end{split}

Answer. See section 14.1.3. Notice the distinction of $$i$$ and $$j$$ indexes.

\begin{split} \begin{align} y_i & \sim Normal(\mu_i, \sigma) \\ \mu_i & = \alpha_{GROUP[j]} + \beta_{GROUP[j]} x_i \\ \begin{bmatrix} \alpha_{GROUP} \\ \beta_{GROUP} \end{bmatrix} & \sim MVNormal( \begin{bmatrix} \alpha \\ \beta \end{bmatrix} , \mathbf{S}) \\ \alpha & \sim Normal(0, 10) \\ \beta & \sim Normal(0, 1) \\ \mathbf{S} & = \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{pmatrix} \mathbf{R} \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{pmatrix} \\ \sigma_{\alpha} & \sim Exponential(1) \\ \sigma_{\beta} & \sim Exponential(1) \end{align} \end{split}

14E2. Think up a context in which varying intercepts will be positively correlated with varying slopes. Provide a mechanistic explanation for the correlation.

Answer. If you believe that education leads to greater wealth, then in the prediction of wealth based on parentās wealth. If you come from family with a lot of money youāll start off with a lot of money (the intercept), and because wealth builds on wealth the addition of education will help more than for someone who doesnāt have as many resources to work with to start (the slope). You could replace education with ambition or something else along that theme, as well.

We could also adapt the cafe example to use an M indicator (for morning) rather than an A indicator for afternoon. The intercepts would have to become the afternoon wait. Long afternoon waits are correlated with even longer morning waits.

14E3. When is it possible for a varying slopes model to have fewer effective parameters (as estimated by WAIC or PSIS) than the corresponding model with fixed (unpooled) slopes? Explain.

Answer. When there is some relationship between the intercepts and slopes that helps the model regularize itself. That is, if there is some correlation between intercepts and slopes, the model should be able to detect this and effectively learn to predict intercepts from slopes, and vice versa. This is similar to how in an intercepts-only multilevel model knowing some intercepts should help the model predict other intercepts; the known intercepts are a regularizing prior for new intercepts.

14M1. Repeat the cafĆ© robot simulation from the beginning of the chapter. This time, set rho to zero, so that there is no correlation between intercepts and slopes. How does the posterior distribution of the correlation reflect this change in the underlying simulation?

a <- 3.5 # average morning wait time
b <- (-1) # average difference afternoon wait time
sigma_a <- 1 # std dev in intercepts
sigma_b <- 0.5 # std dev in slopes
rho <- (-0.0) # correlation between intercepts and slopes

Mu <- c(a, b)

sigmas <- c(sigma_a, sigma_b) # standard deviations
Rho <- matrix(c(1, rho, rho, 1), nrow = 2) # correlation matrix
Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas)

N_cafes <- 20

library(MASS)
set.seed(5) # used to replicate example
vary_effects <- mvrnorm(N_cafes, Mu, Sigma)

a_cafe <- vary_effects[, 1]
b_cafe <- vary_effects[, 2]

iplot(function() {
plot(a_cafe, b_cafe,
col = rangi2,
xlab = "intercepts (a_cafe)", ylab = "slopes (b_cafe)",
main = "Figure 14.2 (updated simulation)"
)

library(ellipse)
for (l in c(0.1, 0.3, 0.5, 0.8, 0.99)) {
lines(ellipse(Sigma, centre = Mu, level = l), col = col.alpha("black", 0.2))
}
})

set.seed(22)
N_visits <- 10
afternoon <- rep(0:1, N_visits * N_cafes / 2)
cafe_id <- rep(1:N_cafes, each = N_visits)
mu <- a_cafe[cafe_id] + b_cafe[cafe_id] * afternoon
sigma <- 0.5 # std dev within cafes
wait <- rnorm(N_visits * N_cafes, mu, sigma)
d <- data.frame(cafe = cafe_id, afternoon = afternoon, wait = wait)

Attaching package: āellipseā

The following object is masked from āpackage:rethinkingā:

pairs

The following object is masked from āpackage:graphicsā:

pairs


Refit the model:

set.seed(867530)
m_cafe_rho_zero <- ulam(
alist(
wait ~ normal(mu, sigma),
mu <- a_cafe[cafe] + b_cafe[cafe] * afternoon,
c(a_cafe, b_cafe)[cafe] ~ multi_normal(c(a, b), Rho, sigma_cafe),
a ~ normal(5, 2),
b ~ normal(-1, 0.5),
sigma_cafe ~ exponential(1),
sigma ~ exponential(1),
Rho ~ lkj_corr(2)
),
data = d, chains = 4, cores = 4
)
display(precis(m_cafe_rho_zero, depth = 3), mimetypes="text/plain")

post <- extract.samples(m_cafe_rho_zero) # posterior
R <- rlkjcorr(1e4, K = 2, eta = 2) # prior

iplot(function() {
dens(postRho[, 1, 2], xlim = c(-1, 1)) dens(R[, 1, 2], add = TRUE, lty = 2) })   mean sd 5.5% 94.5% n_eff Rhat4 b_cafe[1] -0.97493827 2.706508e-01 -1.4053396 -0.5321296 2558.033 0.9985400 b_cafe[2] -1.57678103 2.754448e-01 -2.0302763 -1.1429525 2258.533 0.9987245 b_cafe[3] -1.75208873 2.937009e-01 -2.2155046 -1.2700634 2475.211 1.0004769 b_cafe[4] -1.38179526 2.640233e-01 -1.7984554 -0.9628625 2678.364 0.9997784 b_cafe[5] -0.85552684 2.673436e-01 -1.2846377 -0.4451974 2146.033 0.9990816 b_cafe[6] -1.03052158 2.727011e-01 -1.4647225 -0.5959246 2341.112 1.0009258 b_cafe[7] -1.07318113 2.612501e-01 -1.5018916 -0.6632712 2523.437 0.9998771 b_cafe[8] -1.66622936 2.724191e-01 -2.1166805 -1.2528713 2411.146 1.0007548 b_cafe[9] -1.09643654 2.628003e-01 -1.5073631 -0.6833408 2509.108 0.9986496 b_cafe[10] -0.85150181 2.644251e-01 -1.2527569 -0.4179187 2583.283 0.9994003 b_cafe[11] -0.94436493 2.765405e-01 -1.3854714 -0.4974285 2568.146 0.9985809 b_cafe[12] -1.08532110 2.689619e-01 -1.5141136 -0.6673375 2766.285 1.0006618 b_cafe[13] -1.82885505 2.764952e-01 -2.2718830 -1.3952793 2499.467 0.9990739 b_cafe[14] -1.09088505 2.773739e-01 -1.5262770 -0.6556258 2419.601 0.9986678 b_cafe[15] -2.08013935 2.922605e-01 -2.5285539 -1.6230970 2367.273 1.0003216 b_cafe[16] -1.13964911 2.668720e-01 -1.5497444 -0.7027246 2683.230 0.9985100 b_cafe[17] -0.84962462 2.791103e-01 -1.2910448 -0.4098209 2627.050 1.0004932 b_cafe[18] 0.13329806 2.973738e-01 -0.3546918 0.6137691 2271.025 0.9996235 b_cafe[19] -0.05123218 2.914838e-01 -0.5206416 0.4161668 2227.332 0.9997787 b_cafe[20] -0.94603730 2.643219e-01 -1.3730128 -0.5212390 2253.346 0.9992708 a_cafe[1] 4.32822010 2.020550e-01 3.9961308 4.6447876 2655.898 0.9990083 a_cafe[2] 2.23007232 2.012716e-01 1.9124224 2.5524646 2430.606 0.9988290 a_cafe[3] 4.55863435 2.128146e-01 4.2177601 4.8939822 2441.837 1.0012199 a_cafe[4] 3.31432538 1.956863e-01 2.9999354 3.6215312 2700.695 0.9993147 a_cafe[5] 1.92490907 1.986681e-01 1.6207611 2.2437742 2374.147 0.9986663 a_cafe[6] 4.24067936 1.990592e-01 3.9296951 4.5483649 2351.843 0.9995529 a_cafe[7] 3.77722779 1.908971e-01 3.4746905 4.0842010 2392.385 0.9998529 a_cafe[8] 4.12167102 2.033142e-01 3.8058353 4.4586906 2764.740 1.0017864 a_cafe[9] 3.91765196 1.923953e-01 3.6088796 4.2293424 2563.342 0.9989676 a_cafe[10] 3.46479236 2.000028e-01 3.1430532 3.7814975 2343.924 0.9993257 a_cafe[11] 1.95604246 2.003907e-01 1.6292569 2.2770210 2497.175 0.9988287 a_cafe[12] 3.98535177 1.923491e-01 3.6838910 4.2991258 2487.950 0.9997389 a_cafe[13] 4.14900451 2.017190e-01 3.8299155 4.4610220 2534.936 0.9989087 a_cafe[14] 3.31122820 2.030783e-01 2.9852590 3.6312353 2594.142 0.9994834 a_cafe[15] 4.62640653 2.048657e-01 4.3016374 4.9457133 2585.030 1.0030749 a_cafe[16] 3.48916086 1.981373e-01 3.1825172 3.8082679 2782.811 1.0004568 a_cafe[17] 4.13659311 2.049720e-01 3.8173995 4.4642264 2521.503 0.9996007 a_cafe[18] 5.57907198 2.139529e-01 5.2380410 5.9266490 2419.397 0.9987837 a_cafe[19] 3.07139300 2.068774e-01 2.7288774 3.4107189 2409.187 0.9994801 a_cafe[20] 3.72791143 1.933597e-01 3.4257859 4.0296822 2214.056 0.9992702 a 3.70381547 2.298564e-01 3.3455440 4.0634464 2344.457 1.0003163 b -1.09575262 1.568141e-01 -1.3393037 -0.8450130 2020.744 0.9996859 sigma_cafe[1] 0.98523622 1.841260e-01 0.7384284 1.3121756 2191.504 1.0001043 sigma_cafe[2] 0.63873123 1.387958e-01 0.4412895 0.8792464 1888.167 0.9996571 sigma 0.47388432 2.663314e-02 0.4333083 0.5167588 2110.001 0.9995346 Rho[1,1] 1.00000000 0.000000e+00 1.0000000 1.0000000 NaN NaN Rho[1,2] -0.04711170 2.294524e-01 -0.4003752 0.3325041 1973.052 0.9995424 Rho[2,1] -0.04711170 2.294524e-01 -0.4003752 0.3325041 1973.052 0.9995424 Rho[2,2] 1.00000000 9.264531e-17 1.0000000 1.0000000 1999.376 0.9979980  The model has inferred rho is zero, as expected. 14M2. Fit this multilevel model to the simulated cafĆ© data: \begin{split} \begin{align} W_i & \sim Normal(\mu_i, \sigma) \\ \mu_i & = \alpha_{CAFE[j]} + \beta_{CAFE[j]} A_i \\ \alpha_{CAFE} & \sim Normal(\alpha, \sigma_{\alpha}) \\ \beta_{CAFE} & \sim Normal(\beta, \sigma_{\beta}) \\ \alpha & \sim Normal(0, 10) \\ \sigma, \sigma_{\alpha}, \sigma_{\beta} & \sim Exponential(1) \end{align} \end{split} Use WAIC to compare this model to the model from the chapter, the one that uses a multi-variate Gaussian prior. Explain the result. a <- 3.5 # average morning wait time b <- (-1) # average difference afternoon wait time sigma_a <- 1 # std dev in intercepts sigma_b <- 0.5 # std dev in slopes rho <- (-0.7) # correlation between intercepts and slopes Mu <- c(a, b) sigmas <- c(sigma_a, sigma_b) # standard deviations Rho <- matrix(c(1, rho, rho, 1), nrow = 2) # correlation matrix Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas) N_cafes <- 20 library(MASS) set.seed(5) # used to replicate example vary_effects <- mvrnorm(N_cafes, Mu, Sigma) a_cafe <- vary_effects[, 1] b_cafe <- vary_effects[, 2] set.seed(22) N_visits <- 10 afternoon <- rep(0:1, N_visits * N_cafes / 2) cafe_id <- rep(1:N_cafes, each = N_visits) mu <- a_cafe[cafe_id] + b_cafe[cafe_id] * afternoon sigma <- 0.5 # std dev within cafes wait <- rnorm(N_visits * N_cafes, mu, sigma) d <- data.frame(cafe = cafe_id, afternoon = afternoon, wait = wait) m_cafe_separate_variation <- ulam( alist( wait ~ normal(mu, sigma), mu <- a_cafe[cafe] + b_cafe[cafe] * afternoon, a_cafe[cafe] ~ normal(a, sigma_a), b_cafe[cafe] ~ normal(b, sigma_b), a ~ normal(5, 2), b ~ normal(-1, 0.5), sigma ~ exponential(1), sigma_a ~ exponential(1), sigma_b ~ exponential(1) ), data = d, chains = 4, cores = 4, log_lik = TRUE ) m14.1 <- ulam( alist( wait ~ normal(mu, sigma), mu <- a_cafe[cafe] + b_cafe[cafe] * afternoon, c(a_cafe, b_cafe)[cafe] ~ multi_normal(c(a, b), Rho, sigma_cafe), a ~ normal(5, 2), b ~ normal(-1, 0.5), sigma_cafe ~ exponential(1), sigma ~ exponential(1), Rho ~ lkj_corr(2) ), data = d, chains = 4, cores = 4, log_lik = TRUE )  Answer. A compare plot with its associated raw data: iplot(function() { plot(compare(m14.1, m_cafe_separate_variation)) }, ar=4.5) display_markdown("Raw data (preceding plot):") display(compare(m14.1, m_cafe_separate_variation), mimetypes="text/plain")  Raw data (preceding plot):  WAIC SE dWAIC dSE pWAIC weight m14.1 303.9898 17.71746 0.000000 NA 32.61315 0.731242 m_cafe_separate_variation 305.9916 18.16013 2.001866 2.29854 32.59902 0.268758  The new model fits the training data worse (as expected, having fewer parameters). It doesnāt make up for this in the penalty term, though, because it also does less parameter sharing. 14M3. Re-estimate the varying slopes model for the UCBadmit data, now using a non-centered parameterization. Compare the efficiency of the forms of the model, using n_eff. Which is better? Which chain sampled faster? ERROR. There is no varying slopes model for UCBadmit data. This is the first chapter on varying slopes, and no model in this chapter is based on that dataset. Itās likely the non-centered parameterization does better, if we had what model the author is referring to. 14M4. Use WAIC to compare the Gaussian process model of Oceanic tools to the models fit to the same data in Chapter 11. Pay special attention to the effective numbers of parameters, as estimated by WAIC. Answer. Letās reproduce results from the two chapters: data(islandsDistMatrix) data(Kline2) # load the ordinary data, now with coordinates d <- Kline2 dP <- scale(log(d$population)) d$contact_id <- ifelse(d$contact == "high", 2, 1) dat <- list( T = d$total_tools,
P = d$P, cid = d$contact_id
)

# intercept only
m11.9 <- ulam(
alist(
T ~ dpois(lambda),
log(lambda) <- a,
a ~ dnorm(3, 0.5)
),
data = dat, chains = 4, log_lik = TRUE
)

# interaction model
m11.10 <- ulam(
alist(
T ~ dpois(lambda),
log(lambda) <- a[cid] + b[cid] * P,
a[cid] ~ dnorm(3, 0.5),
b[cid] ~ dnorm(0, 0.2)
),
data = dat, chains = 4, log_lik = TRUE
)

dat2 <- list(T = d$total_tools, P = d$population, cid = d$contact_id) m11.11 <- ulam( alist( T ~ dpois(lambda), lambda <- exp(a[cid]) * P^b[cid] / g, a[cid] ~ dnorm(1, 1), b[cid] ~ dexp(1), g ~ dexp(1) ), data = dat2, chains = 4, log_lik = TRUE ) flush.console() display(precis(m11.11, depth=2), mimetypes="text/plain") flush.console() d$society <- 1:10 # index observations
dat_list <- list(
T = d$total_tools, P = d$population,
society = d$society, Dmat = islandsDistMatrix ) m14.8 <- ulam( alist( T ~ dpois(lambda), lambda <- (a * P^b / g) * exp(k[society]), vector[10]:k ~ multi_normal(0, SIGMA), matrix[10, 10]:SIGMA <- cov_GPL2(Dmat, etasq, rhosq, 0.01), c(a, b, g) ~ dexp(1), etasq ~ dexp(2), rhosq ~ dexp(0.5) ), data = dat_list, chains = 4, cores = 4, iter = 2000, log_lik = TRUE ) flush.console() display(precis(m14.8, depth=2), mimetypes="text/plain") flush.console()  SAMPLING FOR MODEL 'ea958e18e4c604fc5bca60c02cb71b9a' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 6e-06 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.06 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.006376 seconds (Warm-up) Chain 1: 0.00597 seconds (Sampling) Chain 1: 0.012346 seconds (Total) Chain 1: SAMPLING FOR MODEL 'ea958e18e4c604fc5bca60c02cb71b9a' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 3e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 2: Iteration: 1000 / 1000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.00655 seconds (Warm-up) Chain 2: 0.006357 seconds (Sampling) Chain 2: 0.012907 seconds (Total) Chain 2: SAMPLING FOR MODEL 'ea958e18e4c604fc5bca60c02cb71b9a' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 3e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 3: Iteration: 1000 / 1000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.006347 seconds (Warm-up) Chain 3: 0.00622 seconds (Sampling) Chain 3: 0.012567 seconds (Total) Chain 3: SAMPLING FOR MODEL 'ea958e18e4c604fc5bca60c02cb71b9a' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 3e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 4: Iteration: 1000 / 1000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.006358 seconds (Warm-up) Chain 4: 0.00625 seconds (Sampling) Chain 4: 0.012608 seconds (Total) Chain 4:  SAMPLING FOR MODEL '05013cdfcc1a8f52ceb5601ee475a98d' 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.016807 seconds (Warm-up) Chain 1: 0.014315 seconds (Sampling) Chain 1: 0.031122 seconds (Total) Chain 1: SAMPLING FOR MODEL '05013cdfcc1a8f52ceb5601ee475a98d' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 6e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 2: Iteration: 1000 / 1000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.015198 seconds (Warm-up) Chain 2: 0.012859 seconds (Sampling) Chain 2: 0.028057 seconds (Total) Chain 2: SAMPLING FOR MODEL '05013cdfcc1a8f52ceb5601ee475a98d' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 5e-06 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 3: Iteration: 1000 / 1000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.017515 seconds (Warm-up) Chain 3: 0.016001 seconds (Sampling) Chain 3: 0.033516 seconds (Total) Chain 3: SAMPLING FOR MODEL '05013cdfcc1a8f52ceb5601ee475a98d' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 6e-06 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 4: Iteration: 1000 / 1000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.016075 seconds (Warm-up) Chain 4: 0.01535 seconds (Sampling) Chain 4: 0.031425 seconds (Total) Chain 4:  SAMPLING FOR MODEL '58422f20040c774e9740e486280fe76b' NOW (CHAIN 1). Chain 1: Chain 1: Gradient evaluation took 1e-05 seconds Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 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.207328 seconds (Warm-up) Chain 1: 0.159079 seconds (Sampling) Chain 1: 0.366407 seconds (Total) Chain 1: SAMPLING FOR MODEL '58422f20040c774e9740e486280fe76b' NOW (CHAIN 2). Chain 2: Chain 2: Gradient evaluation took 8e-06 seconds Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds. Chain 2: Adjust your expectations accordingly! Chain 2: Chain 2: Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 2: Iteration: 1000 / 1000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 0.179033 seconds (Warm-up) Chain 2: 0.183245 seconds (Sampling) Chain 2: 0.362278 seconds (Total) Chain 2: SAMPLING FOR MODEL '58422f20040c774e9740e486280fe76b' NOW (CHAIN 3). Chain 3: Chain 3: Gradient evaluation took 1.8e-05 seconds Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds. Chain 3: Adjust your expectations accordingly! Chain 3: Chain 3: Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 3: Iteration: 1000 / 1000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 0.20088 seconds (Warm-up) Chain 3: 0.174224 seconds (Sampling) Chain 3: 0.375104 seconds (Total) Chain 3: SAMPLING FOR MODEL '58422f20040c774e9740e486280fe76b' NOW (CHAIN 4). Chain 4: Chain 4: Gradient evaluation took 1e-05 seconds Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds. Chain 4: Adjust your expectations accordingly! Chain 4: Chain 4: Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup) Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup) Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup) Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup) Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup) Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup) Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling) Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling) Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling) Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling) Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling) Chain 4: Iteration: 1000 / 1000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 0.252573 seconds (Warm-up) Chain 4: 0.205241 seconds (Sampling) Chain 4: 0.457814 seconds (Total) Chain 4:   mean sd 5.5% 94.5% n_eff Rhat4 a[1] 0.8800983 0.65054147 -0.1834540 1.9201282 565.6799 1.002317 a[2] 0.9746582 0.86439910 -0.4510622 2.3074925 757.6732 1.003718 b[1] 0.2605599 0.03195418 0.2106948 0.3090734 1138.3743 1.001765 b[2] 0.2831560 0.10590607 0.1137215 0.4526956 622.2586 1.003477 g 1.1062560 0.70797752 0.3024785 2.3540911 629.3185 1.000265   mean sd 5.5% 94.5% n_eff Rhat4 k[1] -0.17161009 0.30414571 -0.66143815 3.048374e-01 926.5023 1.003528 k[2] -0.02551078 0.29169884 -0.47878031 4.402735e-01 753.3845 1.002713 k[3] -0.07561711 0.27855221 -0.52055453 3.612202e-01 739.4431 1.001865 k[4] 0.35279917 0.25374071 -0.02367141 7.671806e-01 847.3504 1.001366 k[5] 0.07363948 0.25145421 -0.31178870 4.764131e-01 870.9978 1.002235 k[6] -0.39117077 0.26439183 -0.82002894 8.650118e-05 981.1092 1.001846 k[7] 0.14151847 0.25239664 -0.25475289 5.382755e-01 968.6247 1.000769 k[8] -0.21500487 0.25418783 -0.62187831 1.696835e-01 897.9545 1.002621 k[9] 0.26435768 0.24242971 -0.09946002 6.444307e-01 989.3732 1.000838 k[10] -0.16908003 0.35172372 -0.73773516 3.745784e-01 1123.8908 1.001359 g 0.61245312 0.59785004 0.07536022 1.714195e+00 1780.2318 1.000548 b 0.27736210 0.08869417 0.13684643 4.203858e-01 993.4723 1.002699 a 1.41967173 1.04841388 0.24911481 3.314933e+00 2475.3185 1.000052 etasq 0.19754352 0.19657566 0.03282469 5.615212e-01 1216.1342 1.000394 rhosq 1.33530555 1.70140632 0.08355641 4.533764e+00 2029.9101 1.000207  A compare plot with its associated raw data: iplot(function() { plot(compare(m11.9, m11.10, m11.11, m14.8)) }, ar=4.5) display_markdown("Raw data (preceding plot):") display(compare(m11.9, m11.10, m11.11, m14.8), mimetypes="text/plain")  Raw data (preceding plot):  WAIC SE dWAIC dSE pWAIC weight m14.8 67.37541 2.22661 0.00000 NA 4.011871 9.980596e-01 m11.11 80.20808 11.40787 12.83267 11.55343 4.972381 1.631467e-03 m11.10 83.53601 12.13926 16.16060 12.13169 6.100473 3.089778e-04 m11.9 141.15520 31.71029 73.77979 32.71638 7.945491 9.507763e-17  The new model does better overall relative to m11.1 because it presumably has a more accurate method of estimating the contact between islands. See the first section of 14.5.1 for how the new model improves upon the old low/high contact indicator. The new model also has the lowest penalty term (effective number of parameters). Notice in the precis output above it has 10 intercept terms, and 5 other terms. This is much more than the 5 parameters in m11.1, also above. It achieves this relative to m11.1 because of the significant regularization provided by parameter sharing. 14M5. Modify the phylogenetic distance example to use group size as the outcome and brain size as a predictor. Assuming brain size influences group size, what is your estimate of the effect? How does phylogeny influence the estimate? Answer. First, letās reproduce the results from the chapter: data(Primates301) data(Primates301_nex) d <- Primates301 d$name <- as.character(d$name) dstan <- d[complete.cases(d$group_size, d$body, d$brain), ]
spp_obs <- dstan$name dat_list <- list( N_spp = nrow(dstan), M = standardize(log(dstan$body)),
B = standardize(log(dstan$brain)), G = standardize(log(dstan$group_size)),
Imat = diag(nrow(dstan))
)
m14.9 <- ulam(
alist(
B ~ multi_normal(mu, SIGMA),
mu <- a + bM * M + bG * G,
matrix[N_spp, N_spp]:SIGMA <- Imat * sigma_sq,
a ~ normal(0, 1),
c(bM, bG) ~ normal(0, 0.5),
sigma_sq ~ exponential(1)
),
data = dat_list, chains = 4, cores = 4
)
display(precis(m14.9, depth=2), mimetypes="text/plain")

library(ape)
tree_trimmed <- keep.tip(Primates301_nex, spp_obs)
Rbm <- corBrownian(phy = tree_trimmed)
V <- vcv(Rbm)

# put species in right order
dat_list$V <- V[spp_obs, spp_obs] # convert to correlation matrix dat_list$R <- dat_list$V / max(V) m14.10 <- ulam( alist( B ~ multi_normal(mu, SIGMA), mu <- a + bM * M + bG * G, matrix[N_spp, N_spp]:SIGMA <- R * sigma_sq, a ~ normal(0, 1), c(bM, bG) ~ normal(0, 0.5), sigma_sq ~ exponential(1) ), data = dat_list, chains = 4, cores = 4 ) display(precis(m14.10, depth=2), mimetypes="text/plain") Dmat <- cophenetic(tree_trimmed) # add scaled and reordered distance matrix dat_list$Dmat <- Dmat[spp_obs, spp_obs] / max(Dmat)

m14.11 <- ulam(
alist(
B ~ multi_normal(mu, SIGMA),
mu <- a + bM * M + bG * G,
matrix[N_spp, N_spp]:SIGMA <- cov_GPL1(Dmat, etasq, rhosq, 0.01),
a ~ normal(0, 1),
c(bM, bG) ~ normal(0, 0.5),
etasq ~ half_normal(1, 0.25),
rhosq ~ half_normal(3, 0.25)
),
data = dat_list, chains = 4, cores = 4
)
display(precis(m14.11, depth=2), mimetypes="text/plain")

         mean          sd          5.5%        94.5%      n_eff    Rhat4
a        -0.0001731167 0.017310086 -0.02790983 0.02609835 1939.636 0.9986765
bG        0.1230980344 0.023572481  0.08502910 0.16036884 1088.188 1.0011412
bM        0.8930430495 0.023047869  0.85791614 0.93003799 1069.554 1.0030455
sigma_sq  0.0473859234 0.005699073  0.03923401 0.05715970 1300.850 1.0026744

Warning message in Initialize.corPhyl(phy, dummy.df):
āNo covariate specified, species will be taken as ordered in the data frame. To avoid this message, specify a covariate containing the species names with the 'form' argument.ā

         mean        sd         5.5%        94.5%      n_eff    Rhat4
a        -0.18873820 0.16496198 -0.44557098 0.07618946 2216.384 1.0027970
bG       -0.01217738 0.01971775 -0.04277674 0.01980523 2350.256 0.9990575
bM        0.70105765 0.03861462  0.63972330 0.76345226 1776.208 1.0010190
sigma_sq  0.16148768 0.01875378  0.13339710 0.19432609 2409.753 1.0005959

      mean        sd          5.5%        94.5%      n_eff    Rhat4
a     -0.06587284 0.076292599 -0.18646383 0.05812055 2012.089 0.9999619
bG     0.04981556 0.023797105  0.01034797 0.08758256 2213.185 1.0003827
bM     0.83321792 0.030898044  0.78324467 0.88101712 2727.300 0.9988425
etasq  0.03494790 0.006854313  0.02494174 0.04695430 1937.816 0.9993151
rhosq  2.79131896 0.250754096  2.39426624 3.20118552 2170.700 0.9996398


Reversing the prediction problem:

dat_list <- list(
N_spp = nrow(dstan),
M = standardize(log(dstan$body)), B = standardize(log(dstan$brain)),
G = standardize(log(dstan$group_size)), Imat = diag(nrow(dstan)) ) m14.9_reversed <- ulam( alist( G ~ multi_normal(mu, SIGMA), mu <- a + bM * M + bB * B, matrix[N_spp, N_spp]:SIGMA <- Imat * sigma_sq, a ~ normal(0, 1), c(bM, bB) ~ normal(0, 0.5), sigma_sq ~ exponential(1) ), data = dat_list, chains = 4, cores = 4 ) display(precis(m14.9_reversed, depth=2), mimetypes="text/plain") # put species in right order dat_list$V <- V[spp_obs, spp_obs]
# convert to correlation matrix
dat_list$R <- dat_list$V / max(V)

m14.10_reversed <- ulam(
alist(
G ~ multi_normal(mu, SIGMA),
mu <- a + bM * M + bB * B,
matrix[N_spp, N_spp]:SIGMA <- R * sigma_sq,
a ~ normal(0, 1),
c(bM, bB) ~ normal(0, 0.5),
sigma_sq ~ exponential(1)
),
data = dat_list, chains = 4, cores = 4
)
display(precis(m14.10_reversed, depth=2), mimetypes="text/plain")

# add scaled and reordered distance matrix
dat_list$Dmat <- Dmat[spp_obs, spp_obs] / max(Dmat) m14.11_reversed <- ulam( alist( G ~ multi_normal(mu, SIGMA), mu <- a + bM * M + bB * B, matrix[N_spp, N_spp]:SIGMA <- cov_GPL1(Dmat, etasq, rhosq, 0.01), a ~ normal(0, 1), c(bM, bB) ~ normal(0, 0.5), etasq ~ half_normal(1, 0.25), rhosq ~ half_normal(3, 0.25) ), data = dat_list, chains = 4, cores = 4 ) display(precis(m14.11_reversed, depth=2), mimetypes="text/plain")   mean sd 5.5% 94.5% n_eff Rhat4 a -0.001417642 0.05685219 -0.09147605 0.08793679 1327.7488 0.9992159 bB 1.004936856 0.20442226 0.68104136 1.34138179 854.5307 1.0005631 bM -0.337334331 0.20456292 -0.67326748 -0.01079579 828.9535 1.0016998 sigma_sq 0.519471576 0.06055068 0.42832399 0.61695123 1429.3501 1.0018326   mean sd 5.5% 94.5% n_eff Rhat4 a -0.47774105 0.5678362 -1.359264822 0.4414133 1496.797 1.000575 bB -0.07150671 0.2583164 -0.476816067 0.3477977 1091.952 1.000832 bM 0.34676958 0.2160707 0.001476129 0.6938241 1003.165 1.001472 sigma_sq 2.68982107 0.3120843 2.234004692 3.2161848 1479.938 1.003012   mean sd 5.5% 94.5% n_eff Rhat4 a -0.5021950 0.3442075 -1.0440793 0.03602275 1511.5889 1.000513 bB 0.1850170 0.2626387 -0.2375381 0.60738516 968.2030 1.000422 bM 0.1880509 0.2259942 -0.1873261 0.54406083 989.1999 1.000826 etasq 0.9309690 0.1206105 0.7537561 1.14232420 1593.3978 1.000191 rhosq 3.0219251 0.2384142 2.6610231 3.40037960 1592.1098 1.003639  Interestingly, the ordinary regression model makes an even more extreme prediction about how group size influences brain size than the reverse. Still, the uncertainties are large. The addition of phylogeny once again removes the presumed effect. The uncertainties are so large in these latter two models that they really say almost nothing about an association. 14H1. Letās revisit the Bangladesh fertility data, data(bangladesh), from the practice problems for Chapter 13. Fit a model with both varying intercepts by district_id and varying slopes of urban by district_id. You are still predicting use.contraception. Inspect the correlation between the intercepts and slopes. Can you interpret this correlation, in terms of what it tells you about the pattern of contraceptive use in the sample? It might help to plot the mean (or median) varying effect estimates for both the intercepts and slopes, by district. Then you can visualize the correlation and maybe more easily think through what it means to have a particular correlation. Plotting predicted proportion of women using contraception, with urban women on one axis and rural on the other, might also help. Answer. To review, the help for the bangladesh data.frame: data(bangladesh) display(help(bangladesh)) bc_df <- bangladesh bc_df$district_id <- as.integer(as.factor(bc_df$district))   bangladesh {rethinking} R Documentation ## Bangladesh contraceptive use data ### Description Contraceptive use data from 1934 Bangladeshi women. ### Usage data(bangladesh)  ### Format 1. woman : ID number for each woman in sample 2. district : Number for each district 3. use.contraception : 0/1 indicator of contraceptive use 4. living.children : Number of living children 5. age.centered : Centered age 6. urban : 0/1 indicator of urban context ### References Bangladesh Fertility Survey, 1989 [Package rethinking version 2.13 ] A head and summary of the bangladesh data.frame, with the new variable suggested by the author in question 13H1: display(head(bc_df)) display(summary(bc_df))  A data.frame: 6 Ć 7 womandistrictuse.contraceptionliving.childrenage.centeredurbandistrict_id <int><int><int><int><dbl><int><int> 11104 18.440011 22101 -5.559911 33103 1.440011 44104 8.440011 55101-13.559011 66101-11.560011  woman district use.contraception living.children Min. : 1.0 Min. : 1.00 Min. :0.0000 Min. :1.000 1st Qu.: 484.2 1st Qu.:14.00 1st Qu.:0.0000 1st Qu.:1.000 Median : 967.5 Median :29.00 Median :0.0000 Median :3.000 Mean : 967.5 Mean :29.35 Mean :0.3925 Mean :2.652 3rd Qu.:1450.8 3rd Qu.:45.00 3rd Qu.:1.0000 3rd Qu.:4.000 Max. :1934.0 Max. :61.00 Max. :1.0000 Max. :4.000 age.centered urban district_id Min. :-13.560000 Min. :0.0000 Min. : 1.00 1st Qu.: -7.559900 1st Qu.:0.0000 1st Qu.:14.00 Median : -1.559900 Median :0.0000 Median :29.00 Mean : 0.002198 Mean :0.2906 Mean :29.25 3rd Qu.: 6.440000 3rd Qu.:1.0000 3rd Qu.:45.00 Max. : 19.440000 Max. :1.0000 Max. :60.00  Sampling from the varying intercepts model: # If you want to reproduce results from the previous chapter: # # bc_dat <- list( # UseContraception = bc_df$use.contraception,
#   DistrictId = bc_df$district_id # ) # m_bc_ve <- ulam( # alist( # UseContraception ~ dbinom(1, p), # logit(p) <- a[DistrictId], # a[DistrictId] ~ dnorm(a_bar, sigma), # a_bar ~ dnorm(0, 1.5), # sigma ~ dexp(1) # ), # data = bc_dat, chains = 4, cores = 4, log_lik = TRUE # ) bc_dat <- list( UseContraception = bc_df$use.contraception,
DistrictId = bc_df$district_id, Urban = bc_df$urban
)
m_bc_vis <- ulam(
alist(
UseContraception ~ dbinom(1, p),
logit(p) <- a_district[DistrictId] + b_district[DistrictId] * Urban,
c(a_district, b_district)[DistrictId] ~ multi_normal(c(a, b), Rho, sigma_intercepts_slopes),
a ~ normal(0, 2),
b ~ normal(0, 0.5),
sigma_intercepts_slopes ~ exponential(1),
Rho ~ lkj_corr(2)
),
data = bc_dat, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m_bc_vis, depth=3), mimetypes="text/plain")
iplot(function() {
plot(precis(m_bc_vis, depth=3), main="m_bc_vis")
}, ar=0.4)

post <- extract.samples(m_bc_vis) # posterior
R <- rlkjcorr(1e4, K = 2, eta = 2) # prior

iplot(function() {
dens(post$Rho[, 1, 2], xlim = c(-1, 1)) dens(R[, 1, 2], add = TRUE, lty = 2) })  Warning message: āBulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#bulk-essā  Warning message: āTail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#tail-essā   mean sd 5.5% 94.5% b_district[1] 1.0670860 0.3782827 0.4859153953 1.6648881 b_district[2] 0.5967072 0.6645532 -0.4330004363 1.6322922 b_district[3] 0.8848223 0.8057362 -0.3856858371 2.2487559 b_district[4] 1.5814457 0.6347096 0.6560098348 2.6659842 b_district[5] 0.6113797 0.5942621 -0.3188659912 1.5578244 b_district[6] 1.2810511 0.5527699 0.4667433362 2.2225933 b_district[7] 0.7745289 0.6689434 -0.2951438126 1.8745062 b_district[8] 0.9065896 0.6014774 0.0002929298 1.8780210 b_district[9] 0.9863604 0.6166353 0.0425672943 1.9837517 b_district[10] 1.1954027 0.7592615 0.0276866328 2.4434961 b_district[11] 1.5756218 0.8116946 0.4444062308 2.9595681 b_district[12] 0.4322015 0.5665931 -0.4999139671 1.3354042 b_district[13] 0.2684581 0.5874672 -0.6925710533 1.2110490 b_district[14] 1.2110843 0.4319248 0.5424204957 1.9125391 b_district[15] 0.4493872 0.5861374 -0.4998182743 1.3609454 b_district[16] 0.5669155 0.6321994 -0.4256460915 1.5818233 b_district[17] 0.7624304 0.6745130 -0.3045817881 1.8874580 b_district[18] 0.8479431 0.4911165 0.0857845437 1.6496457 b_district[19] 0.9680160 0.6254921 0.0303852148 2.0321759 b_district[20] 0.5403988 0.6952316 -0.5619983408 1.6191578 b_district[21] -0.3209662 0.6859289 -1.4403290374 0.7075015 b_district[22] 0.9749350 0.6772639 -0.0584712608 2.1394007 b_district[23] 0.8028469 0.6992094 -0.2879172061 1.9132030 b_district[24] 1.2329838 0.7372550 0.1094543016 2.4527721 b_district[25] 0.1879679 0.4248847 -0.5040083386 0.8567129 b_district[26] 0.5668825 0.7233731 -0.5650370280 1.7357724 b_district[27] 1.1476359 0.6180868 0.2054302074 2.1557858 b_district[28] 0.7334328 0.6005227 -0.2234647593 1.6730947 b_district[29] 1.1410977 0.5579254 0.2976138234 2.0465502 b_district[30] 1.0106844 0.5002790 0.2772004374 1.8338426 ā® ā® ā® ā® ā® a_district[39] -0.24878952 3.446544e-01 -0.8001652 0.30037495 a_district[40] -0.53367224 4.056508e-01 -1.1726077 0.11575554 a_district[41] -0.08785561 3.647116e-01 -0.6613238 0.48234779 a_district[42] -0.05918072 5.055008e-01 -0.8316718 0.77688896 a_district[43] -0.22036656 3.018605e-01 -0.6984655 0.24937437 a_district[44] -1.05308364 3.537547e-01 -1.6495761 -0.50317067 a_district[45] -0.91672877 3.193301e-01 -1.4355843 -0.41527769 a_district[46] -0.09571376 2.228935e-01 -0.4450982 0.26582506 a_district[47] -0.51778990 4.465133e-01 -1.2312732 0.21138532 a_district[48] -0.22569530 3.243286e-01 -0.7331933 0.30002079 a_district[49] -1.04824482 5.298423e-01 -1.9150630 -0.27016833 a_district[50] -0.67188490 4.024037e-01 -1.3064462 -0.02679504 a_district[51] -0.64603156 3.617746e-01 -1.2429197 -0.06510586 a_district[52] -0.11721394 2.784987e-01 -0.5601760 0.32948042 a_district[53] -0.70873036 5.682276e-01 -1.5958060 0.20074087 a_district[54] -0.75671498 6.164120e-01 -1.7673297 0.23632531 a_district[55] -0.10909625 3.251141e-01 -0.6168951 0.42430576 a_district[56] -1.17789470 3.753756e-01 -1.7906242 -0.60817455 a_district[57] -0.17468722 3.562620e-01 -0.7293605 0.39079686 a_district[58] -1.17842367 4.955465e-01 -2.0236499 -0.42535810 a_district[59] -1.15104807 3.848318e-01 -1.7852631 -0.56316102 a_district[60] -1.17134926 3.309090e-01 -1.7302584 -0.68202729 a -0.68697115 9.680213e-02 -0.8431896 -0.53422309 b 0.64241389 1.533845e-01 0.3911194 0.89095418 sigma_intercepts_slopes[1] 0.57516450 9.473592e-02 0.4356859 0.73634803 sigma_intercepts_slopes[2] 0.76059552 1.976562e-01 0.4589757 1.08665688 Rho[1,1] 1.00000000 0.000000e+00 1.0000000 1.00000000 Rho[1,2] -0.65363331 1.671887e-01 -0.8613258 -0.34409840 Rho[2,1] -0.65363331 1.671887e-01 -0.8613258 -0.34409840 Rho[2,2] 1.00000000 6.272117e-17 1.0000000 1.00000000 n_eff Rhat4 b_district[1] 1831.8290 1.0012389 b_district[2] 2883.8034 0.9994656 b_district[3] 2096.1790 0.9996167 b_district[4] 742.8528 1.0039141 b_district[5] 2158.3104 0.9983583 b_district[6] 1916.0610 0.9987753 b_district[7] 2448.2445 0.9993163 b_district[8] 1779.2636 1.0011896 b_district[9] 2200.3913 0.9991816 b_district[10] 1714.6876 1.0003149 b_district[11] 1021.2957 1.0014837 b_district[12] 3222.1393 1.0006625 b_district[13] 2065.5455 0.9998587 b_district[14] 1408.5814 1.0017142 b_district[15] 2473.9455 0.9988708 b_district[16] 1695.4276 1.0019869 b_district[17] 2690.2394 0.9987431 b_district[18] 2195.1137 0.9999696 b_district[19] 1638.6104 1.0022221 b_district[20] 2519.0662 0.9989753 b_district[21] 921.4559 1.0005069 b_district[22] 2082.2386 1.0003840 b_district[23] 2111.5169 1.0019826 b_district[24] 1380.2483 1.0008525 b_district[25] 1650.5842 1.0014639 b_district[26] 2522.6748 1.0012108 b_district[27] 1632.7972 1.0001834 b_district[28] 3012.6085 0.9992483 b_district[29] 1858.0194 0.9996650 b_district[30] 1301.2451 1.0008541 ā® ā® ā® a_district[39] 2481.3080 1.0003322 a_district[40] 1950.4880 1.0025755 a_district[41] 2499.2493 0.9988345 a_district[42] 1573.8640 1.0001231 a_district[43] 2089.1375 0.9988561 a_district[44] 2712.3422 0.9987354 a_district[45] 2734.5355 1.0001725 a_district[46] 2624.5404 0.9995971 a_district[47] 2193.7884 0.9993473 a_district[48] 1930.0847 0.9995102 a_district[49] 2290.7133 0.9989410 a_district[50] 2363.4590 0.9991670 a_district[51] 2288.5704 0.9997117 a_district[52] 1992.2301 0.9997977 a_district[53] 2065.2941 0.9997024 a_district[54] 1532.5275 0.9999255 a_district[55] 2246.6413 0.9994746 a_district[56] 2510.3441 0.9997795 a_district[57] 1348.0636 0.9997778 a_district[58] 2411.5027 0.9985544 a_district[59] 2692.0315 0.9991757 a_district[60] 2315.3140 0.9989207 a 1699.3030 0.9988051 b 1236.1160 1.0003488 sigma_intercepts_slopes[1] 604.1727 1.0028934 sigma_intercepts_slopes[2] 263.0453 1.0080110 Rho[1,1] NaN NaN Rho[1,2] 379.9839 1.0097282 Rho[2,1] 379.9839 1.0097282 Rho[2,2] 1937.6490 0.9979980  Not surprisingly, the b slope parameter and most of the b_district are positive. Normally, living in an urban area makes one more open to ideas like contraception. The a parameters are likely negative because in general many women do not use contraception. The correlation is negative, because a is typically negative and b is typically positive. It is also relatively large in absolute magnitude, implying that a more negative intercept is associated with a more positive slope (and vice versa). We can see this in a plot of mean intercepts vs. mean slopes: iplot(function() { post <- extract.samples(m_bc_vis) a2 <- apply(post$a_district, 2, mean)
b2 <- apply(post$b_district, 2, mean) plot(a2, b2, xlab = "intercept", ylab = "slope", pch = 1, col = rangi2, ylim = c(min(b2) - 0.1, max(b2) + 0.1), xlim = c(min(a2) - 0.1, max(a2) + 0.1) ) Mu_est <- c(mean(post$a), mean(post$b)) rho_est <- mean(post$Rho[, 1, 2])
sa_est <- mean(post$sigma_intercepts_slopes[, 1]) sb_est <- mean(post$sigma_intercepts_slopes[, 2])
cov_ab <- sa_est * sb_est * rho_est
Sigma_est <- matrix(c(sa_est^2, cov_ab, cov_ab, sb_est^2), ncol = 2)
library(ellipse)
for (l in c(0.1, 0.3, 0.5, 0.8, 0.99)) {
lines(ellipse(Sigma_est, centre = Mu_est, level = l),
col = col.alpha("black", 0.2)
)
}
})


This is a similar relationship to what we saw in the cafe example at the start of the chapter; more extreme intercepts (cafes, districts) are associated with more extreme slopes. In this case the districts with especially low contraceptive use are associated with a greater increase in use when a woman āmovesā to an urban area within the district.

14H2. Now consider the predictor variables age.centered and living.children, also contained in data(bangladesh). Suppose that age influences contraceptive use (changing attitudes) and number of children (older people have had more time to have kids). Number of children may also directly influence contraceptive use. Draw a DAG that reflects these hypothetical relationships. Then build models needed to evaluate the DAG. You will need at least two models. Retain district and urban, as in 14H1. What do you conclude about the causal influence of age and children?

Consider the following causal diagram. Is this diagram reasonable? The author asked us to only consider this one DAG, but itās worth thinking whether itās reasonable so we know what to expect from our inferences.

The Age variable only has arrows pointing out of it. As discussed elsewhere (see question 12H7), this is the only way we should put time-based variables on DAGs. We can think about Age as describing the timing of someoneās birth. In terms of a [counterfactual conditional][cc], we would say that if a woman had been born e.g. 10 years earlier she would not have used contraception. A counterfactual conditional like this describes a causal theory, but expressed regarding the past. Casual theories should still be tested, when possible, by actively controlling a variable in the future. In this case we wonāt be able to test the theory in the future (as well as the past) because we canāt control the date of anyoneās birth, unless we ran an experiment with identical twins where we somehow froze the zygote of one twin in time for years. Weāll include a description of all existing and potential arrows in this DAG as counterfactual conditionals below.

Itās reasonable to think that Age could affect Urban. More people are moving to cities every year across the world, for jobs, and so in general weād expect younger women to be living in cities. That is, if a woman was born later she would be more likely to be urban.

Itās possible that Age could affect District if there was a people migration happening across the country at one point in the past and slightly older women were less likely to move. That is, if a woman was born later she could be living in a different district.

Some districts could be more urbanized and therefore District could predict Urban, and vice versa. That is, if a woman was living in a different district she could be more likely to be urban. If a woman was not in an urban area she would be more likely to be in a different district.

Some districts might have governmental programs or tax incentives to encourage or discourage children, implying a relationship between District and Children mediated by the program. This would confound our causal inference about the effect of Children on UseContraception if we didnāt include District in our model. Similarly, for urban programs encouraging or discouraging children. If a woman was living in a different district, she may have had more or fewer children.

Itās also likely that women move to more rural areas for e.g. cheaper housing when they have more children, or to districts with cheaper housing and costs of living. That is, a woman had fewer children, she may be living in a less urban area or in general somewhere else.

Like the controversy over global warming in the United States, there are often many possible confounds we can suggest to make a DAG more complicated. It can be hard to decide which need to be included in every model without getting into a lot of details.

In this question we are predicting UseContraception from Children but weād typically think of the opposite causal path: contraception clearly influences the number of children a woman will have. In this case, we are essentially assuming that the number of children a woman has in the present influences her present decision to use contraception. In general, though, these variables interact, so we have to look at the casual paths similar as a causal time series similar to how the author treated the causal influence of group size on brain size at the beginning of section 14.5.2. Our causal inferences here will only apply to the point in time when these surveys were taken.

Weāll ignore all of this. The DAG expected by the question:

library(dagitty)

expected_dag <- dagitty('
dag {
bb="0,0,1,1"
Age [pos="0.3,0.1"]
Children [exposure,pos="0.5,0.1"]
District [pos="0.65,0.2"]
Urban [pos="0.65,0.1"]
UseContraception [outcome,pos="0.4,0.2"]
Age -> Children
Age -> UseContraception
Children -> UseContraception
District -> UseContraception
Urban -> UseContraception
}')
iplot(function() plot(expected_dag), scale=10)

Attaching package: ādagittyā

The following object is masked from āpackage:apeā:

edges


Implied conditional independencies:

display(impliedConditionalIndependencies(expected_dag))

Age _||_ Dstr
Age _||_ Urbn
Chld _||_ Dstr
Chld _||_ Urbn
Dstr _||_ Urbn


We will build two models, as the question suggests. The first will include both predictors, so we can infer the direct/total causal effect of children on UseContraception:

display(adjustmentSets(expected_dag, exposure="Children", outcome="UseContraception", effect="direct"))

{ Age }

{ Age }


The first model will also predict the direct effect of Age on UseContraception:

display(adjustmentSets(expected_dag, exposure="Age", outcome="UseContraception", effect="direct"))

{ Children }


The second will include only Age so we can infer the total causal effect of Age on UseContraception:

display(adjustmentSets(expected_dag, exposure="Age", outcome="UseContraception", effect="total"))

 {}


Fitting the first model:

data(bangladesh)
bc_df$district_id <- as.integer(as.factor(bc_df$district))

bc_dat <- list(
UseContraception = bc_df$use.contraception, DistrictId = bc_df$district_id,
Urban = bc_df$urban, Age = standardize(bc_df$age.centered),
Children = bc_df$living.children ) m_bc_age_children <- ulam( alist( UseContraception ~ dbinom(1, p), logit(p) <- a_district[DistrictId] + b_district[DistrictId] * Urban + bAge * Age + a_children[Children], c(a_district, b_district)[DistrictId] ~ multi_normal(c(a, b), Rho, sigma_intercepts_slopes), a_children[Children] ~ normal(0, 1), a ~ normal(0, 2), b ~ normal(0, 0.5), bAge ~ normal(0, 0.5), sigma_intercepts_slopes ~ exponential(1), Rho ~ lkj_corr(2) ), data = bc_dat, chains = 4, cores = 4, log_lik = TRUE ) display(precis(m_bc_age_children, depth=3), mimetypes="text/plain") iplot(function() { plot(precis(m_bc_age_children, depth=3), main="m_bc_age_children") }, ar=0.4)  Warning message: āThe largest R-hat is NA, indicating chains have not mixed. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#r-hatā  Warning message: āBulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#bulk-essā  Warning message: āTail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#tail-essā   mean sd 5.5% 94.5% b_district[1] 1.0999167 0.3908649 0.49159563 1.7417916 b_district[2] 0.7624000 0.6713341 -0.29147820 1.8229521 b_district[3] 1.0043785 0.8060123 -0.20528392 2.3405111 b_district[4] 1.6438375 0.6611769 0.66845155 2.7468812 b_district[5] 0.7022189 0.5919901 -0.23480486 1.6426724 b_district[6] 1.4202445 0.5546021 0.60059064 2.3576044 b_district[7] 0.9034358 0.7198666 -0.18233987 2.0596096 b_district[8] 1.0075832 0.6370145 0.03735138 2.0919893 b_district[9] 1.1294339 0.6425540 0.14775703 2.1998940 b_district[10] 1.1811700 0.7393612 0.05655746 2.3937003 b_district[11] 1.5655565 0.8298174 0.33789497 2.9904372 b_district[12] 0.4899962 0.5848951 -0.46789684 1.3872553 b_district[13] 0.3906549 0.5638410 -0.54148595 1.2365869 b_district[14] 1.4123033 0.4611426 0.69587260 2.1682716 b_district[15] 0.4825745 0.6121071 -0.52356196 1.4521823 b_district[16] 0.5243952 0.6783461 -0.53337786 1.6515876 b_district[17] 0.8794557 0.7167730 -0.23753859 2.0253755 b_district[18] 1.0268427 0.5122685 0.21549225 1.8279752 b_district[19] 1.0270314 0.6237950 0.06036288 2.0466071 b_district[20] 0.5158529 0.7126419 -0.66672066 1.6332599 b_district[21] -0.3279448 0.7220129 -1.58727814 0.7005328 b_district[22] 1.1335031 0.7273222 -0.00905651 2.3106420 b_district[23] 0.9190812 0.7432049 -0.20580069 2.0794233 b_district[24] 1.2563269 0.7596747 0.09519452 2.4983379 b_district[25] 0.4333071 0.4399538 -0.26694340 1.1220195 b_district[26] 0.7074407 0.7235723 -0.42221307 1.8767707 b_district[27] 1.1893079 0.5969365 0.25626517 2.1375584 b_district[28] 0.7021765 0.6034404 -0.30269080 1.6450539 b_district[29] 1.2046125 0.5733478 0.32493704 2.1606517 b_district[30] 0.9428385 0.4741211 0.22253355 1.7132574 ā® ā® ā® ā® ā® a_district[44] -1.00757529 6.136324e-01 -1.9911938 -0.02934434 a_district[45] -0.97830892 5.980057e-01 -1.9538494 -0.02719440 a_district[46] -0.03314476 5.566992e-01 -0.9117773 0.86163084 a_district[47] -0.38205712 6.907393e-01 -1.4876866 0.73290678 a_district[48] -0.10582759 6.133921e-01 -1.0728309 0.85300833 a_district[49] -0.94369673 7.508873e-01 -2.1100061 0.23350557 a_district[50] -0.51192605 6.692862e-01 -1.5518202 0.56943445 a_district[51] -0.67770025 6.343153e-01 -1.6873408 0.33168960 a_district[52] -0.06297412 5.825162e-01 -0.9578809 0.84667572 a_district[53] -0.66126462 8.167981e-01 -1.9708303 0.71126875 a_district[54] -0.70715221 8.193827e-01 -2.0321862 0.58303659 a_district[55] 0.04195840 6.380596e-01 -0.9442829 1.07382583 a_district[56] -1.11721014 6.541691e-01 -2.1445767 -0.07994060 a_district[57] -0.01447033 6.520955e-01 -1.0314944 1.03048964 a_district[58] -1.21266444 6.882528e-01 -2.2795963 -0.09667581 a_district[59] -1.08185653 6.427204e-01 -2.0945485 -0.07955301 a_district[60] -1.21806642 6.201448e-01 -2.2460559 -0.27990187 a_children[1] -1.02952192 5.207069e-01 -1.8500820 -0.19564254 a_children[2] 0.08108856 5.209392e-01 -0.7603851 0.89890579 a_children[3] 0.31428682 5.225278e-01 -0.5115419 1.14501164 a_children[4] 0.29755674 5.199256e-01 -0.5168079 1.12623139 a -0.63573454 5.211178e-01 -1.4590586 0.18399513 b 0.71397519 1.589268e-01 0.4585300 0.96581740 bAge -0.22737831 6.922955e-02 -0.3359584 -0.11826217 sigma_intercepts_slopes[1] 0.60466470 1.017295e-01 0.4521114 0.77621481 sigma_intercepts_slopes[2] 0.77828611 2.040382e-01 0.4604246 1.12411203 Rho[1,1] 1.00000000 0.000000e+00 1.0000000 1.00000000 Rho[1,2] -0.64507055 1.684769e-01 -0.8600055 -0.32967889 Rho[2,1] -0.64507055 1.684769e-01 -0.8600055 -0.32967889 Rho[2,2] 1.00000000 6.153016e-17 1.0000000 1.00000000 n_eff Rhat4 b_district[1] 1839.5571 0.9984079 b_district[2] 3409.2984 0.9999534 b_district[3] 1706.9000 0.9987609 b_district[4] 669.9120 1.0055624 b_district[5] 3190.8290 1.0003951 b_district[6] 1044.2280 1.0018416 b_district[7] 2599.6768 1.0010587 b_district[8] 1737.1760 0.9997373 b_district[9] 2613.2646 0.9985529 b_district[10] 1507.0233 1.0003189 b_district[11] 1326.2529 1.0000899 b_district[12] 1944.0288 0.9995990 b_district[13] 2282.0064 0.9993529 b_district[14] 1403.3040 1.0000731 b_district[15] 2218.0897 1.0004139 b_district[16] 2212.6778 1.0010683 b_district[17] 2962.5253 0.9990428 b_district[18] 2344.2862 0.9990850 b_district[19] 2133.9400 0.9998191 b_district[20] 2473.1233 0.9997963 b_district[21] 656.1066 1.0051644 b_district[22] 1823.7037 0.9991960 b_district[23] 2581.8564 0.9991919 b_district[24] 1221.3117 1.0003195 b_district[25] 2773.4273 1.0004201 b_district[26] 2760.2540 0.9990907 b_district[27] 1834.5880 0.9995749 b_district[28] 1917.9901 0.9994573 b_district[29] 1352.9754 0.9996811 b_district[30] 2076.7044 1.0010753 ā® ā® ā® a_district[44] 63.01853 1.0496414 a_district[45] 58.62364 1.0528270 a_district[46] 53.68775 1.0638962 a_district[47] 78.95846 1.0383481 a_district[48] 62.06339 1.0590729 a_district[49] 99.78424 1.0320701 a_district[50] 72.83560 1.0464830 a_district[51] 73.44679 1.0488270 a_district[52] 55.62306 1.0669655 a_district[53] 101.89727 1.0326683 a_district[54] 109.91239 1.0319721 a_district[55] 61.94630 1.0530400 a_district[56] 59.75073 1.0507623 a_district[57] 61.80999 1.0585249 a_district[58] 75.79770 1.0377480 a_district[59] 73.98629 1.0443772 a_district[60] 67.82203 1.0530765 a_children[1] 44.84983 1.0711356 a_children[2] 43.03068 1.0756245 a_children[3] 44.55434 1.0768130 a_children[4] 42.62129 1.0777767 a 42.14154 1.0779772 b 1048.78662 1.0006900 bAge 2819.64745 0.9999493 sigma_intercepts_slopes[1] 695.84499 1.0051581 sigma_intercepts_slopes[2] 210.86961 1.0173553 Rho[1,1] NaN NaN Rho[1,2] 489.22550 1.0078223 Rho[2,1] 489.22550 1.0078223 Rho[2,2] 1753.63666 0.9979980  Fitting the second model: bc_dat <- list( UseContraception = bc_df$use.contraception,
DistrictId = bc_df$district_id, Urban = bc_df$urban,
Age = standardize(bc_df$age.centered) ) m_bc_age <- ulam( alist( UseContraception ~ dbinom(1, p), logit(p) <- a_district[DistrictId] + b_district[DistrictId] * Urban + bAge * Age, c(a_district, b_district)[DistrictId] ~ multi_normal(c(a, b), Rho, sigma_intercepts_slopes), a ~ normal(0, 2), b ~ normal(0, 0.5), bAge ~ normal(0, 0.5), sigma_intercepts_slopes ~ exponential(1), Rho ~ lkj_corr(2) ), data = bc_dat, chains = 4, cores = 4, log_lik = TRUE ) display(precis(m_bc_age, depth=3), mimetypes="text/plain") iplot(function() { plot(precis(m_bc_age, depth=3), main="m_bc_age") }, ar=0.4)  Warning message: āBulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#bulk-essā  Warning message: āTail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#tail-essā   mean sd 5.5% 94.5% b_district[1] 1.0702136 0.4067277 0.43379121 1.7248416 b_district[2] 0.6305582 0.6912087 -0.47099620 1.6988042 b_district[3] 0.9705495 0.7902006 -0.25163575 2.2485909 b_district[4] 1.6429103 0.6505227 0.66302723 2.7696751 b_district[5] 0.6356402 0.6139752 -0.33570188 1.6054755 b_district[6] 1.2899324 0.5596966 0.46314778 2.2128706 b_district[7] 0.7735038 0.7020769 -0.33740407 1.8648556 b_district[8] 0.9321451 0.6491707 -0.06116912 1.9916973 b_district[9] 0.9917004 0.6714883 -0.03799627 2.0999914 b_district[10] 1.2438967 0.7840562 0.07436494 2.5460826 b_district[11] 1.6047628 0.8444329 0.31912278 2.9721018 b_district[12] 0.4267246 0.5973661 -0.52816545 1.3816863 b_district[13] 0.2860661 0.5781811 -0.67536741 1.1610802 b_district[14] 1.2535781 0.4232179 0.61710183 1.9822707 b_district[15] 0.4207730 0.5816601 -0.51392978 1.3405875 b_district[16] 0.5996888 0.6929692 -0.45002673 1.7215217 b_district[17] 0.7783539 0.7053412 -0.30324018 1.9164783 b_district[18] 0.8717355 0.4831137 0.12611949 1.6636872 b_district[19] 0.9747263 0.6338828 0.04406323 2.0289520 b_district[20] 0.5074256 0.7368209 -0.67157506 1.6615544 b_district[21] -0.3972980 0.6662763 -1.51749499 0.6102074 b_district[22] 0.9984948 0.7441407 -0.10027794 2.1627915 b_district[23] 0.8053357 0.7483680 -0.35752667 1.9867326 b_district[24] 1.2397900 0.7926735 0.07207726 2.6201476 b_district[25] 0.1991053 0.4276917 -0.49510261 0.8834087 b_district[26] 0.5621284 0.7225624 -0.58839125 1.6748054 b_district[27] 1.1391534 0.5962353 0.18383570 2.0653965 b_district[28] 0.7071302 0.5805124 -0.21788665 1.6359116 b_district[29] 1.1159222 0.5861940 0.20646296 2.0684271 b_district[30] 1.0274215 0.5023722 0.28346471 1.8659726 ā® ā® ā® ā® ā® a_district[40] -0.56713225 3.942035e-01 -1.196080694 0.05692364 a_district[41] -0.09445060 3.518192e-01 -0.646965536 0.47154724 a_district[42] -0.07645624 4.785484e-01 -0.819215319 0.67868498 a_district[43] -0.20928463 3.245208e-01 -0.706348092 0.30756017 a_district[44] -1.05076852 3.564967e-01 -1.631410062 -0.51042323 a_district[45] -0.90137965 3.195985e-01 -1.430802543 -0.39026654 a_district[46] -0.09572654 2.242291e-01 -0.458310570 0.26297881 a_district[47] -0.52708188 4.322266e-01 -1.204760880 0.15630742 a_district[48] -0.23006057 3.313788e-01 -0.745034819 0.29127181 a_district[49] -1.03023255 5.507298e-01 -1.939590161 -0.17612258 a_district[50] -0.68555663 3.984149e-01 -1.333323326 -0.04795596 a_district[51] -0.63496193 3.707027e-01 -1.231812436 -0.04323402 a_district[52] -0.11045522 2.763375e-01 -0.548731329 0.32194903 a_district[53] -0.71767053 5.517730e-01 -1.635052142 0.15193477 a_district[54] -0.73863665 5.836698e-01 -1.661464909 0.14930034 a_district[55] -0.12110828 3.453037e-01 -0.655014058 0.43152436 a_district[56] -1.18951735 3.893304e-01 -1.803571667 -0.59085226 a_district[57] -0.17772950 3.702615e-01 -0.786031845 0.39128593 a_district[58] -1.15459443 4.728163e-01 -1.948179197 -0.45065063 a_district[59] -1.15298373 3.783223e-01 -1.753136917 -0.57052335 a_district[60] -1.18366333 3.435825e-01 -1.755721227 -0.66582234 a -0.69088692 9.915224e-02 -0.847275102 -0.53470701 b 0.64328244 1.582692e-01 0.390786598 0.89773685 bAge 0.08306864 4.881746e-02 0.006175679 0.16085576 sigma_intercepts_slopes[1] 0.57408849 9.799291e-02 0.424628883 0.73759442 sigma_intercepts_slopes[2] 0.78539789 1.963405e-01 0.483581293 1.10263891 Rho[1,1] 1.00000000 0.000000e+00 1.000000000 1.00000000 Rho[1,2] -0.64207947 1.638724e-01 -0.854575139 -0.34237655 Rho[2,1] -0.64207947 1.638724e-01 -0.854575139 -0.34237655 Rho[2,2] 1.00000000 6.192971e-17 1.000000000 1.00000000 n_eff Rhat4 b_district[1] 1829.9801 1.0010016 b_district[2] 2009.4880 1.0004020 b_district[3] 1741.9526 0.9989266 b_district[4] 743.9184 0.9995581 b_district[5] 2805.5472 0.9998584 b_district[6] 1618.8969 0.9998220 b_district[7] 1825.5960 1.0007015 b_district[8] 1679.4763 0.9997061 b_district[9] 2363.5877 0.9992498 b_district[10] 1563.5257 1.0045377 b_district[11] 1167.8123 1.0002489 b_district[12] 2498.4809 0.9991042 b_district[13] 1756.7462 1.0007052 b_district[14] 1228.6417 1.0001107 b_district[15] 1901.5656 0.9994154 b_district[16] 2039.6668 0.9990231 b_district[17] 2409.0061 0.9997818 b_district[18] 2201.4167 0.9996737 b_district[19] 1558.7646 1.0014130 b_district[20] 1943.8058 1.0000445 b_district[21] 822.0622 1.0015318 b_district[22] 1979.9077 0.9994562 b_district[23] 2290.9894 1.0005602 b_district[24] 1513.2038 0.9984700 b_district[25] 1832.4354 0.9993031 b_district[26] 2059.9716 0.9999822 b_district[27] 1594.2774 1.0006990 b_district[28] 2185.1060 0.9984959 b_district[29] 1720.6701 1.0000656 b_district[30] 1947.2835 0.9995516 ā® ā® ā® a_district[40] 1670.8647 0.9995937 a_district[41] 1792.6930 1.0000349 a_district[42] 1172.3509 0.9993053 a_district[43] 2325.3694 0.9993269 a_district[44] 3272.7154 0.9992261 a_district[45] 2406.2164 0.9993421 a_district[46] 2844.9845 0.9996790 a_district[47] 2090.3225 0.9991635 a_district[48] 2392.3548 1.0011711 a_district[49] 2326.2755 1.0012745 a_district[50] 2666.9974 1.0002669 a_district[51] 2380.2410 0.9985290 a_district[52] 1783.3829 0.9995301 a_district[53] 2138.6105 1.0030239 a_district[54] 1778.1201 0.9996700 a_district[55] 1840.4563 0.9994024 a_district[56] 2513.9127 0.9999035 a_district[57] 1381.4852 1.0009826 a_district[58] 1886.2134 1.0009924 a_district[59] 2061.4345 0.9985233 a_district[60] 1931.9196 1.0007637 a 1396.8717 1.0008768 b 1163.2596 1.0003428 bAge 3897.3416 0.9989784 sigma_intercepts_slopes[1] 645.1497 1.0061509 sigma_intercepts_slopes[2] 240.2612 1.0026197 Rho[1,1] NaN NaN Rho[1,2] 442.8782 1.0064511 Rho[2,1] 442.8782 1.0064511 Rho[2,2] 1692.5187 0.9979980  In the first model, we see that with more children the likelihood of using birth control goes up, as expected. The likelihood of using birth control also goes down with age, as expected if you assume that younger women are more open to newer ideas like birth control. In the second model we see the total causal effect of age on whether a woman uses birth control is positive. That is, the older a woman is the more likely she is to use birth control. That is, the tendency for a woman to use birth control as she gets older because she has more children is stronger than the influence of changing acceptance of birth control. 14H3. Modify any models from 14H2 that contained that children variable and model the variable now as a monotonic ordered category, like education from the week we did ordered categories. Education in that example had 8 categories. Children here will have fewer (no one in the sample had 8 children). So modify the code appropriately. What do you conclude about the causal influence of each additional child on use of contraception? ERROR: The author uses the term week above as if he has a syllabus for the book. Answer. Fitting the model: data(bangladesh) bc_df <- bangladesh bc_df$district_id <- as.integer(as.factor(bc_df$district)) bc_dat <- list( UseContraception = bc_df$use.contraception,
DistrictId = bc_df$district_id, Urban = bc_df$urban,
Age = standardize(bc_df$age.centered), Children = as.integer(bc_df$living.children),
alpha = rep(2, 3)
)
m_bc_ordered_children <- ulam(
alist(
UseContraception ~ dbinom(1, p),
logit(p) <- bC*sum(delta_j[1:Children]) + a_district[DistrictId] + b_district[DistrictId] *
Urban + bAge * Age,
c(a_district, b_district)[DistrictId] ~ multi_normal(c(a, b), Rho, sigma_intercepts_slopes),
bC ~ normal(0, 1),
a ~ normal(0, 2),
b ~ normal(0, 0.5),
bAge ~ normal(0, 0.5),
sigma_intercepts_slopes ~ exponential(1),
Rho ~ lkj_corr(2),
vector[4]: delta_j <<- append_row( 0 , delta ),
simplex[3]: delta ~ dirichlet( alpha )
),
data = bc_dat, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m_bc_ordered_children, depth=3), mimetypes="text/plain")
iplot(function() {
plot(precis(m_bc_ordered_children, depth=3), main="m_bc_ordered_children")
}, ar=0.4)

Warning message:
āBulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-essā

Warning message:
āTail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-essā

                           mean        sd           5.5%        94.5%
b_district[1]               1.1172469  0.3947964     0.48297366 1.7406465
b_district[2]               0.7614204  0.7128541    -0.32013969 1.8925728
b_district[3]               0.9775367  0.8165073    -0.28479200 2.3091791
b_district[4]               1.6304651  0.6217377     0.70311733 2.6565408
b_district[5]               0.6884590  0.6266713    -0.30505475 1.6964227
b_district[6]               1.3976158  0.5506682     0.58861327 2.3132589
b_district[7]               0.8752725  0.6839709    -0.19915593 1.9709745
b_district[8]               0.9775638  0.6223842     0.02468108 2.0054431
b_district[9]               1.1287229  0.6457095     0.15336211 2.2037115
b_district[10]              1.1993028  0.7611886     0.10176672 2.4951891
b_district[11]              1.5533671  0.8312594     0.32903001 2.9187281
b_district[12]              0.5132767  0.6062558    -0.46498870 1.4522622
b_district[13]              0.3750357  0.5881354    -0.59452498 1.2741826
b_district[14]              1.3850316  0.4405043     0.67993847 2.0776020
b_district[15]              0.4313439  0.6274121    -0.61220327 1.4028669
b_district[16]              0.5163426  0.6432999    -0.46993826 1.5422191
b_district[17]              0.8800332  0.7051507    -0.19434518 1.9770823
b_district[18]              1.0208456  0.4937756     0.23166464 1.8259519
b_district[19]              1.0040505  0.6176643     0.07376052 2.0267407
b_district[20]              0.4920475  0.7194264    -0.65689430 1.6608345
b_district[21]             -0.3435313  0.7086118    -1.55182460 0.6896701
b_district[22]              1.0893223  0.7290794    -0.05315765 2.2403613
b_district[23]              0.8678881  0.7129492    -0.25082946 2.0184497
b_district[24]              1.2810500  0.7522196     0.11836281 2.5333473
b_district[25]              0.4349194  0.4393901    -0.24534788 1.1284563
b_district[26]              0.6664842  0.7000972    -0.46075819 1.7295660
b_district[27]              1.1723796  0.6142847     0.21989934 2.1488261
b_district[28]              0.6975674  0.6067147    -0.28778595 1.6661818
b_district[29]              1.1837080  0.6063302     0.27756238 2.1866475
b_district[30]              0.9411080  0.4801367     0.20010201 1.7490492
ā®                          ā®           ā®            ā®           ā®
a_district[44]             -2.01565316 3.889580e-01 -2.65596603 -1.4197053
a_district[45]             -1.97231042 3.553071e-01 -2.54360962 -1.4135556
a_district[46]             -1.05993405 2.479984e-01 -1.47103244 -0.6728968
a_district[47]             -1.40711734 4.719714e-01 -2.15791705 -0.6477676
a_district[48]             -1.09022149 3.433844e-01 -1.62885343 -0.5490830
a_district[49]             -1.97443974 5.736895e-01 -2.92491602 -1.0750045
a_district[50]             -1.52572109 4.364849e-01 -2.21756447 -0.8397019
a_district[51]             -1.69507483 3.930730e-01 -2.33437232 -1.0770527
a_district[52]             -1.09322637 3.058946e-01 -1.58148714 -0.5976245
a_district[53]             -1.69162223 6.504140e-01 -2.71149771 -0.6315280
a_district[54]             -1.72704061 6.375521e-01 -2.75649355 -0.7129678
a_district[55]             -0.96697391 3.694604e-01 -1.54392148 -0.3801536
a_district[56]             -2.13901583 4.255826e-01 -2.83738719 -1.4699294
a_district[57]             -1.00697529 3.769271e-01 -1.59853190 -0.4072379
a_district[58]             -2.24133056 5.284252e-01 -3.12622964 -1.4170601
a_district[59]             -2.12211217 4.168315e-01 -2.81841840 -1.4745107
a_district[60]             -2.23019533 3.914657e-01 -2.89008046 -1.6357820
bC                          1.38506932 1.591415e-01  1.13804534  1.6396411
a                          -1.65171232 1.529888e-01 -1.90285604 -1.4009028
b                           0.70576856 1.594368e-01  0.45896996  0.9545962
bAge                       -0.25808498 6.446993e-02 -0.36040392 -0.1548705
sigma_intercepts_slopes[1]  0.60851857 1.033183e-01  0.45612360  0.7832236
sigma_intercepts_slopes[2]  0.77719933 1.977456e-01  0.48570605  1.1170490
Rho[1,1]                    1.00000000 0.000000e+00  1.00000000  1.0000000
Rho[1,2]                   -0.63287420 1.694468e-01 -0.85518047 -0.3378406
Rho[2,1]                   -0.63287420 1.694468e-01 -0.85518047 -0.3378406
Rho[2,2]                    1.00000000 6.513255e-17  1.00000000  1.0000000
delta[1]                    0.73575546 7.864344e-02  0.60302275  0.8571040
delta[2]                    0.16514339 7.786641e-02  0.05346976  0.3011820
delta[3]                    0.09910115 5.303689e-02  0.02646703  0.1976305
n_eff     Rhat4
b_district[1]              1858.6535 1.0021939
b_district[2]              2422.5678 0.9998569
b_district[3]              2364.3998 1.0004502
b_district[4]               799.8742 1.0005938
b_district[5]              2641.5894 0.9984578
b_district[6]              1016.8648 1.0012057
b_district[7]              2286.4364 0.9990444
b_district[8]              1947.3531 1.0010813
b_district[9]              1700.2795 1.0002677
b_district[10]             1220.2683 1.0032249
b_district[11]              908.4975 1.0035830
b_district[12]             1876.7732 0.9991535
b_district[13]             2190.9224 1.0006425
b_district[14]             1216.8354 1.0014568
b_district[15]             1807.2520 0.9985229
b_district[16]             1998.1333 1.0005741
b_district[17]             2069.1264 1.0008744
b_district[18]             2138.1279 0.9998421
b_district[19]             2460.2234 0.9988916
b_district[20]             2318.7812 0.9983729
b_district[21]              726.9583 1.0003159
b_district[22]             1815.5166 0.9988166
b_district[23]             2638.9181 1.0000360
b_district[24]             1296.3829 1.0019795
b_district[25]             2110.8114 1.0005550
b_district[26]             2391.6174 0.9994238
b_district[27]             1657.3948 1.0023487
b_district[28]             1419.3083 1.0008453
b_district[29]             1603.2549 0.9994156
b_district[30]             2471.6563 0.9992102
ā®                          ā®         ā®
a_district[44]             1409.2554 1.0008529
a_district[45]             1064.7414 1.0004311
a_district[46]             1543.1175 0.9996723
a_district[47]             1993.7202 1.0010994
a_district[48]             1475.3781 0.9997769
a_district[49]             1532.6337 0.9992893
a_district[50]             1927.3143 1.0005260
a_district[51]             1398.2887 1.0005994
a_district[52]             1627.6985 0.9995263
a_district[53]             1621.6553 1.0004567
a_district[54]             1710.0057 0.9999061
a_district[55]             2001.8494 1.0008422
a_district[56]             1499.8706 1.0055149
a_district[57]             1729.1973 1.0018148
a_district[58]             1560.1178 1.0038867
a_district[59]             1365.4777 1.0006630
a_district[60]             1717.9668 1.0004389
bC                          555.2256 1.0100292
a                           566.5249 1.0060587
b                           965.3262 1.0032480
bAge                       1144.8212 1.0040346
sigma_intercepts_slopes[1]  587.3616 1.0079381
sigma_intercepts_slopes[2]  248.2555 1.0146701
Rho[1,1]                         NaN       NaN
Rho[1,2]                    463.0337 1.0047042
Rho[2,1]                    463.0337 1.0047042
Rho[2,2]                   2144.0760 0.9979980
delta[1]                   3036.3896 1.0011396
delta[2]                   2264.5731 1.0002415
delta[3]                   2666.5015 0.9993974


Notice in the summary of this data frame, above, that the women surveyed had a minimum of one and a maximum of four children. It appears the maximum jump in likelihood of a woman using contraception occurs after having the second child.

14H4. Varying effects models are useful for modeling time series, as well as spatial clustering. In a time series, the observations cluster by entities that have continuity through time, such as individuals. Since observations within individuals are likely highly correlated, the multilevel structure can help quite a lot. Youāll use the data in data(Oxboys), which is 234 height measurements on 26 boys from an Oxford Boys Club (I think these were like youth athletic leagues?), at 9 different ages (centered and standardized) per boy. Youāll be interested in predicting height, using age, clustered by Subject (individual boy). Fit a model with varying intercepts and slopes (on age), clustered by Subject. Present and interpret the parameter estimates. Which varying effect contributes more variation to the heights, the intercept or the slope?

data(Oxboys)
ox <- Oxboys


Answer. A head and summary of the Oxboys data.frame (no help exists):

display(head(ox))
display(summary(ox))

A data.frame: 6 Ć 4
SubjectageheightOccasion
<int><dbl><dbl><int>
11-1.0000140.51
21-0.7479143.42
31-0.4630144.83
41-0.1643147.14
51-0.0027147.75
61 0.2466150.26
    Subject          age               height         Occasion
Min.   : 1.0   Min.   :-1.00000   Min.   :126.2   Min.   :1
1st Qu.: 7.0   1st Qu.:-0.46300   1st Qu.:143.8   1st Qu.:3
Median :13.5   Median :-0.00270   Median :149.5   Median :5
Mean   :13.5   Mean   : 0.02263   Mean   :149.5   Mean   :5
3rd Qu.:20.0   3rd Qu.: 0.55620   3rd Qu.:155.5   3rd Qu.:7
Max.   :26.0   Max.   : 1.00550   Max.   :174.8   Max.   :9


Fitting the model:

ox_dat <- list(
Subject = ox$Subject, Age = ox$age,
Height = ox$height ) m_ox <- ulam( alist( Height ~ dnorm(mu, sigma), mu <- a_boy[Subject] + b_boy[Subject] * Age, c(a_boy, b_boy)[Subject] ~ multi_normal(c(a, b), Rho, sigma_intercepts_slopes), a ~ normal(0, 1), b ~ normal(0, 0.5), sigma_intercepts_slopes ~ exponential(1), Rho ~ lkj_corr(2), sigma ~ normal(0, 2) ), data = ox_dat, chains = 4, cores = 4, log_lik = TRUE ) display(precis(m_ox, depth=3), mimetypes="text/plain") iplot(function() { plot(precis(m_ox, depth=3), main="m_ox") }, ar=0.7)   mean sd 5.5% 94.5% b_boy[1] 7.137421 0.3470028 6.593290 7.685841 b_boy[2] 5.481388 0.3354717 4.953505 6.018759 b_boy[3] 4.871198 0.3242252 4.349729 5.384686 b_boy[4] 9.248308 0.3286225 8.723941 9.775897 b_boy[5] 6.262325 0.3430345 5.720164 6.794605 b_boy[6] 4.065094 0.3414588 3.507843 4.625468 b_boy[7] 5.057129 0.3254772 4.550538 5.600158 b_boy[8] 6.472143 0.3409617 5.918703 7.011299 b_boy[9] 6.008275 0.3482750 5.451821 6.546127 b_boy[10] 3.822563 0.3308904 3.288179 4.352849 b_boy[11] 8.414984 0.3420993 7.896069 8.966929 b_boy[12] 7.003539 0.3169374 6.496528 7.504444 b_boy[13] 8.411397 0.3244596 7.903736 8.921488 b_boy[14] 8.584270 0.3479208 8.024445 9.140723 b_boy[15] 7.084154 0.3311591 6.576340 7.612390 b_boy[16] 4.632670 0.3351358 4.111953 5.170819 b_boy[17] 8.503973 0.3483309 7.942851 9.053097 b_boy[18] 5.993513 0.3326001 5.483076 6.526629 b_boy[19] 8.976750 0.3399502 8.432994 9.514783 b_boy[20] 4.472405 0.3421088 3.928483 5.025561 b_boy[21] 7.459345 0.3350197 6.935254 7.996470 b_boy[22] 8.026042 0.3425854 7.475904 8.563561 b_boy[23] 7.155082 0.3293630 6.633925 7.688765 b_boy[24] 6.758130 0.3489079 6.217386 7.324170 b_boy[25] 4.119345 0.3261452 3.618808 4.628240 b_boy[26] 5.611484 0.3198873 5.097121 6.124518 a_boy[1] 148.113666 0.2214899 147.760040 148.470348 a_boy[2] 142.852459 0.2277650 142.483513 143.220769 a_boy[3] 155.647841 0.2273460 155.290838 156.013241 a_boy[4] 165.077054 0.2245467 164.719447 165.438460 ā® ā® ā® ā® ā® a_boy[6] 146.7816740 0.21651718 146.4285482 147.1306175 a_boy[7] 146.1256625 0.22946990 145.7590553 146.4807607 a_boy[8] 148.2952714 0.22680539 147.9407225 148.6682984 a_boy[9] 138.1430095 0.21399152 137.8103545 138.4762010 a_boy[10] 130.2641576 0.22281089 129.9045046 130.6156562 a_boy[11] 150.0505071 0.22596823 149.6941118 150.4049100 a_boy[12] 156.8038612 0.21954019 156.4661561 157.1615807 a_boy[13] 156.0790796 0.22703441 155.7234655 156.4528918 a_boy[14] 159.4761872 0.22282774 159.1283867 159.8305711 a_boy[15] 144.2824583 0.22543143 143.9169227 144.6406060 a_boy[16] 147.5473591 0.22410248 147.1885303 147.8914143 a_boy[17] 142.9797313 0.22542690 142.6231427 143.3483465 a_boy[18] 151.1778984 0.21780214 150.8322288 151.5216472 a_boy[19] 164.5871053 0.23456883 164.2127483 164.9628855 a_boy[20] 151.4691104 0.21856394 151.1150465 151.8145878 a_boy[21] 150.5169563 0.21928624 150.1524631 150.8591893 a_boy[22] 154.5666657 0.22340787 154.2099223 154.9231753 a_boy[23] 151.0604161 0.23528629 150.6893835 151.4446547 a_boy[24] 153.1466540 0.23028317 152.7808444 153.4993456 a_boy[25] 139.2048731 0.21329264 138.8642392 139.5451884 a_boy[26] 137.9953264 0.21980365 137.6500644 138.3334579 a 0.6537408 1.01336819 -0.9283681 2.2715824 b 0.1333895 0.48955902 -0.6369618 0.9309867 sigma_intercepts_slopes[1] 74.6934922 4.68032793 67.4601565 82.1614231 sigma_intercepts_slopes[2] 3.5748334 0.32779123 3.0554130 4.1134228 Rho[1,1] 1.0000000 0.00000000 1.0000000 1.0000000 Rho[1,2] 0.8880278 0.03607893 0.8246749 0.9352457 Rho[2,1] 0.8880278 0.03607893 0.8246749 0.9352457 Rho[2,2] 1.0000000 0.00000000 1.0000000 1.0000000 sigma 0.6647896 0.03639783 0.6084315 0.7252927 n_eff Rhat4 b_boy[1] 3528.749 0.9999860 b_boy[2] 4168.444 0.9996024 b_boy[3] 2986.183 0.9986559 b_boy[4] 3184.555 0.9988569 b_boy[5] 4628.295 0.9986748 b_boy[6] 3157.599 0.9986585 b_boy[7] 3067.210 1.0018310 b_boy[8] 3605.736 0.9991257 b_boy[9] 3235.929 0.9989471 b_boy[10] 3524.432 0.9991593 b_boy[11] 3555.813 0.9985748 b_boy[12] 3295.006 0.9983287 b_boy[13] 3261.745 0.9985648 b_boy[14] 3638.736 1.0001942 b_boy[15] 3922.221 0.9984418 b_boy[16] 3387.989 0.9991188 b_boy[17] 3695.833 0.9999522 b_boy[18] 4232.337 0.9991695 b_boy[19] 3671.879 0.9997550 b_boy[20] 2721.981 0.9992497 b_boy[21] 3176.938 0.9983887 b_boy[22] 3561.369 0.9994862 b_boy[23] 3902.230 0.9987558 b_boy[24] 3227.228 1.0002230 b_boy[25] 3818.300 0.9996149 b_boy[26] 3410.292 1.0004828 a_boy[1] 4561.809 0.9987926 a_boy[2] 3535.504 0.9998455 a_boy[3] 3005.947 0.9996865 a_boy[4] 3959.527 0.9986625 ā® ā® ā® a_boy[6] 3586.519 0.9993449 a_boy[7] 3621.095 0.9999323 a_boy[8] 2529.734 0.9983566 a_boy[9] 2600.322 0.9996230 a_boy[10] 3124.062 0.9989691 a_boy[11] 3161.117 0.9988722 a_boy[12] 3231.492 0.9985778 a_boy[13] 3953.531 0.9989110 a_boy[14] 2868.897 0.9992057 a_boy[15] 3993.774 1.0000635 a_boy[16] 3126.013 0.9988499 a_boy[17] 3317.890 0.9990095 a_boy[18] 3539.491 0.9989600 a_boy[19] 2805.662 0.9986506 a_boy[20] 3454.852 0.9985530 a_boy[21] 3114.616 0.9984869 a_boy[22] 3719.990 1.0003180 a_boy[23] 3721.650 0.9988570 a_boy[24] 3812.427 0.9991344 a_boy[25] 3265.778 0.9990980 a_boy[26] 3590.703 0.9994189 a 4519.370 0.9981847 b 2040.379 0.9995646 sigma_intercepts_slopes[1] 2474.403 0.9992004 sigma_intercepts_slopes[2] 2001.709 0.9986325 Rho[1,1] NaN NaN Rho[1,2] 2216.757 0.9990147 Rho[2,1] 2216.757 0.9990147 Rho[2,2] NaN NaN sigma 2050.545 1.0026862  The intercepts term contributes much more to the heights than the slopes. Notice the b_boy terms are the maximum that the slopes can contribute to mu since Age is at most equal one. 14H5. Now consider the correlation between the varying intercepts and slopes. Can you explain its value? How would this estimated correlation influence your predictions about a new sample of boys? Answer. The correlation is positive, because more positive intercepts are associated with more positive slopes. That is, boys who are on average taller also grow faster. 14H6. Use mvrnorm (in library(MASS)) or rmvnorm (in library(mvtnorm)) to simulate a new sample of boys, based upon the posterior mean values of the parameters. That is, try to simulate varying intercepts and slopes, using the relevant parameter estimates, and then plot the predicted trends of height on age, one trend for each simulated boy you produce. A sample of 10 simulated boys is plenty, to illustrate the lesson. You can ignore uncertainty in the posterior, just to make the problem a little easier. But if you want to include the uncertainty about the parameters, go for it. Note that you can construct an arbitrary variance-covariance matrix to pass to either mvrnorm or rmvnorm with something like: R code 14.54 S <- matrix( c( sa^2 , sa*sb*rho , sa*sb*rho , sb^2 ) , nrow=2 )  where sa is the standard deviation of the first variable, sb is the standard deviation of the second variable, and rho is the correlation between them. Answer. The expected plot: post <- extract.samples(m_ox) Mu_est <- c(mean(post$a), mean(post$b)) rho_est <- mean(post$Rho[,1,2])
sa_est <- mean(post$sigma_intercepts_slopes[,1]) sb_est <- mean(post$sigma_intercepts_slopes[,2])
cov_ab <- sa_est*sb_est*rho_est
Sigma_est <- matrix(c(sa_est^2, cov_ab, cov_ab, sb_est^2), nrow=2)

library(MASS)
set.seed(5)
N_boys <- 10
vary_effects <- mvrnorm(N_boys, Mu_est, Sigma_est)

iplot(function() {
plot(NULL, main="Posterior predictive plot, 10 simulated boys",
xlab="age", ylab="height",
xlim=c(-1.1, 1.1),
ylim=c(min(vary_effects[,1] - vary_effects[,2]) - 10, max(vary_effects[,1] + vary_effects[,2]) + 10)
)
for (idx in 1:N_boys)
abline(a=vary_effects[idx, 1], b=vary_effects[idx, 2])
})


The lesson here is that vague priors can kill predictive power, and that in situations like this one the priors have a strong influence on the results. In the previous plot the priors on a and b were copied and pasted from another model:

\begin{split} \begin{align} a & \sim Normal(0, 1) \\ b & \sim Normal(0, 0.5) \end{align} \end{split}

These are clearly unreasonable. If we fix these priors we get a much better result:

\begin{split} \begin{align} a & \sim Normal(150, 40) \\ b & \sim Normal(40, 40) \end{align} \end{split}
m_ox2 <- ulam(
alist(
Height ~ dnorm(mu, sigma),
mu <- a_boy[Subject] + b_boy[Subject] * Age,
c(a_boy, b_boy)[Subject] ~ multi_normal(c(a, b), Rho, sigma_intercepts_slopes),
a ~ normal(150, 40),
b ~ normal(40, 40),
sigma_intercepts_slopes ~ exponential(1),
Rho ~ lkj_corr(2),
sigma ~ normal(0, 2)
),
data = ox_dat, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m_ox2, depth=3), mimetypes="text/plain")
iplot(function() {
plot(precis(m_ox2, depth=3), main="m_ox2")
}, ar=0.7)

post <- extract.samples(m_ox2)
Mu_est <- c(mean(post$a), mean(post$b))
rho_est <- mean(post$Rho[,1,2]) sa_est <- mean(post$sigma_intercepts_slopes[,1])
sb_est <- mean(post\$sigma_intercepts_slopes[,2])
cov_ab <- sa_est*sb_est*rho_est
Sigma_est <- matrix(c(sa_est^2, cov_ab, cov_ab, sb_est^2), nrow=2)

library(MASS)
set.seed(5)
N_boys <- 10
vary_effects <- mvrnorm(N_boys, Mu_est, Sigma_est)

iplot(function() {
plot(NULL, main="Posterior predictive plot, 10 simulated boys",
xlab="age", ylab="height",
xlim=c(-1.1, 1.1),
ylim=c(min(vary_effects[,1] - vary_effects[,2]) - 10, max(vary_effects[,1] + vary_effects[,2]) + 10)
)
for (idx in 1:N_boys)
abline(a=vary_effects[idx, 1], b=vary_effects[idx, 2])
})

                           mean        sd           5.5%        94.5%
b_boy[1]                     7.122827  0.3322087      6.583848    7.655954
b_boy[2]                     5.460675  0.3309538      4.937802    6.000382
b_boy[3]                     4.934086  0.3372113      4.401751    5.461362
b_boy[4]                     9.290755  0.3408733      8.753078    9.830426
b_boy[5]                     6.279010  0.3414299      5.738504    6.831309
b_boy[6]                     4.107758  0.3332389      3.586919    4.640147
b_boy[7]                     5.064438  0.3433725      4.520589    5.604039
b_boy[8]                     6.463974  0.3402083      5.927326    6.991294
b_boy[9]                     5.958746  0.3234855      5.461362    6.483401
b_boy[10]                    3.754424  0.3422636      3.209548    4.298276
b_boy[11]                    8.376141  0.3459315      7.827744    8.901756
b_boy[12]                    7.037235  0.3251566      6.529788    7.552236
b_boy[13]                    8.421781  0.3320623      7.887876    8.969216
b_boy[14]                    8.606473  0.3324009      8.088091    9.125618
b_boy[15]                    7.052871  0.3329311      6.496938    7.582998
b_boy[16]                    4.648309  0.3358858      4.131960    5.187072
b_boy[17]                    8.427887  0.3242495      7.909494    8.945488
b_boy[18]                    6.003476  0.3500593      5.438516    6.563315
b_boy[19]                    9.016830  0.3303680      8.477545    9.540540
b_boy[20]                    4.527208  0.3365359      3.984102    5.066829
b_boy[21]                    7.447827  0.3292831      6.908307    7.986375
b_boy[22]                    8.023659  0.3309146      7.499065    8.543218
b_boy[23]                    7.152485  0.3387705      6.633455    7.682342
b_boy[24]                    6.784337  0.3493435      6.229581    7.355831
b_boy[25]                    4.106523  0.3317459      3.575951    4.634723
b_boy[26]                    5.557906  0.3539407      4.985938    6.127773
a_boy[1]                   148.125102  0.2240055    147.756752  148.485698
a_boy[2]                   142.868576  0.2203292    142.518234  143.210920
a_boy[3]                   155.632299  0.2254481    155.266723  155.978131
a_boy[4]                   165.064997  0.2285159    164.709111  165.443212
ā®                          ā®           ā®            ā®           ā®
a_boy[6]                   146.7861098 2.194073e-01 146.4384452 147.1318658
a_boy[7]                   146.1210404 2.337751e-01 145.7343224 146.4824160
a_boy[8]                   148.2930374 2.243746e-01 147.9392885 148.6591341
a_boy[9]                   138.1513979 2.313018e-01 137.7888748 138.5179245
a_boy[10]                  130.2768959 2.161344e-01 129.9233319 130.6063412
a_boy[11]                  150.0537839 2.192917e-01 149.7113235 150.4125345
a_boy[12]                  156.7994271 2.247133e-01 156.4529092 157.1583169
a_boy[13]                  156.0717378 2.231343e-01 155.6971782 156.4154992
a_boy[14]                  159.4622504 2.268403e-01 159.1051252 159.8277467
a_boy[15]                  144.2840014 2.288367e-01 143.9157850 144.6475227
a_boy[16]                  147.5418819 2.208381e-01 147.1855594 147.8994849
a_boy[17]                  142.9979601 2.268262e-01 142.6356783 143.3596771
a_boy[18]                  151.1794610 2.169724e-01 150.8242030 151.5177369
a_boy[19]                  164.5677562 2.193259e-01 164.2165483 164.9211439
a_boy[20]                  151.4573911 2.289793e-01 151.1095460 151.8243764
a_boy[21]                  150.5201332 2.205926e-01 150.1666056 150.8674592
a_boy[22]                  154.5665399 2.197841e-01 154.2112713 154.9120035
a_boy[23]                  151.0684985 2.245058e-01 150.7159832 151.4270030
a_boy[24]                  153.1327253 2.242584e-01 152.7759679 153.4905381
a_boy[25]                  139.2209405 2.212194e-01 138.8686246 139.5646712
a_boy[26]                  138.0036029 2.302090e-01 137.6423519 138.3753935
a                          149.3923945 1.455555e+00 147.0672831 151.7613879
b                            6.5311084 3.284126e-01   6.0277505   7.0514299
sigma_intercepts_slopes[1]   7.3280618 8.909815e-01   6.0736747   8.8576264
sigma_intercepts_slopes[2]   1.6304761 2.332820e-01   1.2976833   2.0271578
Rho[1,1]                     1.0000000 0.000000e+00   1.0000000   1.0000000
Rho[1,2]                     0.5257282 1.299665e-01   0.2859021   0.7119550
Rho[2,1]                     0.5257282 1.299665e-01   0.2859021   0.7119550
Rho[2,2]                     1.0000000 7.984801e-17   1.0000000   1.0000000
sigma                        0.6652505 3.387909e-02   0.6144083   0.7236621
n_eff    Rhat4
b_boy[1]                   3397.241 0.9993882
b_boy[2]                   3824.896 0.9991548
b_boy[3]                   4551.947 0.9986551
b_boy[4]                   4846.730 0.9984812
b_boy[5]                   3796.926 0.9981290
b_boy[6]                   3969.628 0.9994277
b_boy[7]                   3717.169 0.9995144
b_boy[8]                   3885.725 0.9988696
b_boy[9]                   3130.566 0.9982011
b_boy[10]                  3829.525 0.9996160
b_boy[11]                  4275.056 0.9990806
b_boy[12]                  3486.493 0.9986922
b_boy[13]                  3565.861 0.9988748
b_boy[14]                  5191.336 0.9989458
b_boy[15]                  3540.452 0.9982362
b_boy[16]                  3274.932 0.9995209
b_boy[17]                  3718.099 0.9982985
b_boy[18]                  3729.188 0.9992369
b_boy[19]                  4254.105 0.9986004
b_boy[20]                  3826.369 0.9985832
b_boy[21]                  3864.268 0.9989389
b_boy[22]                  3111.481 0.9993075
b_boy[23]                  3595.851 0.9989247
b_boy[24]                  3738.130 0.9983798
b_boy[25]                  3236.979 0.9988405
b_boy[26]                  4879.405 0.9984281
a_boy[1]                   3237.550 0.9987787
a_boy[2]                   4049.226 0.9986528
a_boy[3]                   3117.776 0.9985072
a_boy[4]                   5100.197 0.9987063
ā®                          ā®        ā®
a_boy[6]                   3308.713 1.0001489
a_boy[7]                   3880.149 0.9990667
a_boy[8]                   3317.791 0.9991302
a_boy[9]                   5297.428 0.9994021
a_boy[10]                  4105.215 1.0000839
a_boy[11]                  3896.860 0.9988511
a_boy[12]                  3777.433 0.9989543
a_boy[13]                  4259.095 0.9990450
a_boy[14]                  4502.372 0.9996217
a_boy[15]                  4221.052 0.9987526
a_boy[16]                  4264.674 0.9988805
a_boy[17]                  3706.006 0.9986395
a_boy[18]                  3493.254 0.9986243
a_boy[19]                  3833.874 0.9992270
a_boy[20]                  3961.773 0.9986611
a_boy[21]                  4278.342 0.9984816
a_boy[22]                  4313.655 0.9987957
a_boy[23]                  5199.071 0.9988071
a_boy[24]                  4171.203 0.9989638
a_boy[25]                  4217.568 0.9995400
a_boy[26]                  4517.137 0.9988017
a                          2719.812 0.9993785
b                          2536.216 1.0007571
sigma_intercepts_slopes[1] 2971.258 0.9986441
sigma_intercepts_slopes[2] 2646.944 0.9995006
Rho[1,1]                        NaN       NaN
Rho[1,2]                   2607.412 0.9985915
Rho[2,1]                   2607.412 0.9985915
Rho[2,2]                   1956.060 0.9979980
sigma                      2380.168 0.9985868