Statistical Rethinking 2: Chapter 6
Vincent Arel-Bundock
Status
Estimated and checked against the book:
- m6.1
- m6.2
- m6.3
- m6.4
- m6.5
- m6.6
- m6.7
- m6.1
- m6.8
- m6.9
- m6.10
- m6.11
- m6.12
Libraries
library(tidyverse)
library(tidybayes)
library(rstan)
library(patchwork)
options(mc.cores = 4)Section 6.1: Multicollinearity
set.seed(909)
N <- 100
stan_data <- tibble(height = rnorm(N,10,2),
leg_prop = runif(N,0.4,0.5),
leg_left = leg_prop*height + rnorm( N , 0 , 0.02 ),
leg_right = leg_prop*height + rnorm( N , 0 , 0.02 )) %>%
compose_data
stan_program <- '
data {
int n;
vector[n] height;
vector[n] leg_left;
vector[n] leg_right;
}
parameters {
real a;
real bl;
real br;
real<lower=0> sigma;
}
model {
vector[n] mu;
mu = a + bl * leg_left + br * leg_right;
height ~ normal(mu, sigma);
a ~ normal(10, 100);
bl ~ normal(2, 10);
br ~ normal(2, 10);
sigma ~ exponential(1);
}
'
m6.1 <- stan(model_code = stan_program, data = stan_data)m6.1## Inference for Stan model: 6a21a0fb5e106f7ac61d55274e661840.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a 0.97 0.01 0.30 0.36 0.78 0.97 1.17 1.55 2004 1
## bl 0.25 0.06 2.57 -4.85 -1.45 0.29 2.00 5.23 1650 1
## br 1.75 0.06 2.58 -3.25 -0.01 1.70 3.44 6.82 1648 1
## sigma 0.64 0.00 0.05 0.55 0.60 0.63 0.66 0.73 2293 1
## lp__ -5.22 0.04 1.52 -9.08 -5.91 -4.85 -4.11 -3.41 1304 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:50:15 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).# Figure 6.2
datplot <- m6.1 %>%
spread_draws(bl, br) %>%
mutate(bl_br_sum = bl + br)
p1 <- ggplot(datplot, aes(br, bl)) + geom_point(alpha = .1)
p2 <- ggplot(datplot, aes(bl_br_sum)) + geom_density()
p1 + p2
stan_program <- '
data {
int n;
vector[n] height;
vector[n] leg_left;
vector[n] leg_right;
}
parameters {
real a;
real bl;
real<lower=0> sigma;
}
model {
vector[n] mu;
mu = a + bl * leg_left;
height ~ normal(mu, sigma);
a ~ normal(10, 100);
bl ~ normal(2, 10);
sigma ~ exponential(1);
}
'
m6.2 <- stan(model_code = stan_program, data = stan_data)m6.2## Inference for Stan model: 9894da73a0f125d7189dd60c63773ccc.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a 0.98 0.01 0.29 0.42 0.79 0.99 1.18 1.56 1234 1
## bl 2.00 0.00 0.06 1.87 1.95 1.99 2.04 2.12 1276 1
## sigma 0.63 0.00 0.05 0.55 0.60 0.63 0.66 0.73 1722 1
## lp__ -4.87 0.03 1.20 -7.91 -5.41 -4.58 -3.99 -3.48 1585 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:50:40 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).dat <- read.csv('data/milk.csv', sep = ';') %>%
mutate(K = as.vector(scale(kcal.per.g)),
F = as.vector(scale(perc.fat)),
L = as.vector(scale(perc.lactose)))
stan_program <- '
data {
int n;
int k;
matrix[n, k] X;
vector[n] y;
}
parameters {
vector[k] b;
real<lower=0> sigma;
}
model {
vector[n] mu;
mu = X * b;
y ~ normal(mu, sigma);
b ~ normal(0, 0.5);
sigma ~ exponential(1);
}
'
stan_data <- compose_data(dat,
y = K,
X = model.matrix(~F, dat),
k = ncol(X))
m6.3 <- stan(model_code = stan_program, data = stan_data)
stan_data <- compose_data(dat,
y = K,
X = model.matrix(~L, dat),
k = ncol(X))
m6.4 <- stan(model_code = stan_program, data = stan_data)
stan_data <- compose_data(dat,
y = K,
X = model.matrix(~F + L, dat),
k = ncol(X))
m6.5 <- stan(model_code = stan_program, data = stan_data)m6.3## Inference for Stan model: 57d226675d4bd9c6ae3ce51f63a1d8e3.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## b[1] 0.00 0.00 0.09 -0.19 -0.06 0.00 0.06 0.18 3570 1
## b[2] 0.86 0.00 0.09 0.66 0.80 0.86 0.92 1.04 3227 1
## sigma 0.49 0.00 0.07 0.37 0.44 0.48 0.53 0.65 2732 1
## lp__ 4.02 0.04 1.38 0.36 3.45 4.38 5.02 5.54 1558 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:51:07 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).m6.4## Inference for Stan model: 57d226675d4bd9c6ae3ce51f63a1d8e3.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## b[1] 0.00 0.00 0.08 -0.15 -0.05 0.00 0.05 0.15 3649 1
## b[2] -0.90 0.00 0.08 -1.05 -0.95 -0.90 -0.85 -0.75 3482 1
## sigma 0.41 0.00 0.06 0.32 0.37 0.41 0.45 0.54 2707 1
## lp__ 8.81 0.03 1.25 5.58 8.25 9.13 9.73 10.26 1673 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:51:07 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).m6.5## Inference for Stan model: 57d226675d4bd9c6ae3ce51f63a1d8e3.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## b[1] 0.00 0.00 0.08 -0.14 -0.05 0.00 0.05 0.15 2911 1
## b[2] 0.25 0.00 0.20 -0.13 0.12 0.25 0.38 0.65 1893 1
## b[3] -0.67 0.00 0.20 -1.06 -0.80 -0.67 -0.54 -0.28 1939 1
## sigma 0.41 0.00 0.06 0.32 0.37 0.41 0.45 0.54 2282 1
## lp__ 9.15 0.04 1.47 5.51 8.39 9.48 10.22 11.01 1437 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:51:08 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Section 6.2: Post-treatment bias
## R code 6.13
set.seed(71)
# number of plants
N <- 100
# simulate initial heights
h0 <- rnorm(N,10,2)
# assign treatments and simulate fungus and growth
treatment <- rep( 0:1 , each=N/2 )
fungus <- rbinom( N , size=1 , prob=0.5 - treatment*0.4 )
h1 <- h0 + rnorm(N, 5 - 3*fungus)
# compose a clean data frame
d <- data.frame( h0=h0 , h1=h1 , treatment=treatment , fungus=fungus )
stan_data <- compose_data(d)
stan_program <- '
data {
int n;
vector[n] h1;
vector[n] h0;
}
parameters {
real p;
real<lower=0> sigma;
}
model {
vector[n] mu;
mu = h0 * p;
h1 ~ normal(mu, sigma);
sigma ~ exponential(1);
p ~ lognormal(0, 0.25);
}
'
m6.6 <- stan(model_code = stan_program, data = stan_data)m6.6## Inference for Stan model: c3b555e32b2b63ee697970fdcfef4113.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## p 1.43 0.00 0.02 1.39 1.41 1.43 1.44 1.46 3335 1
## sigma 1.83 0.00 0.13 1.61 1.74 1.82 1.91 2.09 2690 1
## lp__ -112.83 0.02 0.98 -115.47 -113.20 -112.53 -112.14 -111.88 2137 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:51:31 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).stan_program <- '
data {
int n;
vector[n] h1;
vector[n] h0;
vector[n] treatment;
vector[n] fungus;
}
parameters {
real bt;
real bf;
real a;
real<lower=0> sigma;
}
model {
vector[n] mu;
vector[n] p;
p = a + bt * treatment + bf * fungus;
mu = h0 .* p;
h1 ~ normal(mu, sigma);
a ~ lognormal(0, 0.25);
bt ~ normal(0, 0.5);
bf ~ normal(0, 0.5);
sigma ~ exponential(1);
}
'
m6.7 <- stan(model_code = stan_program, data = stan_data)m6.7## Inference for Stan model: ad0eb54949ddc84fc0c13bafb5bcabf9.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## bt 0.00 0.00 0.03 -0.06 -0.02 0.00 0.02 0.06 1736 1
## bf -0.27 0.00 0.04 -0.34 -0.29 -0.27 -0.24 -0.20 2137 1
## a 1.48 0.00 0.03 1.43 1.47 1.48 1.50 1.53 1654 1
## sigma 1.45 0.00 0.11 1.25 1.37 1.44 1.51 1.68 2826 1
## lp__ -89.85 0.03 1.41 -93.39 -90.56 -89.53 -88.80 -88.09 1822 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:51:54 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).stan_program <- '
data {
int n;
vector[n] h1;
vector[n] h0;
vector[n] treatment;
vector[n] fungus;
}
parameters {
real bt;
real a;
real<lower=0> sigma;
}
model {
vector[n] mu;
vector[n] p;
p = a + bt * treatment;
mu = h0 .* p;
h1 ~ normal(mu, sigma);
a ~ lognormal(0, 0.25);
bt ~ normal(0, 0.5);
sigma ~ exponential(1);
}
'
m6.8 <- stan(model_code = stan_program, data = stan_data)m6.8## Inference for Stan model: 6184ea605a71870fc61fbbc528db9e32.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## bt 0.08 0.00 0.03 0.02 0.06 0.08 0.11 0.15 1936 1
## a 1.38 0.00 0.03 1.33 1.36 1.38 1.40 1.43 2020 1
## sigma 1.79 0.00 0.13 1.56 1.70 1.78 1.87 2.06 2508 1
## lp__ -110.47 0.03 1.19 -113.54 -111.02 -110.17 -109.59 -109.09 1673 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:52:30 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Section 6.3: Collider bias
# remotes::install_github('rmcelreath/rethinking')
library(rethinking)
stan_data <- sim_happiness(seed = 1977, N_years = 1000) %>%
filter(age > 17) %>%
mutate(age = (age - 18) / (65 - 18),
married = married + 1) %>%
compose_data
stan_program <- '
data {
int n;
vector[n] happiness;
vector[n] age;
int married[n];
}
parameters {
real a[2];
real bA;
real sigma;
}
model {
vector[n] mu;
for (i in 1:n) {
mu[i] = a[married[i]] + bA * age[i];
}
happiness ~ normal(mu, sigma);
a ~ normal(0, 1);
bA ~ normal(0, 2);
sigma ~ exponential(1);
}
'
m6.9 <- stan(model_code = stan_program, data = stan_data)m6.9## Inference for Stan model: e94e7c77147afa14e7669ca2febb0e9f.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a[1] -0.23 0.00 0.06 -0.36 -0.28 -0.24 -0.19 -0.11 1475 1
## a[2] 1.26 0.00 0.08 1.09 1.20 1.26 1.32 1.43 1500 1
## bA -0.75 0.00 0.11 -0.97 -0.83 -0.75 -0.67 -0.53 1332 1
## sigma 0.99 0.00 0.02 0.95 0.98 0.99 1.01 1.04 2630 1
## lp__ -474.82 0.04 1.46 -478.48 -475.53 -474.48 -473.77 -473.01 1705 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:52:59 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).stan_program <- '
data {
int n;
vector[n] happiness;
vector[n] age;
int married[n];
}
parameters {
real a;
real bA;
real sigma;
}
model {
vector[n] mu;
for (i in 1:n) {
mu[i] = a + bA * age[i];
}
happiness ~ normal(mu, sigma);
a ~ normal(0, 1);
bA ~ normal(0, 2);
sigma ~ exponential(1);
}
'
m6.10 <- stan(model_code = stan_program, data = stan_data)m6.10## Inference for Stan model: 66129038c498a75608907a5db88dfcd0.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a 0.00 0.00 0.07 -0.14 -0.05 0.00 0.05 0.15 1583 1
## bA -0.01 0.00 0.13 -0.26 -0.09 -0.01 0.08 0.24 1515 1
## sigma 1.22 0.00 0.03 1.16 1.20 1.22 1.24 1.27 2371 1
## lp__ -668.86 0.03 1.19 -671.89 -669.46 -668.57 -667.97 -667.47 1414 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:53:21 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).## R code 6.25
N <- 200 # number of grandparent-parent-child triads
b_GP <- 1 # direct effect of G on P
b_GC <- 0 # direct effect of G on C
b_PC <- 1 # direct effect of P on C
b_U <- 2 # direct effect of U on P and C
## R code 6.26
set.seed(1)
U <- 2*rbern( N , 0.5 ) - 1
G <- rnorm( N )
P <- rnorm( N , b_GP*G + b_U*U )
C <- rnorm( N , b_PC*P + b_GC*G + b_U*U )
d <- data.frame( C=C , P=P , G=G , U=U )
stan_data <- compose_data(d)
stan_program <- '
data {
int n;
vector[n] C;
vector[n] P;
vector[n] G;
}
parameters {
real a;
real b_PC;
real b_GC;
real<lower=0> sigma;
}
model {
vector[n] mu;
mu = a + b_PC * P + b_GC * G;
C ~ normal(mu, sigma);
sigma ~ exponential(1);
b_PC ~ normal(0, 1);
b_GC ~ normal(0, 1);
}
'
m6.11 <- stan(model_code = stan_program, data = stan_data)
stan_program <- '
data {
int n;
vector[n] C;
vector[n] P;
vector[n] G;
vector[n] U;
}
parameters {
real a;
real b_PC;
real b_GC;
real b_U;
real<lower=0> sigma;
}
model {
vector[n] mu;
mu = a + b_PC * P + b_GC * G + b_U * U;
C ~ normal(mu, sigma);
sigma ~ exponential(1);
b_PC ~ normal(0, 1);
b_GC ~ normal(0, 1);
b_U ~ normal(0, 1);
}
'
m6.12 <- stan(model_code = stan_program, data = stan_data)m6.11## Inference for Stan model: 05d80ceefb44ad92d033a1a4fe0aa6be.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a -0.12 0.00 0.10 -0.32 -0.19 -0.12 -0.05 0.08 3880 1
## b_PC 1.79 0.00 0.05 1.70 1.76 1.79 1.82 1.88 3722 1
## b_GC -0.84 0.00 0.11 -1.04 -0.91 -0.84 -0.77 -0.62 3860 1
## sigma 1.43 0.00 0.07 1.30 1.38 1.43 1.48 1.58 3717 1
## lp__ -174.38 0.03 1.44 -177.94 -175.08 -174.07 -173.33 -172.60 1847 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:53:42 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).m6.12## Inference for Stan model: 63f14fe8921051e80ae787d3b155dee6.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a -0.12 0.00 0.07 -0.27 -0.17 -0.12 -0.07 0.02 2946 1
## b_PC 1.01 0.00 0.07 0.88 0.97 1.01 1.06 1.15 1453 1
## b_GC -0.04 0.00 0.10 -0.24 -0.11 -0.04 0.03 0.15 1757 1
## b_U 1.99 0.00 0.16 1.68 1.89 2.00 2.10 2.29 1619 1
## sigma 1.04 0.00 0.05 0.94 1.00 1.03 1.07 1.14 3151 1
## lp__ -110.47 0.04 1.67 -114.49 -111.34 -110.10 -109.25 -108.30 1588 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 11:54:05 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).