カリフォルニア大学バークレー校の学部別合格率
参考資料
Richard McElreath, Statistical Rethinking: A Bayesian Course with Examples in R and Stan, CRC Press/Taylor & Francis Group, 2016. Chapter13
データの概要
カリフォルニア大学バークレー校の入学試験のデータ。 学部別・性別の応募者数と合格者数が入っている。 学部数は6, 性別はmale/femaleの2カテゴリ, 合計12レコード。
- dept: 学部
- aplicant.gender: 性別
- admit: 合格者数
- reject: 不合格者数
- appllications: 応募者数
data(UCBadmit)
d <- UCBadmit
d$male <- ifelse(d$applicant.gender=="male", 1, 0)
d$dept_id <- coerce_index(d$dept)
head(d)分析の目的
学部ごとに性別によって合格率に差異があるかを明らかにする。
モデル式
学部\(i\)の合格者数は, 応募者数\(N_i\), 合格率\(p_i\)の二項分布に従うと考える。合格率\(p_i\)は学部による効果\(\alpha_{DEPT}\)と性別による効果\(\beta\)によって決まる。学部による効果\(\alpha_{DEPT}\)は平均\(\alpha\), 標準偏差\(\sigma\)の正規分布に従うとする。性別による効果は全ての学部で同じであるとする。
\[ \begin{align} A_i &\sim {\rm Binomial}(N_i, p_i)\\ p_i &\sim {\rm inv\_logit}(\alpha_{{\rm DEPT}[i]} + \beta M_i)\\ \alpha_{\rm DEPT} &\sim {\rm Normal}(\alpha, \sigma)\\ \alpha &\sim {\rm Normal}(0, 10)\\ \beta &\sim {\rm Normal}(0, 1)\\ \sigma &\sim{\rm HalfCauchy}(0,2) \end{align} \]
パラメータ推定実行
推定結果
fit## Inference for Stan model: AdmissionModel1.
## 4 chains, each with iter=4500; warmup=500; thin=1;
## post-warmup draws per chain=4000, total post-warmup draws=16000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## a -0.59 0.01 0.65 -1.90 -0.96 -0.59 -0.22
## a_dept[1] 0.68 0.00 0.10 0.48 0.61 0.67 0.74
## a_dept[2] 0.63 0.00 0.12 0.40 0.55 0.63 0.71
## a_dept[3] -0.58 0.00 0.07 -0.73 -0.63 -0.58 -0.53
## a_dept[4] -0.62 0.00 0.09 -0.78 -0.67 -0.62 -0.56
## a_dept[5] -1.06 0.00 0.10 -1.26 -1.13 -1.06 -0.99
## a_dept[6] -2.61 0.00 0.16 -2.92 -2.71 -2.60 -2.50
## b -0.10 0.00 0.08 -0.25 -0.15 -0.10 -0.04
## sigma_dept 1.48 0.01 0.57 0.79 1.09 1.35 1.72
## p[1] 0.64 0.00 0.02 0.61 0.63 0.64 0.65
## p[2] 0.66 0.00 0.02 0.62 0.65 0.66 0.68
## p[3] 0.63 0.00 0.02 0.59 0.62 0.63 0.64
## p[4] 0.65 0.00 0.03 0.60 0.63 0.65 0.67
## p[5] 0.34 0.00 0.02 0.30 0.32 0.34 0.35
## p[6] 0.36 0.00 0.02 0.33 0.35 0.36 0.37
## p[7] 0.33 0.00 0.02 0.29 0.32 0.33 0.34
## p[8] 0.35 0.00 0.02 0.31 0.34 0.35 0.36
## p[9] 0.24 0.00 0.02 0.20 0.23 0.24 0.25
## p[10] 0.26 0.00 0.02 0.22 0.25 0.26 0.27
## p[11] 0.06 0.00 0.01 0.05 0.06 0.06 0.07
## p[12] 0.07 0.00 0.01 0.05 0.06 0.07 0.08
## lp__ -2602.21 0.03 2.18 -2607.26 -2603.50 -2601.90 -2600.62
## 97.5% n_eff Rhat
## a 0.72 11005 1
## a_dept[1] 0.87 8940 1
## a_dept[2] 0.86 9474 1
## a_dept[3] -0.44 16000 1
## a_dept[4] -0.45 16000 1
## a_dept[5] -0.86 16000 1
## a_dept[6] -2.31 16000 1
## b 0.06 7291 1
## sigma_dept 2.89 9603 1
## p[1] 0.67 16000 1
## p[2] 0.70 8923 1
## p[3] 0.67 16000 1
## p[4] 0.70 9479 1
## p[5] 0.37 16000 1
## p[6] 0.39 16000 1
## p[7] 0.37 16000 1
## p[8] 0.39 16000 1
## p[9] 0.28 16000 1
## p[10] 0.30 16000 1
## p[11] 0.08 16000 1
## p[12] 0.09 16000 1
## lp__ -2598.89 6839 1
##
## Samples were drawn using NUTS(diag_e) at Mon Sep 17 01:21:01 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).MCMCの収束判定
(1){bayeseplot}のmcmc_trace関数
- MCMCサンプル系列のトレースプロットを表示。
- parsでプロットするパラメータを指定する。
- regex_parsで配列・ベクトル・マトリック型のパラメータ(今回のa_deptやp)をまとめて指定する。
- facet_argsにlistとしてnrowやncolなどのオプションを渡すとプロットの並べ方を指定できる。
posterior <- as.matrix(fit)
mcmc_trace(posterior,
pars = c("a", "b"),
regex_pars = c("a_dept"),
facet_args = list(nrow = 3))
パラメータの事後分布の可視化
(1){bayeseplot}のmcmc_areas関数
- パラメータごとの事後分布をプロットする。
- parsでプロットするパラメータを指定する。
- regex_parsで配列・ベクトル・マトリック型のパラメータ(今回のa_deptやp)をまとめて指定する。
- 何%の信頼区間にするかはprobで指定する。
- rhatを指定するとパラメータのRhat値に基づいてプロットに着色する。
posterior <- as.matrix(fit)
mcmc_areas(posterior,
pars=c("a", "b", "sigma_dept"),
regex_pars=c("a_dept"),
prob = 0.8,
rhat = rep(1, 9))
(2){bayeseplot}のmcmc_areas関数
- パラメータごとの事後分布をプロットする。
- 可視化している内容もオプションもmcmc_areasとほぼ同じ。
posterior <- as.matrix(fit)
mcmc_intervals(posterior,
pars=c("a", "b", "sigma_dept"),
regex_pars=c("a_dept"),
prob = 0.5,
rhat = rep(1, 9))
(3){rstan}のstan_dens関数
- パラメータごとの事後分布をプロットする。
- 引数にはstanfitオブジェクトがそのまま渡せる。
- parsでプロットするパラメータを指定する。配列・ベクトル・マトリックス型のパラメータはまとめて指定することができる(“a_dept”と指定するとa_deptの6要素が全てプロットされる)。
- n_rowとn_colでプロットの並べ方を指定できる。
stan_dens(fit,
pars=c("a", "a_dept", "b", "sigma_dept"),
nrow=5, ncol=2)
(4){bayeseplot}のmcmc_scatter関数
- 二つのパラメータの事後分布の相関関係を散布図で可視化する。
- parsでプロットするパラメータの組を指定する。
posterior <- as.matrix(fit)
mcmc_scatter(posterior, pars = c("a", "b"), alpha = 0.2)