This chapter extended the basic multilevel strategy of partial pooling to slopes as well as intercepts. Accomplishing this meant modeling covariation in the statistical population of parameters. The LKJcorr prior was introduced as a convenient family of priors for correlation matrices. You saw how covariance models can be applied to causal inference, using instrumental variables and the front-door criterion. Gaussian processes represent a practical method of extending the varying effects strategy to continuous dimensions of similarity, such as spatial, network, phylogenetic, or any other abstract distance between entities in the data.
Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them. Problems are labeled Easy (E), Medium (M), and Hard(H).
Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points.
14E1. Add to the following model varying slopes on the predictor x. \[\begin{align} \ y_i ∼ Normal(μ_i, σ) \\ \ μ_i = α_{group[i]} + βx_i \\ \ α_{group} ∼ Normal(α, σ_α) \\ \ α ∼ Normal(0, 10) \\ \ β ∼ Normal(0, 1) \\ \ σ ∼ Exponential(1) \\ \ σ_α ∼ Exponential(1) \\ \end{align}\]
# μ_i = α_{group[i]} + β_{group[i]}x_i14E2. Think up a context in which varying intercepts will be positively correlated with varying slopes. Provide a mechanistic explanation for the correlation.
# Let’s say we are interested studying ant hill size as a function of food availability. We can reasonably expect larger ant hills (higher intercepts) to benefit more strongly from increased food availability as their foraging will be much more efficient (steeper slope).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.
# When there is little or next-to-no variation among clusters. The absence of this among-cluster variation induces very strong shrinkage.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. Plot the posterior distribution. How does the posterior distribution of the correlation reflect this change in the underlying simulation?
a <- 3.5
b <- (-1)
sigma_a <- 1
sigma_b <- 0.5
rho <- 0
Mu <- c(a, b)
cov_ab <- sigma_a * sigma_b * rho
Sigma <- matrix(c(sigma_a^2, cov_ab, cov_ab, sigma_b^2), ncol = 2)
N_cafes <- 20
set.seed(6)
vary_effects <- mvrnorm(N_cafes, Mu, Sigma)
a_cafe <- vary_effects[, 1]
b_cafe <- vary_effects[, 2]
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
wait <- rnorm(N_visits * N_cafes, mu, sigma)
d <- data.frame(cafe = cafe_id, afternoon = afternoon, wait = wait)
#m_M1 <- 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 = 6, cores = 6
#)
#post <- extract.samples(m_M1)
#ggplot() +
# stat_halfeye(aes(x = post$Rho[, 1, 2])) +
#theme_bw() +
#labs(x = "Rho")
#precis(m_M1, pars = "Rho[1,2]")Due error I had to comment out the code, the model code was same as in the book
14M2. Fit this multilevel model to the simulated café data: \[\begin{align} \ W_i ∼ Normal(μ_i, σ) \\ \ μ_i = α_{CAFE[i]} + β_{CAFE[i]}A_i \\ \ α_{CAFE} ∼ Normal(α, σ_α) \\ \ β_{CAFE} ∼ Normal(β, σ_β) \\ \ α ∼ Normal(0, 10) \\ \ β ∼ Normal(0, 10) \\ \ σ, σ_α, σ_β ∼ Exponential(1) \\ \end{align}\] Use WAIC to compare this model to the model from the chapter, the one that uses a multi-variate Gaussian prior. Visualize the differences between the posterior means. Explain the result.
a <- 3.5
b <- -1
sigma_a <- 1
sigma_b <- 0.5
rho <- -0.7
Mu <- c(a, b)
cov_ab <- sigma_a * sigma_b * rho
Sigma <- matrix(c(sigma_a^2, cov_ab, cov_ab, sigma_b^2), ncol = 2)
N_cafes <- 20
set.seed(42)
vary_effects <- mvrnorm(N_cafes, Mu, Sigma)
a_cafe <- vary_effects[, 1]
b_cafe <- vary_effects[, 2]
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
wait <- rnorm(N_visits * N_cafes, mu, sigma)
d <- data.frame(cafe = cafe_id, afternoon = afternoon, wait = wait)
#m_M2Baseline <- 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
#)
#m_M2 <- ulam(
# alist(
# wait ~ dnorm(mu, sigma),
# mu <- a_cafe[cafe] + b_cafe[cafe] * afternoon,
#a_cafe[cafe] ~ dnorm(a, sigma_alpha),
#b_cafe[cafe] ~ dnorm(b, sigma_beta),
#a ~ dnorm(0, 10),
#b ~ dnorm(0, 10),
#sigma ~ dexp(1),
#sigma_alpha ~ dexp(1),
#sigma_beta ~ dexp(1)
#),
#data = d, chains = 4, cores = 4
#)
#post_Base <- extract.samples(m_M2Baseline)
#a_Base <- apply(post_Base$a_cafe, 2, mean)
#b_Base <- apply(post_Base$b_cafe, 2, mean)
#post_M2 <- extract.samples(m_M2)
#a_M2 <- apply(post_M2$a_cafe, 2, mean)
#b_M2 <- apply(post_M2$b_cafe, 2, mean)
#plot(a_M2, b_M2,
# xlab = "intercept", ylab = "slope",
#pch = 16, col = rangi2, ylim = c(min(b_M2) - 0.05, max(b_M2) + 0.05),
#xlim = c(min(a_M2) - 0.1, max(a_M2) + 0.1), cex = 2
#)
#points(a_Base, b_Base, pch = 1, cex = 2)Due error I had to comment out the code
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?
data(UCBadmit)
d <- UCBadmit
dat_list <- list(
admit = d$admit,
applications = d$applications,
male = ifelse(d$applicant.gender == "male", 1, 0),
dept_id = rep(1:6, each = 2)
)
str(dat_list)## List of 4
## $ admit : int [1:12] 512 89 353 17 120 202 138 131 53 94 ...
## $ applications: int [1:12] 825 108 560 25 325 593 417 375 191 393 ...
## $ male : num [1:12] 1 0 1 0 1 0 1 0 1 0 ...
## $ dept_id : int [1:12] 1 1 2 2 3 3 4 4 5 5 ...#Begin_C <- Sys.time()
#m_M3 <- ulam(
# alist(
# admit ~ dbinom(applications, p),
# logit(p) <- a[dept_id] + bm[dept_id] * male,
#c(a, bm)[dept_id] ~ multi_normal(c(a_bar, bm_bar), Rho, sigma_dept),
#a_bar ~ dnorm(0, 1),
#bm_bar ~ dnorm(0, 1),
#sigma_dept ~ dexp(1),
#Rho ~ dlkjcorr(2)
#),
#data = dat_list, chains = 4, cores = 4
#)
#End_C <- Sys.time()
#precis(m_M3, 3)
#Begin_NC <- Sys.time()
#m_M3NonCent <- ulam(
# alist(
# admit ~ dbinom(applications, p),
# logit(p) <- a_bar + v[dept_id, 1] + (bm_bar + v[dept_id, 2]) * male,
#transpars > matrix[dept_id, 2]:v <- compose_noncentered(sigma_dept, L_Rho, z),
#matrix[2, dept_id]:z ~ dnorm(0, 1),
#a_bar ~ dnorm(0, 1.5),
#bm_bar ~ dnorm(0, 1),
#vector[2]:sigma_dept ~ dexp(1),
#cholesky_factor_corr[2]:L_Rho ~ lkj_corr_cholesky(2)
#),
#data = dat_list, chains = 4, cores = 4
#)
#End_NC <- Sys.time()
#precis(m_M3NonCent, 3)
#End_C - Begin_C
#End_NC - Begin_NCDue error I had to comment out the code
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.
data(Kline)
d <- Kline
d$P <- scale(log(d$population))
d$contact_id <- ifelse(d$contact == "high", 2, 1)
d$society <- 1:10
data(islandsDistMatrix)
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, cores = 4, log_lik = TRUE
)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Warning in .local(object, ...): some chains had errors; consider specifying
## chains = 1 to debug## [[1]]
## Stan model '58422f20040c774e9740e486280fe76b' does not contain samples.
##
## [[2]]
## Stan model '58422f20040c774e9740e486280fe76b' does not contain samples.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
)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/3.6/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Warning in .local(object, ...): some chains had errors; consider specifying
## chains = 1 to debug## [[1]]
## Stan model '88f2ad7a384c819633238280368dab9f' does not contain samples.
##
## [[2]]
## Stan model '88f2ad7a384c819633238280368dab9f' does not contain samples.compare(m11.11, m14.8)## WAIC SE dWAIC dSE pWAIC weight
## m14.8 67.55575 2.327581 0.00000 NA 4.045784 0.998485131
## m11.11 80.53757 11.065026 12.98182 10.89504 5.096207 0.00151486914M5. Modify the phylogenetic distance example to use group size as the outcome and brain size as a predictor. Plot the estimates. Assuming brain size influences group size, what is your estimate of the effect? How does phylogeny influence the estimate?
data(Primates301)
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))
)
data(Primates301_nex)
tree_trimmed <- keep.tip(Primates301_nex, spp_obs)
Rbm <- corBrownian(phy = tree_trimmed)
V <- vcv(Rbm) ## Warning 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.Dmat <- cophenetic(tree_trimmed)
dat_list$V <- V[spp_obs, spp_obs]
dat_list$R <- dat_list$V / max(V)
#m_M5Ordi <- 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
#)
#m_M5Brown <- 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
#)
#dat_list$Dmat <- Dmat[spp_obs, spp_obs] / max(Dmat)
#m_M5GP <- 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
#)
#plot(coeftab(m_M5Ordi, m_M5Brown, m_M5GP), pars = c("bM", "bB"))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.
data(bangladesh)
d <- bangladesh
dat_list <- list(
C = d$use.contraception,
did = as.integer(as.factor(d$district)),
urban = d$urban
)
#m_H1 <- ulam(
# alist(
# C ~ bernoulli(p),
# logit(p) <- a[did] + b[did] * urban,
#c(a, b)[did] ~ multi_normal(c(abar, bbar), Rho, Sigma),
#abar ~ normal(0, 1),
#bbar ~ normal(0, 0.5),
#Rho ~ lkj_corr(2),
#Sigma ~ exponential(1)
#),
#data = dat_list, chains = 4, cores = 4, iter = 4000
#)
#precis(m_H1)
#precis(m_H1, depth = 3, pars = c("Rho", "Sigma"))
#post <- extract.samples(m_H1)
#a <- apply(post$a, 2, mean)
#b <- apply(post$b, 2, mean)
#plot(a, b, xlab = "a (intercept)", ylab = "b (urban slope)")
#abline(h = 0, lty = 2)
#abline(v = 0, lty = 2)
#R <- apply(post$Rho, 2:3, mean)
#s <- apply(post$Sigma, 2, mean)
#S <- diag(s) %*% R %*% diag(s)
#ll <- c(0.5, 0.67, 0.89, 0.97)
#for (l in ll) {
# el <- ellipse(S, centre = c(mean(post$abar), mean(post$bbar)), level = l)
#lines(el, col = "black", lwd = 0.5)
#}
#u0 <- inv_logit(a)
#u1 <- inv_logit(a + b)
#plot(u0, u1, xlim = c(0, 1), ylim = c(0, 1), xlab = "urban = 0", ylab = "urban = 1")
#abline(h = 0.5, lty = 2)
#abline(v = 0.5, lty = 2)Due error I had to comment out the code