CPS1985データを使い、wageを従属変数、genderとageと独立変数とした重回帰分析を行う。 RStudioでの実行。R.4.3.3+Stools4.3を使用。
data {
int<lower=0> N; // サンプルサイズ
int<lower=0, upper=1> X_1[N]; // 独立変数（性別、0: 男性、1: 女性）
vector[N] X_2; // 独立変数: 年齢
vector[N] Y; // 従属変数（賃金）
}
parameters {
real alpha; // 切片
real beta_gender; // 性別の影響
real beta_age; // 年齢の影響
real<lower=0> sigma; // 誤差項の標準偏差
}
model {
vector[N] gender_effect = to_vector(X_1); // 性別が1の場合は1、それ以外（つまり0）は0
Y ~ normal(alpha + beta_gender * gender_effect + beta_age * X_2, sigma);
}
library(rstan) library(AER) # データの準備 N <- nrow(CPS1985) Y <- CPS1985$wage X_1 <- ifelse(CPS1985$gender == "male", 0, 1) X_2 <- CPS1985$age # データをリストにまとめる stan_data <- list(N = N, Y = Y, X_1 = X_1, X_2 = X_2) # Stanモデルのフィット fit <- sampling(stan_model, data = stan_data, iter = 2000, chains = 4, seed=1234) # 結果の出力 print(fit)
Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha 6.92 0.02 0.71 5.50 6.44 6.92 7.38 8.33 2068 1
beta_gender -2.29 0.01 0.44 -3.16 -2.59 -2.30 -1.99 -1.41 2818 1
beta_age 0.09 0.00 0.02 0.05 0.07 0.09 0.10 0.12 2050 1
sigma 4.95 0.00 0.16 4.66 4.84 4.95 5.06 5.27 3265 1
lp__ -1118.87 0.04 1.48 -1122.47 -1119.61 -1118.52 -1117.76 -1117.06 1524 1
Samples were drawn using NUTS(diag_e) at Thu Apr 18 02:40:46 2024.
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)
bayesplot::mcmc_trace(as.array(fit), pars = c("alpha", "beta_gender", "beta_age", "sigma"))
bayesplot::mcmc_pairs(
as.array(fit),
pars = c("alpha", "beta_gender", "beta_age", "sigma"),
off_diag_args = list(size = 0.5) # オフダイアゴナルのプロットの点の大きさ
)
事後予測チェック
library(bayesplot)
posterior_samples <- extract(fit)
# 予測値の生成
num_draws <- length(posterior_samples$alpha) # 事後サンプルの数
predicted_Y <- matrix(NA, nrow = num_draws, ncol = N) # 初期化
for (i in 1:num_draws) {
predicted_Y[i, ] <- rnorm(N,
mean = posterior_samples$alpha[i] +
posterior_samples$beta_gender[i] * X_1 +
posterior_samples$beta_age[i] * X_2,
sd = posterior_samples$sigma[i])
}
# 事後予測チェック
ppc_check <- bayesplot::pp_check(object = Y, yrep = predicted_Y, fun="dens_overlay")
print(ppc_check)
実効サンプルサイズ
fit_summary <- monitor(fit) print(fit_summary)
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat valid
alpha 6.921184e+00 0.0155193212 0.70682223 5.502206e+00 6.440878e+00 6.924570e+00 7.381204e+00 8.3312505 2070 1.000772 1
beta_gender -2.291384e+00 0.0083537210 0.44279597 -3.161696e+00 -2.592320e+00 -2.298043e+00 -1.992746e+00 -1.4104326 2798 1.001946 1
beta_age 8.559054e-02 0.0003989017 0.01808944 4.954774e-02 7.318514e-02 8.540466e-02 9.828155e-02 0.1209101 2051 1.000292 1
sigma 4.951900e+00 0.0027256369 0.15593607 4.658480e+00 4.843548e+00 4.946649e+00 5.058095e+00 5.2701769 3265 1.000555 1
lp__ -1.118869e+03 0.0376995835 1.47783101 -1.122466e+03 -1.119609e+03 -1.118521e+03 -1.117764e+03 -1117.0557380 1526 1.000905 1
Q5 Q50 Q95 MCSE_Q2.5 MCSE_Q25 MCSE_Q50 MCSE_Q75 MCSE_Q97.5 MCSE_SD Bulk_ESS Tail_ESS
alpha 5.755831e+00 6.924570e+00 8.051954 0.0341940210 0.0179667943 0.0179815608 0.0174125228 0.041111513 0.0109835488 2076 1973
beta_gender -3.012165e+00 -2.298043e+00 -1.551494 0.0235019727 0.0141219316 0.0087052647 0.0147062847 0.027382883 0.0060181400 2834 2237
beta_age 5.623037e-02 8.540466e-02 0.114979 0.0009738498 0.0005416239 0.0005417372 0.0005132177 0.001034017 0.0002830609 2062 2106
sigma 4.704350e+00 4.946649e+00 5.220389 0.0075309437 0.0034408536 0.0036507362 0.0048615927 0.008442741 0.0019301027 3288 2363
lp__ -1.121690e+03 -1.118521e+03 -1117.174445 0.1265368331 0.0521949723 0.0342846729 0.0308008386 0.013573997 0.0266666435 1620 2162
感度分析
data {
int<lower=0> N;
int<lower=0, upper=1> X_1[N];
vector[N] X_2;
vector[N] Y;
}
parameters {
real alpha;
real beta_gender;
real beta_age;
real<lower=0> sigma;
}
model {
// 変更した事前分布
alpha ~ normal(7, 1); // 事後分布の95%確信区間から取りうる値を指定
beta_gender ~ normal(-2.3, 0.45);
beta_age ~ normal(0.085, 0.018);
sigma ~ inv_gamma(2, 1/(4.97^2)); // 逆ガンマ分布
// データに基づくモデル
vector[N] gender_effect = to_vector(X_1);
// 観測データの尤度
Y ~ normal(alpha + beta_gender * gender_effect + beta_age * X_2, sigma);
}
fit_sensitivity <- sampling(sensitivity_analysis, data = stan_data, iter = 2000, chains = 4, seed=1234) print(fit_sensitivity )
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha 6.95 0.01 0.46 6.06 6.64 6.94 7.26 7.87 2365 1 beta_gender -2.28 0.01 0.31 -2.88 -2.50 -2.29 -2.07 -1.68 3136 1 beta_age 0.08 0.00 0.01 0.06 0.08 0.08 0.09 0.11 2249 1 sigma 4.93 0.00 0.15 4.64 4.83 4.93 5.04 5.24 2997 1 lp__ -1123.62 0.03 1.35 -1126.96 -1124.32 -1123.35 -1122.60 -1121.85 1770 1
