はじめに
具体的な問題としては，Stanの出力結果で配列化したパラメータをRで受け取った時の問題です。
配列化されたstanfitオブジェクトは次のように出力されます。
> fit
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -492.44 -492.03 11.36 11.31 -511.69 -474.75 1.00 760 1384
theta[1,1] 0.34 0.29 0.24 0.26 0.02 0.78 1.00 6872 2172
theta[2,1] 0.33 0.28 0.24 0.27 0.02 0.78 1.00 7548 2512
theta[3,1] 0.34 0.30 0.23 0.27 0.02 0.77 1.00 6163 2511
theta[4,1] 0.33 0.29 0.24 0.27 0.02 0.78 1.00 6515 2574
theta[5,1] 0.33 0.30 0.24 0.27 0.03 0.78 1.00 5928 2481
theta[6,1] 0.34 0.30 0.24 0.27 0.02 0.79 1.00 7291 1701
theta[7,1] 0.45 0.43 0.28 0.35 0.04 0.93 1.00 3659 2281
theta[8,1] 0.59 0.63 0.27 0.32 0.09 0.97 1.00 4358 2337
theta[9,1] 0.59 0.62 0.27 0.33 0.10 0.97 1.00 4639 2536
これをextractして,tibble型にしてもこんな感じです。
> df <- fit.stanfit %>% rstan::extract() %>% as.data.frame %>% as_tibble
> df
# A tibble: 4,000 × 1,013
theta.1.1 theta.2.1 theta.3.1 theta.4.1 theta.5.1 theta.6.1 theta.7.1 theta.8.1 theta.9.1 theta.10.1 theta.11.1 theta.12.1
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.120 0.0796 0.584 0.264 0.00837 0.592 0.484 0.279 0.530 0.814 0.969 0.510
2 0.293 0.309 0.0590 0.0541 0.0484 0.0178 0.183 0.163 0.292 0.979 0.949 0.633
3 0.406 0.200 0.0668 0.158 0.244 0.433 0.446 0.0483 0.262 0.179 0.954 0.146
4 0.133 0.583 0.216 0.480 0.400 0.0277 0.749 0.499 0.982 0.288 0.708 0.145
5 0.0257 0.242 0.384 0.0892 0.0582 0.261 0.904 0.636 0.868 0.965 0.624 0.974
6 0.382 0.575 0.362 0.585 0.0867 0.138 0.0604 0.0633 0.192 0.678 0.614 0.513
7 0.881 0.372 0.476 0.575 0.607 0.582 0.183 0.647 0.860 0.749 0.558 0.967
8 0.371 0.166 0.304 0.289 0.440 0.499 0.490 0.746 0.621 0.575 0.308 0.903
9 0.642 0.356 0.241 0.239 0.0686 0.336 0.128 0.738 0.590 0.984 0.813 0.817
これを縦長にして，なんとかここまできたとします。
> df %>% rowid_to_column("iter") %>% pivot_longer(-iter)
# A tibble: 4,052,000 × 3
iter name value
<int> <chr> <dbl>
1 1 theta.1.1 0.120
2 1 theta.2.1 0.0796
3 1 theta.3.1 0.584
4 1 theta.4.1 0.264
5 1 theta.5.1 0.00837
6 1 theta.6.1 0.592
7 1 theta.7.1 0.484
8 1 theta.8.1 0.279
9 1 theta.9.1 0.530
10 1 theta.10.1 0.814
あとは名前でgroup_byしてsummariseすればEAPとかは出るんですが・・・
> df %>% rowid_to_column("iter") %>% pivot_longer(-iter) %>%
+ group_by(name) %>% summarise(EAP = mean(value))
# A tibble: 1,013 × 2
name EAP
<chr> <dbl>
1 lp__ -492.
2 lp.1.1 -12.6
3 lp.1.2 -1.82
4 lp.10.1 -1.50
5 lp.10.2 -5.68
6 lp.100.1 -1.14
7 lp.100.2 -10.3
8 lp.101.1 -1.25
9 lp.101.2 -11.0
10 lp.102.1 -1.37
# … with 1,003 more rows
ここでふと，「あ，ここから配列[i,j]を配列名，i, jの変数名にしたいな，と。
それが今日のお題。
Q.変数の中にある数字を使ってグループ化したい
次のデータフレームsampleを例にします。
sample <- data.frame(
name = c("theta.1.1","gamma.1.2","theta.2.1","lambda.2.2"),
value = c(0.120,0.0796,0.584,0.264)
)
これの中身。
> sample
name value
1 theta.1.1 0.1200
2 gamma.1.2 0.0796
3 theta.2.1 0.5840
4 lambda.2.2 0.2640
目標はこちら。
name value variables val1 val2
1 theta.1.1 0.1200 theta 1 1
2 gamma.1.2 0.0796 gamma 1 2
3 theta.2.1 0.5840 theta 2 1
4 lambda.2.2 0.2640 lambda 2 2
さあ，どうしましょう。
解答編はこちら。