Data
import CSV
import Random
import TuringModels
using DataFrames
Random.seed!(1)
data_path = joinpath(TuringModels.project_root, "data", "UCBadmit.csv")
df = CSV.read(data_path, DataFrame; delim=';')
dept_map = Dict(key => idx for (idx, key) in enumerate(unique(df.dept)))
df.male = [g == "male" ? 1 : 0 for g in df.gender]
df.dept_id = [dept_map[de] for de in df.dept]
df12×7 DataFrame
Row │ dept gender admit reject applications male dept_id
│ String1 String7 Int64 Int64 Int64 Int64 Int64
─────┼───────────────────────────────────────────────────────────────
1 │ A male 512 313 825 1 1
2 │ A female 89 19 108 0 1
3 │ B male 353 207 560 1 2
4 │ B female 17 8 25 0 2
5 │ C male 120 205 325 1 3
6 │ C female 202 391 593 0 3
7 │ D male 138 279 417 1 4
8 │ D female 131 244 375 0 4
9 │ E male 53 138 191 1 5
10 │ E female 94 299 393 0 5
11 │ F male 22 351 373 1 6
12 │ F female 24 317 341 0 6Model
using Turing
@model function m13_3(applications, dept_id, male, admit)
Rho ~ LKJ(2, 2.)
sigma_dept ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 2)
bm ~ Normal(0, 1)
a ~ Normal(0, 1)
dist_mu = [a, bm]
dist_Sigma = sigma_dept .* Rho .* sigma_dept'
dist_Sigma = (dist_Sigma' + dist_Sigma) / 2
a_bm_dept ~ filldist(MvNormal(dist_mu, dist_Sigma), 6)
a_dept, bm_dept = a_bm_dept[1, :], a_bm_dept[2, :]
logit_p = a_dept[dept_id] + bm_dept[dept_id] .* male
admit .~ BinomialLogit.(applications, logit_p)
end;Output
chns = sample(
m13_3(df.applications, df.dept_id, df.male, df.admit),
NUTS(0.95),
1000
)Chains MCMC chain (1000×32×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 79.51 seconds
Compute duration = 79.51 seconds
parameters = Rho[1,1], Rho[2,1], Rho[1,2], Rho[2,2], sigma_dept[1], sigma_dept[2], bm, a, a_bm_dept[1,1], a_bm_dept[2,1], a_bm_dept[1,2], a_bm_dept[2,2], a_bm_dept[1,3], a_bm_dept[2,3], a_bm_dept[1,4], a_bm_dept[2,4], a_bm_dept[1,5], a_bm_dept[2,5], a_bm_dept[1,6], a_bm_dept[2,6]
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
Rho[1,1] 1.0000 0.0000 0.0000 0.0000 NaN NaN NaN
Rho[2,1] -0.3295 0.3516 0.0111 0.0109 1096.3894 0.9992 13.7898
Rho[1,2] -0.3295 0.3516 0.0111 0.0109 1096.3894 0.9992 13.7898
Rho[2,2] 1.0000 0.0000 0.0000 0.0000 840.1684 0.9990 10.5672
sigma_dept[1] 1.6026 0.5714 0.0181 0.0232 849.4314 1.0052 10.6837
sigma_dept[2] 0.4926 0.2569 0.0081 0.0112 489.0014 1.0049 6.1504
bm -0.1773 0.2115 0.0067 0.0066 847.1247 1.0037 10.6547
a -0.3708 0.5502 0.0174 0.0155 863.3647 0.9995 10.8590
a_bm_dept[1,1] 1.3006 0.2496 0.0079 0.0071 835.2350 1.0056 10.5052
a_bm_dept[2,1] -0.7874 0.2661 0.0084 0.0080 800.1063 1.0075 10.0633
a_bm_dept[1,2] 0.7334 0.3216 0.0102 0.0114 885.6146 1.0028 11.1388
a_bm_dept[2,2] -0.1994 0.3248 0.0103 0.0108 888.7963 1.0020 11.1788
a_bm_dept[1,3] -0.6445 0.0838 0.0026 0.0025 960.8583 0.9992 12.0852
a_bm_dept[2,3] 0.0756 0.1394 0.0044 0.0038 1023.2940 0.9992 12.8705
a_bm_dept[1,4] -0.6208 0.1057 0.0033 0.0037 879.5548 0.9991 11.0626
a_bm_dept[2,4] -0.0871 0.1416 0.0045 0.0049 960.8747 0.9990 12.0854
a_bm_dept[1,5] -1.1315 0.1114 0.0035 0.0040 1090.1563 1.0006 13.7115
a_bm_dept[2,5] 0.1196 0.1872 0.0059 0.0063 1171.5534 1.0002 14.7352
a_bm_dept[1,6] -2.6062 0.1989 0.0063 0.0060 1078.6007 0.9996 13.5661
a_bm_dept[2,6] -0.1113 0.2622 0.0083 0.0084 1033.4306 1.0022 12.9980
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Rho[1,1] 1.0000 1.0000 1.0000 1.0000 1.0000
Rho[2,1] -0.8919 -0.5949 -0.3796 -0.0801 0.4383
Rho[1,2] -0.8919 -0.5949 -0.3796 -0.0801 0.4383
Rho[2,2] 1.0000 1.0000 1.0000 1.0000 1.0000
sigma_dept[1] 0.8746 1.2040 1.4763 1.8792 3.0716
sigma_dept[2] 0.1712 0.3295 0.4364 0.5989 1.1279
bm -0.6193 -0.2950 -0.1709 -0.0507 0.2259
a -1.4768 -0.7121 -0.3825 -0.0478 0.7984
a_bm_dept[1,1] 0.8021 1.1325 1.2949 1.4697 1.7812
a_bm_dept[2,1] -1.2934 -0.9645 -0.7873 -0.6045 -0.2592
a_bm_dept[1,2] 0.0926 0.5263 0.7311 0.9415 1.3598
a_bm_dept[2,2] -0.8271 -0.4087 -0.2045 0.0032 0.4698
a_bm_dept[1,3] -0.8163 -0.6986 -0.6469 -0.5866 -0.4777
a_bm_dept[2,3] -0.2171 -0.0156 0.0799 0.1665 0.3466
a_bm_dept[1,4] -0.8242 -0.6901 -0.6212 -0.5506 -0.4060
a_bm_dept[2,4] -0.3617 -0.1864 -0.0829 0.0059 0.1825
a_bm_dept[1,5] -1.3665 -1.2082 -1.1308 -1.0521 -0.9293
a_bm_dept[2,5] -0.2412 -0.0084 0.1173 0.2368 0.4879
a_bm_dept[1,6] -3.0090 -2.7346 -2.5999 -2.4719 -2.2297
a_bm_dept[2,6] -0.6587 -0.2760 -0.1133 0.0684 0.3981
using StatsPlots
StatsPlots.plot(chns)"/home/runner/work/TuringModels.jl/TuringModels.jl/__site/assets/models/varying-intercepts-admission/code/output/chns.svg" 
Original output
m_13_3_rethinking = """
Inference for Stan model: f0d86ec689cbf7921aab4fc0f55616d2.
4 chains, each with iter=5000; warmup=1000; thin=1;
post-warmup draws per chain=4000, total post-warmup draws=16000.
mean se_mean sd 2.5% 25% 50% 75%
bm_dept[1] -0.79 0.00 0.27 -1.32 -0.97 -0.79 -0.61
bm_dept[2] -0.21 0.00 0.33 -0.87 -0.42 -0.20 0.00
bm_dept[3] 0.08 0.00 0.14 -0.18 -0.01 0.08 0.18
bm_dept[4] -0.09 0.00 0.14 -0.37 -0.19 -0.09 0.00
bm_dept[5] 0.12 0.00 0.19 -0.24 -0.01 0.12 0.25
bm_dept[6] -0.12 0.00 0.27 -0.67 -0.29 -0.12 0.05
a_dept[1] 1.30 0.00 0.25 0.82 1.13 1.30 1.47
a_dept[2] 0.74 0.00 0.33 0.08 0.52 0.73 0.95
a_dept[3] -0.65 0.00 0.09 -0.82 -0.71 -0.65 -0.59
a_dept[4] -0.62 0.00 0.10 -0.83 -0.69 -0.62 -0.55
a_dept[5] -1.13 0.00 0.11 -1.36 -1.21 -1.13 -1.05
a_dept[6] -2.60 0.00 0.20 -3.01 -2.73 -2.60 -2.46
a -0.49 0.01 0.73 -1.96 -0.91 -0.49 -0.07
bm -0.16 0.00 0.24 -0.65 -0.29 -0.16 -0.02
sigma_dept[1] 1.67 0.01 0.63 0.88 1.25 1.53 1.94
sigma_dept[2] 0.50 0.00 0.26 0.16 0.33 0.45 0.61
Rho[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00
Rho[1,2] -0.31 0.00 0.35 -0.87 -0.59 -0.35 -0.08
Rho[2,1] -0.31 0.00 0.35 -0.87 -0.59 -0.35 -0.08
Rho[2,2] 1.00 0.00 0.00 1.00 1.00 1.00 1.00
lp__ -2593.92 0.04 3.17 -2601.03 -2595.85 -2593.57 -2591.65
97.5% n_eff Rhat
bm_dept[1] -0.27 6918 1
bm_dept[2] 0.46 9133 1
bm_dept[3] 0.36 13016 1
bm_dept[4] 0.19 12258 1
bm_dept[5] 0.49 12953 1
bm_dept[6] 0.40 12011 1
a_dept[1] 1.81 6993 1
a_dept[2] 1.39 9335 1
a_dept[3] -0.48 13489 1
a_dept[4] -0.42 12595 1
a_dept[5] -0.91 14659 1
a_dept[6] -2.23 13168 1
a 0.96 9541 1
bm 0.31 8946 1
sigma_dept[1] 3.27 9058 1
sigma_dept[2] 1.15 6626 1
Rho[1,1] 1.00 NaN NaN
Rho[1,2] 0.44 11358 1
Rho[2,1] 0.44 11358 1
Rho[2,2] 1.00 16069 1
lp__ -2588.70 5310 1
""";
CC BY-SA 4.0
Rob Goedman, Richard Torkar, Rik Huijzer, Martin Trapp and contributors. Last modified: December 15, 2022. Website built with Franklin.jl and the Julia programming language.