TuringModels

Varying intercepts admission decisions

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

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]
df
12×7 DataFrame
 Row │ dept       gender     admit  reject  applications  male   dept_id
     │ InlineSt…  InlineSt…  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        6

Model

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     = 57.01 seconds
Compute duration  = 57.01 seconds
parameters        = Rho[1,2], a_bm_dept[2,5], sigma_dept[2], sigma_dept[1], a_bm_dept[1,2], a_bm_dept[2,6], a_bm_dept[1,3], a_bm_dept[2,2], Rho[2,2], a_bm_dept[1,4], a_bm_dept[2,3], bm, a_bm_dept[1,6], a_bm_dept[2,4], a_bm_dept[2,1], Rho[2,1], a_bm_dept[1,1], a, Rho[1,1], a_bm_dept[1,5]
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.3242    0.3421     0.0108    0.0097    951.5501    0.9990       16.6906
        Rho[1,2]   -0.3242    0.3421     0.0108    0.0097    951.5501    0.9990       16.6906
        Rho[2,2]    1.0000    0.0000     0.0000    0.0000   1076.6366    0.9990       18.8847
   sigma_dept[1]    1.6552    0.5986     0.0189    0.0153    963.9273    1.0000       16.9077
   sigma_dept[2]    0.5102    0.2673     0.0085    0.0126    457.2282    1.0059        8.0200
              bm   -0.1751    0.2347     0.0074    0.0066    969.2411    0.9996       17.0009
               a   -0.3658    0.5502     0.0174    0.0180    951.9519    0.9994       16.6977
  a_bm_dept[1,1]    1.3132    0.2482     0.0078    0.0091    905.9137    0.9991       15.8902
  a_bm_dept[2,1]   -0.8069    0.2643     0.0084    0.0090    921.2907    0.9990       16.1599
  a_bm_dept[1,2]    0.7656    0.3364     0.0106    0.0121    989.5525    0.9992       17.3572
  a_bm_dept[2,2]   -0.2391    0.3371     0.0107    0.0119    971.9879    0.9991       17.0491
  a_bm_dept[1,3]   -0.6477    0.0871     0.0028    0.0031    984.9351    1.0021       17.2762
  a_bm_dept[2,3]    0.0810    0.1422     0.0045    0.0048   1103.1011    0.9990       19.3489
  a_bm_dept[1,4]   -0.6204    0.1061     0.0034    0.0035    840.3995    0.9990       14.7410
  a_bm_dept[2,4]   -0.0865    0.1417     0.0045    0.0047   1055.3392    0.9999       18.5112
  a_bm_dept[1,5]   -1.1303    0.1132     0.0036    0.0038   1228.2274    1.0001       21.5437
  a_bm_dept[2,5]    0.1232    0.1851     0.0059    0.0067   1014.5512    1.0023       17.7957
  a_bm_dept[1,6]   -2.5920    0.1999     0.0063    0.0052   1103.5988    0.9996       19.3576
  a_bm_dept[2,6]   -0.1313    0.2666     0.0084    0.0088    947.4611    0.9999       16.6189

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.8599   -0.5974   -0.3570   -0.0889    0.4294
        Rho[1,2]   -0.8599   -0.5974   -0.3570   -0.0889    0.4294
        Rho[2,2]    1.0000    1.0000    1.0000    1.0000    1.0000
   sigma_dept[1]    0.8804    1.2407    1.5198    1.9164    3.2900
   sigma_dept[2]    0.1771    0.3352    0.4635    0.6105    1.2524
              bm   -0.6632   -0.3045   -0.1734   -0.0405    0.2931
               a   -1.5158   -0.7031   -0.3637   -0.0035    0.7241
  a_bm_dept[1,1]    0.8513    1.1399    1.3111    1.4784    1.8206
  a_bm_dept[2,1]   -1.3446   -0.9923   -0.8072   -0.6185   -0.3174
  a_bm_dept[1,2]    0.0846    0.5605    0.7703    0.9691    1.4337
  a_bm_dept[2,2]   -0.9312   -0.4516   -0.2300   -0.0230    0.4067
  a_bm_dept[1,3]   -0.8162   -0.7055   -0.6452   -0.5881   -0.4815
  a_bm_dept[2,3]   -0.1859   -0.0175    0.0809    0.1754    0.3727
  a_bm_dept[1,4]   -0.8185   -0.6945   -0.6199   -0.5486   -0.4220
  a_bm_dept[2,4]   -0.3506   -0.1875   -0.0954    0.0095    0.1860
  a_bm_dept[1,5]   -1.3466   -1.2023   -1.1283   -1.0530   -0.9091
  a_bm_dept[2,5]   -0.2238   -0.0021    0.1187    0.2470    0.5051
  a_bm_dept[1,6]   -3.0067   -2.7276   -2.5769   -2.4488   -2.2390
  a_bm_dept[2,6]   -0.6838   -0.3156   -0.1214    0.0407    0.3866

using StatsPlots

StatsPlots.plot(chns)

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
""";