TuringModels

Estimate handedness for each Chimpanzee

This is model m10.4 in Statistical Rethinking Edition 1.

  1. Data
  2. Model
  3. Output
  4. Original output

Data

import CSV
import Random
import TuringModels

using DataFrames
using StatsFuns

Random.seed!(1)

data_path = joinpath(TuringModels.project_root, "data", "chimpanzees.csv")
df = CSV.read(data_path, DataFrame; delim=';')
first(df, 10)
10×8 DataFrame
 Row │ actor  recipient  condition  block  trial  prosoc_left  chose_prosoc  pulled_left
     │ Int64  InlineSt…  Int64      Int64  Int64  Int64        Int64         Int64
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     1  NA                 0      1      2            0             1            0
   2 │     1  NA                 0      1      4            0             0            1
   3 │     1  NA                 0      1      6            1             0            0
   4 │     1  NA                 0      1      8            0             1            0
   5 │     1  NA                 0      1     10            1             1            1
   6 │     1  NA                 0      1     12            1             1            1
   7 │     1  NA                 0      2     14            1             0            0
   8 │     1  NA                 0      2     16            1             0            0
   9 │     1  NA                 0      2     18            0             1            0
  10 │     1  NA                 0      2     20            0             1            0

Model

using Turing

@model m10_4(y, actors, x₁, x₂) = begin
    # Number of unique actors in the data set
    N_actor = length(unique(actors))

    # Set an TArray for the priors/param
    α ~ filldist(Normal(0, 10), N_actor)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logits = α[actors] .+ (βp .+ βpC * x₁) .* x₂
    y .~ BinomialLogit.(1, logits)
end

model = m10_4(df.pulled_left, df.actor, df.condition, df.prosoc_left);

Output

chns = sample(model, NUTS(), 1000)
Chains MCMC chain (1000×21×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 9.7 seconds
Compute duration  = 9.7 seconds
parameters        = βpC, α[5], α[1], α[7], α[2], α[6], α[3], βp, α[4]
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]   -0.7498    0.2578     0.0082    0.0089    939.6271    0.9995       96.8388
        α[2]   10.7558    5.5484     0.1755    0.2573    348.3636    1.0009       35.9027
        α[3]   -1.0478    0.2882     0.0091    0.0106    863.7634    0.9992       89.0202
        α[4]   -1.0596    0.2783     0.0088    0.0116    701.8918    0.9991       72.3376
        α[5]   -0.7458    0.2816     0.0089    0.0096    866.2697    0.9990       89.2785
        α[6]    0.2262    0.2810     0.0089    0.0103   1059.6272    0.9998      109.2061
        α[7]    1.7900    0.4084     0.0129    0.0133    998.6316    1.0002      102.9199
          βp    0.8536    0.2690     0.0085    0.0120    463.3350    1.0011       47.7517
         βpC   -0.1559    0.3072     0.0097    0.0112    743.4991    0.9995       76.6257

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

        α[1]   -1.2464   -0.9215   -0.7576   -0.5729   -0.2393
        α[2]    3.6650    6.5200    9.2989   13.7972   24.1867
        α[3]   -1.6345   -1.2293   -1.0499   -0.8612   -0.5014
        α[4]   -1.6078   -1.2424   -1.0646   -0.8763   -0.4969
        α[5]   -1.2867   -0.9452   -0.7551   -0.5536   -0.1999
        α[6]   -0.3196    0.0360    0.2302    0.4202    0.7655
        α[7]    1.0326    1.5107    1.7936    2.0513    2.6275
          βp    0.3225    0.6794    0.8513    1.0334    1.3902
         βpC   -0.7461   -0.3623   -0.1669    0.0490    0.4916

using StatsPlots

StatsPlots.plot(chns)

Original output

m_10_04s_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
        Mean        SD       Naive SE       MCSE      ESS
a.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000
a.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000
a.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000
a.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000
a.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000
a.6  0.21599365 0.26307574 0.0041595927 0.0045153523 1000
a.7  1.81090866 0.39318577 0.0062168129 0.0071483527 1000
bp  0.83979926 0.26284676 0.0041559722 0.0059795826 1000
bpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000
";