TuringModels

Varying intercepts Chimpanzees

In Statistical Rethinking Edition 1, the model below is called m10.4.

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

Data

import CSV

using TuringModels: project_root
using DataFrames

path = joinpath(project_root, "data", "chimpanzees.csv")
df = CSV.read(path, DataFrame; delim=';');

Model

using Turing

@model m12_4(pulled_left, actor, condition, prosoc_left) = begin
    # Total num of y
    N = length(pulled_left)

    # Separate σ priors for each actor
    σ_actor ~ truncated(Cauchy(0, 1), 0, Inf)

    # Number of unique actors in the data set
    N_actor = length(unique(actor)) #7

    # Vector of actors (1,..,7) which we'll set priors on
    α_actor ~ filldist(Normal(0, σ_actor), N_actor)

    # Prior for intercept, prosoc_left, and the interaction
    α ~ Normal(0, 10)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logitp = α .+ α_actor[actor] .+
            (βp .+ βpC * condition) .* prosoc_left

    pulled_left .~ BinomialLogit.(1, logitp)
end;

Output

chns = sample(
    m12_4(df.pulled_left, df.actor, df.condition, df.prosoc_left),
    NUTS(),
    1000
)
Chains MCMC chain (1000×23×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 14.39 seconds
Compute duration  = 14.39 seconds
parameters        = α, α_actor[3], α_actor[7], α_actor[6], α_actor[1], α_actor[4], α_actor[2], σ_actor, α_actor[5], βp, βpC
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

     σ_actor    2.3283    0.9745     0.0308    0.0581   245.6740    0.9991       17.0773
  α_actor[1]   -1.2855    1.1520     0.0364    0.0773   147.0738    0.9990       10.2234
  α_actor[2]    4.1384    1.6451     0.0520    0.0937   384.4813    0.9990       26.7261
  α_actor[3]   -1.5974    1.1568     0.0366    0.0780   146.8539    0.9990       10.2081
  α_actor[4]   -1.5776    1.1561     0.0366    0.0770   148.2549    0.9991       10.3055
  α_actor[5]   -1.2923    1.1538     0.0365    0.0745   151.5998    0.9991       10.5380
  α_actor[6]   -0.3550    1.1604     0.0367    0.0778   148.9696    0.9991       10.3552
  α_actor[7]    1.1924    1.1689     0.0370    0.0761   158.0406    0.9990       10.9857
           α    0.5698    1.1367     0.0359    0.0762   146.5986    0.9992       10.1904
          βp    0.8303    0.2579     0.0082    0.0103   662.6818    1.0068       46.0644
         βpC   -0.1429    0.2992     0.0095    0.0124   641.0421    0.9998       44.5601

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

     σ_actor    1.1449    1.6583    2.1156    2.7348    4.8391
  α_actor[1]   -4.1733   -1.8091   -1.1774   -0.5477    0.7028
  α_actor[2]    1.4814    3.0265    3.9880    5.1005    7.9054
  α_actor[3]   -4.4279   -2.1524   -1.4758   -0.8923    0.2923
  α_actor[4]   -4.5354   -2.0799   -1.4515   -0.8297    0.4032
  α_actor[5]   -4.1551   -1.8282   -1.1719   -0.5754    0.6383
  α_actor[6]   -3.3066   -0.9104   -0.2197    0.3618    1.6376
  α_actor[7]   -1.7522    0.6478    1.2977    1.9200    3.1474
           α   -1.3467   -0.1549    0.4393    1.1132    3.3263
          βp    0.3387    0.6406    0.8301    1.0093    1.3005
         βpC   -0.7312   -0.3448   -0.1495    0.0536    0.4606

using StatsPlots

StatsPlots.plot(chns)

Original output

m124rethinking = "
             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
 a_actor[1]  -1.13   0.95      -2.62       0.27  2739    1
 a_actor[2]   4.17   1.66       1.80       6.39  3958    1
 a_actor[3]  -1.44   0.95      -2.90      -0.02  2720    1
 a_actor[4]  -1.44   0.94      -2.92      -0.04  2690    1
 a_actor[5]  -1.13   0.94      -2.58       0.31  2727    1
 a_actor[6]  -0.19   0.94      -1.64       1.25  2738    1
 a_actor[7]   1.34   0.97      -0.09       2.87  2889    1
 a            0.42   0.93      -1.00       1.81  2622    1
 bp           0.83   0.26       0.41       1.25  8594    1
 bpC         -0.13   0.30      -0.62       0.34  8403    1
 sigma_actor  2.26   0.94       1.07       3.46  4155    1
";