TuringModels

Multi-multilevel Chimpanzees

In Statistical Rethinking Edition 1, McElreath adds yet another varying intercept, adaptive prior and standard deviation parameter to the Chimpanzees model.

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

Data

import CSV
import Random

using TuringModels: project_root
using DataFrames

Random.seed!(1)

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

Model

using Turing

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

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

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

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

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

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

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

Output

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

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 29.17 seconds
Compute duration  = 29.17 seconds
parameters        = α_block[2], α, α_actor[3], α_actor[7], α_actor[6], α_block[5], α_block[6], α_actor[1], α_actor[4], α_actor[2], σ_actor, α_actor[5], α_block[3], α_block[4], βp, α_block[1], σ_block, β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.2063    0.8334     0.0264    0.0341   386.9044    0.9991       13.2647
     σ_block    0.2146    0.1659     0.0052    0.0132   116.7199    0.9999        4.0016
  α_actor[1]   -1.1679    0.9599     0.0304    0.0726   188.5818    1.0016        6.4654
  α_actor[2]    3.9924    1.3788     0.0436    0.0593   496.6002    1.0019       17.0255
  α_actor[3]   -1.4702    0.9563     0.0302    0.0748   175.1971    1.0018        6.0065
  α_actor[4]   -1.4671    0.9488     0.0300    0.0742   177.6456    1.0012        6.0904
  α_actor[5]   -1.1680    0.9488     0.0300    0.0726   182.2749    1.0026        6.2491
  α_actor[6]   -0.2080    0.9682     0.0306    0.0739   178.3788    1.0020        6.1156
  α_actor[7]    1.3083    0.9661     0.0306    0.0696   199.2704    1.0009        6.8318
  α_block[1]   -0.1769    0.2224     0.0070    0.0153   163.5154    0.9992        5.6060
  α_block[2]    0.0365    0.1775     0.0056    0.0068   539.6219    1.0002       18.5005
  α_block[3]    0.0474    0.1702     0.0054    0.0054   702.6173    1.0004       24.0886
  α_block[4]    0.0057    0.1745     0.0055    0.0052   681.0657    0.9990       23.3498
  α_block[5]   -0.0429    0.1823     0.0058    0.0077   490.2252    0.9990       16.8070
  α_block[6]    0.0991    0.1769     0.0056    0.0083   384.0458    1.0004       13.1667
           α    0.4499    0.9450     0.0299    0.0739   178.4373    1.0017        6.1176
          βp    0.8373    0.2742     0.0087    0.0108   576.6973    0.9990       19.7716
         βpC   -0.1504    0.2990     0.0095    0.0111   697.9459    0.9991       23.9285

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

     σ_actor    1.0856    1.6356    2.0416    2.5556    4.3793
     σ_block    0.0290    0.0909    0.1745    0.2873    0.6438
  α_actor[1]   -3.3331   -1.6389   -1.0905   -0.6026    0.4973
  α_actor[2]    1.7741    3.0285    3.8251    4.7674    7.1154
  α_actor[3]   -3.7439   -1.9875   -1.3871   -0.8809    0.1595
  α_actor[4]   -3.5677   -1.9728   -1.3856   -0.8786    0.1741
  α_actor[5]   -3.3915   -1.6409   -1.0903   -0.6009    0.5588
  α_actor[6]   -2.3817   -0.7160   -0.1269    0.3977    1.5071
  α_actor[7]   -0.8393    0.7927    1.3498    1.9035    3.0654
  α_block[1]   -0.7243   -0.2842   -0.1200   -0.0235    0.1089
  α_block[2]   -0.3261   -0.0516    0.0226    0.1257    0.4158
  α_block[3]   -0.2804   -0.0482    0.0311    0.1318    0.4348
  α_block[4]   -0.3866   -0.0761    0.0026    0.0929    0.3591
  α_block[5]   -0.4616   -0.1365   -0.0262    0.0578    0.2908
  α_block[6]   -0.1994   -0.0060    0.0648    0.1933    0.5114
           α   -1.1867   -0.1241    0.3824    0.9472    2.6300
          βp    0.3102    0.6468    0.8314    1.0225    1.4223
         βpC   -0.7350   -0.3546   -0.1429    0.0487    0.4644

using StatsPlots

StatsPlots.plot(chns)

Original output

m125rethinking = "
             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
a_actor[1]  -1.19   0.98      -2.66       0.35  3451    1
a_actor[2]   4.14   1.59       1.87       6.38  5514    1
a_actor[3]  -1.49   0.99      -2.96       0.05  3493    1
a_actor[4]  -1.50   0.98      -3.00       0.01  3340    1
a_actor[5]  -1.19   0.99      -2.68       0.34  3447    1
a_actor[6]  -0.25   0.99      -1.79       1.25  3361    1
a_actor[7]   1.30   1.01      -0.21       2.89  3673    1
a_block[1]  -0.19   0.23      -0.56       0.10  3825    1
a_block[2]   0.04   0.19      -0.24       0.34  9346    1
a_block[3]   0.06   0.19      -0.22       0.37  9202    1
a_block[4]   0.00   0.18      -0.29       0.29 11314    1
a_block[5]  -0.03   0.19      -0.33       0.25 10596    1
a_block[6]   0.11   0.21      -0.16       0.45  6040    1
a            0.47   0.97      -1.01       1.99  3273    1
bp           0.83   0.26       0.40       1.24 13225    1
bpc         -0.15   0.30      -0.61       0.36  8492    1
sigma_actor  2.27   0.91       1.03       3.35  5677    1
sigma_block  0.23   0.18       0.01       0.44  2269    1
";