Estimated and checked against the book:
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstan)
library(patchwork)
options(mc.cores = parallel::detectCores())
rugged <- read.csv('data/rugged.csv', sep = ';') %>%
mutate(log_gdp = log(rgdppc_2000),
log_gdp_std = log_gdp / mean(log_gdp, na.rm = TRUE),
rugged_std = rugged / max(rugged, na.rm = TRUE),
region = ifelse(cont_africa == 1, 'Africa', 'Not Africa')) %>%
select(country, log_gdp_std, rugged_std, region) %>%
drop_na
stan_data <- compose_data(rugged)
stan_program <- '
data {
int<lower=1> n; // number of observations
vector[n] log_gdp_std; // outcome
vector[n] rugged_std; // regressor
int region[n]; // africa indicator
}
parameters {
real<lower=0> sigma;
vector[2] a;
vector[2] b;
}
model {
vector[n] mu;
for (i in 1:n) {
mu[i] = a[region[i]] + b[region[i]] * (rugged_std[i] - 0.215);
}
a ~ normal(1, 0.1);
b ~ normal(0, 0.3);
sigma ~ exponential(1);
log_gdp_std ~ normal(mu, sigma);
}
'
m9.1 <- stan(model_code = stan_program, data = stan_data)
m9.1
## Inference for Stan model: b01cd93ce6fd8681789fca5db2e21b98.
## 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
## sigma 0.11 0.00 0.01 0.10 0.11 0.11 0.12 0.12 3988 1
## a[1] 0.89 0.00 0.02 0.86 0.88 0.89 0.91 0.93 3931 1
## a[2] 1.04 0.00 0.01 1.02 1.04 1.04 1.05 1.06 4456 1
## b[1] 0.17 0.00 0.09 -0.02 0.11 0.17 0.23 0.34 3965 1
## b[2] -0.18 0.00 0.07 -0.32 -0.23 -0.18 -0.13 -0.03 3969 1
## lp__ 285.08 0.04 1.65 281.06 284.19 285.42 286.33 287.25 1849 1
##
## Samples were drawn using NUTS(diag_e) at Sat Aug 1 13:47:34 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).
mcmc_trace(m9.1)