Africa first candidate model

In Edition 2 on page 244, McElreath defines the model as

To define this in Turing.jl, we change it to

where $y_i$ is the GDPs for nation $i$, $r_i$ is terrain ruggedness for nation $i$ and $\overline{r}$ is the mean of the ruggedness in the whole sample.

McElreath calls two models m8.1. We have tried to get close to the precis output on page 245, but both m8.1 seem to give a different output for $\sigma$. For \sigma and \beta, our results on this page correspond to the resutls by McElreath.

Data

using DataFrames
using Turing
using TuringModels

import CSV

data_path = joinpath(TuringModels.project_root, "data", "rugged.csv")
df = CSV.read(data_path, DataFrame)

df.log_gdp = log.(df.rgdppc_2000)
dropmissing!(df)

df = select(df, :log_gdp, :rugged, :cont_africa);

df.log_gdp_std = df.log_gdp ./ mean(df.log_gdp)
df.rugged_std = df.rugged ./ maximum(df.rugged)

first(df, 8)

8×5 DataFrame
 Row │ log_gdp  rugged   cont_africa  log_gdp_std  rugged_std
     │ Float64  Float64  Int64        Float64      Float64
─────┼────────────────────────────────────────────────────────
   1 │ 7.49261    0.858            1     1.00276    0.138342
   2 │ 6.43238    1.78             1     0.860869   0.287004
   3 │ 6.86612    0.141            1     0.918918   0.0227346
   4 │ 6.90617    0.236            1     0.924278   0.0380522
   5 │ 8.9493     0.181            1     1.19772    0.0291841
   6 │ 7.04582    0.197            1     0.942968   0.0317639
   7 │ 7.36241    0.224            1     0.985338   0.0361174
   8 │ 7.54046    0.515            1     1.00917    0.0830377

Model

@model function model_fn(log_gdp_std, rugged_std, mean_rugged)
    α ~ Normal(1, 0.1)
    β ~ Normal(0, 0.3)
    σ ~ Exponential(1)

    μ = α .+ β * (rugged_std .- mean_rugged)
    log_gdp_std ~ MvNormal(μ, σ)
end

model = model_fn(df.log_gdp, df.rugged_std, mean(df.rugged_std));

Output

chns = sample(model, NUTS(), MCMCThreads(), 1000, 3)

Chains MCMC chain (1000×15×3 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 3
Samples per chain = 1000
Wall duration     = 6.07 seconds
Compute duration  = 5.56 seconds
parameters        = α, β, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           α    1.0768    0.0972     0.0018    0.0012   4781.0065    0.9993      859.7386
           β    0.0067    0.2991     0.0055    0.0053   4069.8837    1.0005      731.8618
           σ    6.1882    0.6058     0.0111    0.0093   3881.9884    0.9992      698.0738

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

           α    0.8866    1.0081    1.0770    1.1441    1.2659
           β   -0.5781   -0.2036    0.0003    0.2146    0.5859
           σ    5.1186    5.7773    6.1448    6.5704    7.5015

using StatsPlots

StatsPlots.plot(chns)

"/home/runner/work/TuringModels.jl/TuringModels.jl/__site/assets/models/africa-first-candidate/code/output/chns.svg"

Original output

From page 245: