Statistical Rethinking 2: Chapter 12
Vincent Arel-Bundock
Status
Estimated and checked against the book:
- m12.2
- m12.3
- m12.4
- m12.5
- m12.6
Different results:
- m12.1
Libraries
library(tidyverse)
library(tidybayes)
library(rstan)
library(patchwork)
library(rethinking)
options(mc.cores = 4)Section 12.1: Over-dispersed counts
Warning: I have not yet been able to replicate the ulam results for m12.1 in rstan.
stan_program <- '
data {
int n;
int applications[n];
int admit[n];
int gender[n];
}
parameters {
vector[2] a;
real<lower=0> phi;
}
transformed parameters {
real theta;
theta = phi + 2;
}
model {
vector[n] pbar;
phi ~ exponential(1);
a ~ normal(0, 1.5);
for (i in 1:n) {
pbar[i] = a[gender[i]];
pbar[i] = inv_logit(pbar[i]);
}
admit ~ beta_binomial(applications, pbar * theta, (1 - pbar) * theta);
}
generated quantities {
real da;
da = a[1] - a[2];
}
'
stan_data <- read.csv('data/UCBadmit.csv', sep = ';') %>%
rename(gender = applicant.gender) %>%
mutate(gender = ifelse(gender == 'female', 2, 1)) %>%
compose_data
m12.1 <- stan(model_code = stan_program, data = stan_data, iter = 5000,
control = list(adapt_delta = .95))summary(m12.1, c('a', 'phi'))$summary## mean se_mean sd 2.5% 25% 50% 75%
## a[1] 0.03838894 0.9139069 1.3252213 -2.4734278 -0.5980206 0.23252638 0.9699306
## a[2] -0.24229634 1.0034895 1.4513773 -2.5466493 -1.4998133 0.06501208 1.1753202
## phi 1.35550266 0.4571877 0.9687508 0.2375496 0.5852012 0.99577804 2.0409857
## 97.5% n_eff Rhat
## a[1] 1.880505 2.102679 7.213213
## a[2] 1.409781 2.091871 5.906272
## phi 3.445442 4.489883 1.788069stan_data <- read.csv('data/Kline.csv', sep = ';') %>%
mutate(logpop = standardize(log(population))) %>%
compose_data(n_contact = n_distinct(contact))
stan_program <- "
data{
int total_tools[10];
int population[10];
int contact[10];
}
parameters{
vector[2] a;
vector<lower=0>[2] b;
real<lower=0> g;
real<lower=0> phi;
}
model{
vector[10] lambda;
phi ~ exponential( 1 );
g ~ exponential( 1 );
b ~ exponential( 1 );
a ~ normal( 1 , 1 );
for ( i in 1:10 ) {
lambda[i] = exp(a[contact[i]]) * population[i]^b[contact[i]]/g;
}
total_tools ~ neg_binomial_2( lambda , phi );
}
"
m12.2 <- stan(model_code = stan_program, data = stan_data)summary(m12.2, c('a', 'b', 'phi', 'g'))$summary## mean se_mean sd 2.5% 25% 50% 75%
## a[1] 1.0473810 0.021443773 0.94511985 -0.80883614 0.4340131 1.0365129 1.6850529
## a[2] 0.8963278 0.018600123 0.83486702 -0.71463409 0.3480698 0.9090939 1.4549253
## b[1] 0.2604337 0.003157677 0.13075554 0.02941905 0.1605860 0.2596623 0.3534109
## b[2] 0.2475842 0.002473933 0.09984088 0.05146194 0.1811474 0.2474376 0.3100577
## phi 3.6667090 0.039209148 1.67825227 1.18066075 2.4716135 3.4023163 4.6141558
## g 1.0611297 0.019256332 0.87058932 0.13563939 0.4560126 0.8137217 1.4010581
## 97.5% n_eff Rhat
## a[1] 2.9075997 1942.546 1.003003
## a[2] 2.5134734 2014.666 1.001967
## b[1] 0.5213364 1714.687 1.002796
## b[2] 0.4537079 1628.700 1.001702
## phi 7.5075099 1832.060 1.003350
## g 3.4023922 2043.994 1.000075Section 12.2: Zero-inflated outcomes
set.seed(365)
prob_drink <- .2
rate_work <- 1
n <- 365
drink <- rbinom(n, 1, prob_drink)
y <- (1 - drink) * rpois(n, rate_work)
stan_data <- list(y = y, n = length(y))
stan_program <- "
data {
int n;
int y[n];
}
parameters {
real ap;
real al;
}
model {
real p;
real lambda;
p = inv_logit(ap);
lambda = exp(al);
ap ~ normal(-1.5, 1);
al ~ normal(1, .5);
for ( i in 1:n ) {
if (y[i] == 0)
target += log_mix( p , 0 , poisson_lpmf(0 | lambda) );
if (y[i] > 0)
target += log1m( p ) + poisson_lpmf(y[i] | lambda );
}
}
"
m12.3 <- stan(model_code = stan_program, data = stan_data)summary(m12.3, c('al', 'ap'))$summary## mean se_mean sd 2.5% 25% 50%
## al 0.0109337 0.002489699 0.08774037 -0.1645018 -0.04919428 0.01217841
## ap -1.2783034 0.009979309 0.35100669 -2.0992802 -1.46883873 -1.23371035
## 75% 97.5% n_eff Rhat
## al 0.0720189 0.1778382 1241.953 1.0005837
## ap -1.0311268 -0.7270489 1237.171 0.9995458Section 12.3: Ordered categorical outcomes
stan_data <- read.csv('data/Trolley.csv', sep = ';') %>%
compose_data(k = n_distinct(response) - 1)
stan_program <- "
data {
int n;
int k;
int response[n];
}
parameters {
ordered[k] cutpoints;
}
model {
for (i in 1:n) {
response[i] ~ ordered_logistic(0, cutpoints);
}
cutpoints ~ normal(0, 15);
}
"
m12.4 <- stan(model_code = stan_program, data = stan_data)stan_data <- read.csv('data/Trolley.csv', sep = ';') %>%
compose_data(k = n_distinct(response) - 1)
stan_program <- "
data {
int n;
int k;
int response[n];
int action[n];
int intention[n];
int contact[n];
}
parameters {
ordered[k] cutpoints;
real bA;
real bI;
real bC;
real bIA;
real bIC;
}
model {
vector[n] BI;
vector[n] phi;
for (i in 1:n) {
BI[i] = bI + bIA * action[i] + bIC * contact[i];
phi[i] = bA * action[i] + bC * contact[i] + BI[i] * intention[i];
response[i] ~ ordered_logistic(phi[i], cutpoints);
}
cutpoints ~ normal(0, 15);
bA ~ normal(0, .5);
bI ~ normal(0, .5);
bC ~ normal(0, .5);
bIA ~ normal(0, .5);
bIC ~ normal(0, .5);
}
"
m12.5 <- stan(model_code = stan_program, data = stan_data)datplot <- m12.5 %>%
gather_draws(bA, bI, bC, bIA, bIC) %>%
mean_qi()
ggplot(datplot, aes(y = .variable, x = .value, xmin = .lower, xmax = .upper)) +
geom_pointrange() +
geom_vline(xintercept = 0, linetype = 'dotted')
Section 12.4: Ordered categorical predictors
stan_program <- '
data{
int n;
int response[n];
int contact[n];
int intention[n];
int action[n];
int edu_new[n];
vector[7] alpha;
}
parameters{
ordered[6] kappa;
real bE;
real bC;
real bI;
real bA;
simplex[7] delta;
}
model{
vector[n] phi;
vector[8] delta_j;
delta ~ dirichlet( alpha );
delta_j = append_row(0, delta);
bA ~ normal( 0 , 1 );
bI ~ normal( 0 , 1 );
bC ~ normal( 0 , 1 );
bE ~ normal( 0 , 1 );
kappa ~ normal( 0 , 1.5 );
for ( i in 1:n ) {
phi[i] = bE * sum(delta_j[1:edu_new[i]]) +
bA * action[i] +
bI * intention[i] +
bC * contact[i];
}
for ( i in 1:n ) response[i] ~ ordered_logistic( phi[i] , kappa );
}
'
stan_data <- read.csv('data/Trolley.csv', sep = ';') %>%
# ordered education levels
mutate(edu_new = factor(edu, levels = c("Elementary School",
"Middle School",
"Some High School",
"High School Graduate",
"Some College",
"Bachelor's Degree",
"Master's Degree",
"Graduate Degree"))) %>%
compose_data(alpha = rep(2, 7))
m12.6 <- stan(model_code = stan_program, data = stan_data)summary(m12.6, c('bE', 'bC', 'bI', 'bA', 'delta'))$summary## mean se_mean sd 2.5% 25%
## bE -0.31466345 0.0043757177 0.16834789 -0.675920028 -0.41764204
## bC -0.95724293 0.0008955266 0.05068350 -1.058917578 -0.99092317
## bI -0.71761851 0.0006484496 0.03751921 -0.792554903 -0.74289110
## bA -0.70586725 0.0007248147 0.04056351 -0.784698353 -0.73361235
## delta[1] 0.22572949 0.0029095910 0.13761260 0.026230428 0.11903511
## delta[2] 0.14284329 0.0014032959 0.08746129 0.020022979 0.07680229
## delta[3] 0.19475894 0.0018826557 0.10841821 0.033478261 0.11027181
## delta[4] 0.16942198 0.0016805074 0.09647420 0.024985339 0.09668394
## delta[5] 0.04352034 0.0012670003 0.05185546 0.003628401 0.01593960
## delta[6] 0.09814597 0.0011301656 0.06529187 0.014291132 0.05059393
## delta[7] 0.12557999 0.0011845177 0.07546800 0.019498693 0.06898606
## 50% 75% 97.5% n_eff Rhat
## bE -0.30516717 -0.20576328 0.006266416 1480.188 1.0017980
## bC -0.95692838 -0.92196319 -0.861506461 3203.142 0.9998047
## bI -0.71695203 -0.69218772 -0.645612145 3347.765 0.9997348
## bA -0.70596323 -0.67902477 -0.626449932 3131.966 0.9997095
## delta[1] 0.20266165 0.31337551 0.542845855 2236.930 0.9999760
## delta[2] 0.12641848 0.19081885 0.351299209 3884.484 1.0003031
## delta[3] 0.18036704 0.26367469 0.437345817 3316.368 0.9991919
## delta[4] 0.15650149 0.22718468 0.393230762 3295.654 1.0005179
## delta[5] 0.02888079 0.05198047 0.188774507 1675.080 1.0003732
## delta[6] 0.08472289 0.13159884 0.264602857 3337.598 0.9995975
## delta[7] 0.11090345 0.16833402 0.308565547 4059.220 0.9994046