StanとRでベイズ統計モデリング11
UPDATE: 2023-12-29 21:58:15.311843
はじめに
このノートは「StanとRでベイズ統計モデリング」の内容を写経することで、ベイズ統計への理解を深めていくために作成している。
基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、「StanとRでベイズ統計モデリング」を読み進めるための自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。
今回は第9章「一歩進んだ文法」のチャプターを写経していく。
9.3.1 多変量正規分布
32人分の50m走と走り幅跳びの距離が記録されているデータを利用する。50m走と走り幅跳びの距離は相関することが想定でき、データを生成している確率分布は多変量正規分布が妥当だと考えられる。Stanでは、多変量正規分布を利用する際は、ベクトルを上手く利用する必要がある。
library(dplyr)
library(rstan)
library(ggplot2)
d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap09/input/data-mvn.txt')
ggplot(data = d, aes(x = Y1, y = Y2)) +
theme_bw(base_size = 18) +
geom_point(shape = 1, size = 3)
多変量正規分布を仮定したモデル式は下記の通り。これまでの表記と異なる部分は、\(Y\)がベクトル化されている点であり、Stanファイルで指定する時もvector[K] Y[N]と指定する。他にも、多変量正規分布ではパラメタもベクトルや行列であるため、Stanファイルもそれにあわせる必要がある。
\[ \begin{eqnarray} \overrightarrow{Y[n]} &\sim& MultiNormal(\overrightarrow{\mu}, \Sigma)\\ \end{eqnarray} \]
先程の点を注意して、データを渡す必要がある。
## $N
## [1] 32
##
## $D
## [1] 2
##
## $Y
## Y1 Y2
## 1 9.2 2.56
## 2 9.8 1.99
## 3 9.4 2.40
## 4 9.2 2.27
## 5 8.1 3.68
## 6 9.2 2.81
## 7 9.2 2.50
## 8 8.5 2.65
## 9 8.7 3.36
## 10 9.3 2.56
## 11 9.2 2.54
## 12 9.2 2.72
## 13 8.3 3.46
## 14 9.8 2.50
## 15 8.0 3.33
## 16 9.3 2.30
## 17 10.3 1.76
## 18 8.6 3.26
## 19 9.6 2.48
## 20 8.2 3.14
## 21 9.6 1.76
## 22 9.6 2.17
## 23 9.9 2.20
## 24 8.9 2.90
## 25 9.3 2.49
## 26 8.7 2.42
## 27 10.2 1.92
## 28 9.6 2.07
## 29 9.7 2.53
## 30 9.6 2.32
## 31 9.2 2.22
## 32 8.8 2.71Stanのファイルは下記の通り。vector[D] Y[N]は、\(Y\)は33個(=N)の「長さ2(=K)のベクトル」という意味である。cov_matrix[D] covは、共分散行列を扱える特殊な型で、D×Dサイズの行列が作られる。
data {
int N;
int D;
vector[D] Y[N];
}
parameters {
vector[D] mu;
cov_matrix[D] cov;
}
model {
for (n in 1:N){
Y[n] ~ multi_normal(mu, cov);
}
}\(Y\)は33個の長さ2のベクトルとは、下記のようなイメージである。下記はStanのサンプリング中にprint()関数でプリントデバッグした結果である。
[9.2,2.56]
[9.8,1.99]
[9.4,2.4]
[9.2,2.27]
[8.1,3.68]
[9.2,2.81]
[9.2,2.5]
[8.5,2.65]
[8.7,3.36]
[9.3,2.56]
[9.2,2.54]
[9.2,2.72]
[8.3,3.46]
[9.8,2.5]
[8,3.33]
[9.3,2.3]
[10.3,1.76]
[8.6,3.26]
[9.6,2.48]
[8.2,3.14]
[9.6,1.76]
[9.6,2.17]
[9.9,2.2]
[8.9,2.9]
[9.3,2.49]
[8.7,2.42]
[10.2,1.92]
[9.6,2.07]
[9.7,2.53]
[9.6,2.32]
[9.2,2.22]
[8.8,2.71]ここでは、stan_model()関数で最初にコンパイルしておいてから、
sampling()関数でサンプリングする。
推定結果は下記の通り。
## 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% 50% 97.5% n_eff Rhat
## mu[1] 9.2 0 0.1 9.0 9.2 9.4 2147 1
## mu[2] 2.6 0 0.1 2.4 2.6 2.8 2077 1
## cov[1,1] 0.4 0 0.1 0.3 0.4 0.7 1466 1
## cov[1,2] -0.3 0 0.1 -0.5 -0.3 -0.2 1348 1
## cov[2,1] -0.3 0 0.1 -0.5 -0.3 -0.2 1348 1
## cov[2,2] 0.3 0 0.1 0.2 0.3 0.5 1453 1
## lp__ 25.2 0 1.7 21.0 25.6 27.4 1222 1
##
## Samples were drawn using NUTS(diag_e) at Fri Dec 29 21:58:42 2023.
## 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).9.4.2 simplex型
simplex型はvector型の特別な場合の型で、要素が0-1の範囲で合計が1という条件を満たすことができる型。確率として解釈できるので、カテゴリカル分布や多項分布を利用する際に、あわせて利用されることが多い。
\[ \begin{eqnarray} Y[n] &\sim& Categorical(\overrightarrow{\theta}) \\ \end{eqnarray} \]
サイコロを200回振ったデータを利用する。このサイコロの各目の確率を推定したい。
K <- 6
d <- read.csv(file = 'https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap09/input/data-dice.txt')
d_count <- table(factor(d$Face, levels = 1:K))
data <- list(N = nrow(d), K = K, Y = d$Face)
data## $N
## [1] 200
##
## $K
## [1] 6
##
## $Y
## [1] 1 2 6 5 4 4 6 4 4 3 3 2 2 5 5 4 4 1 2 4 2 2 2 2 2 6 2 2 2 2 2 2 4 1 2 4 1
## [38] 2 2 1 4 4 2 3 1 4 2 6 2 4 4 4 6 3 4 3 1 2 2 2 2 2 1 4 2 5 3 1 2 4 4 3 2 1
## [75] 2 4 2 5 5 2 1 2 2 2 5 4 3 6 3 3 2 5 6 3 4 1 3 2 4 4 2 4 4 2 2 2 4 2 6 5 6
## [112] 1 4 4 2 6 2 2 6 3 3 2 2 2 4 4 4 2 6 2 1 4 2 4 2 5 3 4 4 2 4 5 2 2 4 2 2 1
## [149] 6 2 4 5 2 4 2 3 4 2 5 2 2 2 4 2 4 5 1 3 2 2 4 2 5 4 3 3 4 3 1 2 4 1 2 5 2
## [186] 6 4 2 1 2 2 1 4 1 4 2 5 5 4 5simplex[K] thetaと指定している部分に注意。このように指定することで、要素が0-1の範囲で合計が1という条件を満たす。
data {
int N;
int K;
int<lower=1, upper=K> Y[N];
}
parameters {
simplex[K] theta;
}
model {
for (n in 1:N)
Y[n] ~ categorical(theta);
}ここでは、stan_model()関数で最初にコンパイルしておいてから、
sampling()関数でサンプリングする。
推定結果は下記の通り。
## 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% 50% 97.5% n_eff Rhat
## theta[1] 0.1 0 0.0 0.1 0.1 0.2 5761 1
## theta[2] 0.4 0 0.0 0.3 0.4 0.4 5103 1
## theta[3] 0.1 0 0.0 0.1 0.1 0.1 5131 1
## theta[4] 0.3 0 0.0 0.2 0.3 0.3 5062 1
## theta[5] 0.1 0 0.0 0.1 0.1 0.1 4918 1
## theta[6] 0.1 0 0.0 0.0 0.1 0.1 5056 1
## lp__ -333.0 0 1.6 -336.9 -332.6 -330.9 2200 1
##
## Samples were drawn using NUTS(diag_e) at Fri Dec 29 21:59:03 2023.
## 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).似たようなケースとして、サイコロの目が出た数を集計したデータがあったとする。このデータは多項分布から生成されたと想定することができるので、下記のモデルを利用できる。ここで\(N\)は試行回数で、\(\overrightarrow{Y}\)は長さ\(K\)のベクトルで、各要素は\(N\)回中に何回の目が出たのかを記録している。
\[ \begin{eqnarray} \overrightarrow{Y} &\sim& Multinomial(N, \overrightarrow{\theta}) \\ \end{eqnarray} \]
# K <- 6
# d <- read.csv(file = 'https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap09/input/data-dice.txt')
# Y <- table(factor(d$Face, levels = 1:K))
# data <- list(K = K, Y = Y)
# model95 <- stan_model('note_ahirubayes11-95.stan')
# fit <- sampling(object = model95, data = data, seed = 1989)
# print(fit, prob = c(0.025, 0.5, 0.975), digits_summary = 1)
data {
int K;
int<lower=0> Y[K];
}
parameters {
simplex[K] theta;
}
model {
Y ~ multinomial(theta);
}9.5.2 欠損値
8.3で扱った分析データに欠損値が含まれる場合の対応を考える。欠損値はデータ分析にはよくあることで、何かしらの方法で処理する必要がある。このデータでは、PersonID=1,2,3,16の患者で欠損値が発生している。
d_wide <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap09/input/data-conc-2-NA-wide.txt')
apply(d_wide, 1, function(x) {sum(is.na(x))})## [1] 2 1 3 0 0 0 0 0 0 0 0 0 0 0 0 2欠損値が含まれたままの状態で、8.3で使用したN行T列の2次元配列を渡すとNAのインデックスのところでエラーになってしまう。そのため、wide型からlong型にテーブルを変形して、データを渡すのがよい。つまり、long型に変換した際に、欠損値のレコードを削除する。
\[ \begin{eqnarray} Y[n,t] &\sim& Normal(a[n]\{1 -\exp(-b[n] Time[t]) \}, \sigma_{Y}) \\ log(a[n]) &\sim& Normal(a_{全体平均}, \sigma_{a}) \\ log(b[n]) &\sim& Normal(b_{全体平均}, \sigma_{b}) \end{eqnarray} \]
long型に変更したことで、モデル式自体も変更する必要があるので注意。わかりにくいが、\(\mu[PersonID[i],TimeID[i]]\)となっており、PersonID,TimeIDのインデックスにあたる数値が正規分布の平均パラメタとして渡される。
\[ \begin{eqnarray} Y[i] &\sim& Normal(\mu[PersonID[i], TimeID[i]], \sigma_{Y}) \\ \mu[n,t] &=& a[n]\{1 -\exp(-b[n] Time[t]) \}\\ log(a[n]) &\sim& Normal(a_{全体平均}, \sigma_{a}) \\ log(b[n]) &\sim& Normal(b_{全体平均}, \sigma_{b}) \end{eqnarray} \]
モデル式はこちら。
data {
int I;
int N;
int T;
real Time[T];
int<lower=1, upper=N> PersonID[I];
int<lower=1, upper=T> TimeID[I];
vector[I] Y;
}
parameters {
real a0;
real b0;
vector[N] log_a;
vector[N] log_b;
real<lower=0> s_a;
real<lower=0> s_b;
real<lower=0> s_Y;
}
transformed parameters {
vector[N] a;
vector[N] b;
row_vector[T] mu[N];
a = exp(log_a);
b = exp(log_b);
for (n in 1:N){
for (t in 1:T){
mu[n,t] = a[n]*(1 - exp(-b[n]*Time[t]));
}
}
}
model {
log_a ~ normal(a0, s_a);
log_b ~ normal(b0, s_b);
for (i in 1:I){
Y[i] ~ normal(mu[PersonID[i], TimeID[i]], s_Y);
}
}モデルの挙動を確認しておく。muは観測ノイズを加える前の値であり、測定値が欠損かどうかに関わらず、すべての患者のすべての時点で存在すると考えられるため、欠損値の部分を推定されている。
// N=16
// T=6
transformed parameters {
for (n in 1:N){
for (t in 1:T){
mu[n,t] = a[n]*(1 - exp(-b[n]*Time[t]));
}
}
}
// row_vector[T] mu[N]の中身
| [
PersonID= 1 | [2.86439,5.12237,8.30545,11.5126,12.751,13.4852],
PersonID= 2 | [1.77624,3.24415,5.45976,8.00636,9.19418,10.1272],
PersonID= 3 | [6.6993,11.692,18.1857,23.7954,25.5258,26.275],
PersonID= 4 | [7.72858,13.2546,20.031,25.2666,26.635,27.1106],
PersonID= 5 | [2.7115,4.82007,7.73486,10.5634,11.5977,12.1649],
PersonID= 6 | [4.69957,7.88456,11.506,13.9333,14.4454,14.581],
PersonID= 7 | [2.20191,4.10481,7.17052,11.1702,13.4012,15.7269],
PersonID= 8 | [3.84458,6.54562,9.77645,12.1583,12.7385,12.9228],
PersonID= 9 | [11.5332,18.9006,26.6133,31.0449,31.7828,31.9296],
PersonID=10 | [6.22664,10.3172,14.7697,17.5206,18.0329,18.1495],
PersonID=11 | [6.65016,10.9271,15.447,18.0898,18.542,18.6348],
PersonID=12 | [5.93907,9.92335,14.3894,17.3039,17.8943,18.043],
PersonID=13 | [2.61715,4.46431,6.68817,8.34781,8.75965,8.89349],
PersonID=14 | [6.34325,10.9556,16.7481,21.4301,22.7389,23.2356],
PersonID=15 | [4.69129,8.15539,12.6021,16.3488,17.4626,17.9215],
PersonID=16 | [14.4284,22.0674,28.2533,30.4734,30.6478,30.6627]
| ]
// I=88
model {
for (i in 1:I){
Y[i] ~ normal(mu[PersonID[i], TimeID[i]], s_Y);
}
// Y[1] ~ normal(mu[PersonID[1], TimeID[1]], s_Y);
// Y[1] ~ normal(mu[1,1], s_Y);
// PersonIDが2のTime1はNAなので飛んで
// Y[2] ~ normal(mu[PersonID[2], TimeID[2]], s_Y);
// Y[2] ~ normal(mu[3,1], s_Y);
// Y[15] ~ normal(mu[PersonID[15], TimeID[15]], s_Y);
// Y[15] ~ normal(mu[16,1], s_Y);
// PersonIDが1のTime2はNAなので飛んで
// Y[16] ~ normal(mu[PersonID[16], TimeID[16]], s_Y);
// Y[16] ~ normal(mu[2,2], s_Y);Stanに渡すデータはこのようになる。
N <- nrow(d_wide)
Time <- c(1, 2, 4, 8, 12, 24)
colnames(d_wide) <- c('PersonID', 1:length(Time))
d <- d_wide %>%
tidyr::pivot_longer(cols = -PersonID, values_to = 'Y') %>%
mutate(TimeID = readr::parse_number(name)) %>%
select(-name) %>%
na.omit() %>%
arrange(TimeID)
data <- list(
I = nrow(d), N = N, T = length(Time), Time = Time,
PersonID = d$PersonID, TimeID = d$TimeID, Y = d$Y
)
data## $I
## [1] 88
##
## $N
## [1] 16
##
## $T
## [1] 6
##
## $Time
## [1] 1 2 4 8 12 24
##
## $PersonID
## [1] 1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 2 3 4 5 6 7 8 9 10 11
## [26] 12 13 14 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 4 5 6 7
## [51] 8 9 10 11 12 13 14 15 16 1 2 4 5 6 7 8 9 10 11 12 13 14 15 16 2
## [76] 4 5 6 7 8 9 10 11 12 13 14 15 16
##
## $TimeID
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3
## [39] 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6
## [77] 6 6 6 6 6 6 6 6 6 6 6 6
##
## $Y
## [1] 2.4 5.2 6.7 0.3 6.3 3.0 6.2 14.4 7.7 8.8 6.8 0.7 5.8 1.8 14.3
## [16] 3.9 9.4 12.6 4.7 3.8 4.2 6.8 17.0 10.0 11.3 12.6 4.6 9.4 7.9 7.5
## [31] 4.4 19.4 19.1 7.0 11.8 8.8 9.4 22.7 14.8 13.3 13.4 5.1 14.3 13.2 11.9
## [46] 7.7 23.4 10.2 9.2 15.4 11.3 29.8 15.3 19.2 14.9 9.9 23.1 15.6 31.9 12.5
## [61] 6.4 25.8 12.9 13.9 10.7 12.4 33.0 18.0 17.4 16.8 11.4 22.9 16.5 31.9 8.3
## [76] 26.1 14.8 18.2 16.2 14.7 32.2 18.7 21.8 16.0 11.2 25.2 15.9 30.6ここでは、stan_model()関数で最初にコンパイルしておいてから、
sampling()関数でサンプリングする。
推定結果は下記の通り。
## 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% 50% 97.5% n_eff Rhat
## a0 2.9 0.0 0.1 2.6 2.9 3.1 4043 1
## b0 -1.2 0.0 0.1 -1.5 -1.2 -0.9 3257 1
## log_a[1] 2.6 0.0 0.1 2.3 2.6 2.9 2877 1
## log_a[2] 2.1 0.0 0.1 1.8 2.1 2.4 3951 1
## log_a[3] 3.3 0.0 0.2 2.9 3.3 3.7 2734 1
## log_a[4] 3.3 0.0 0.0 3.2 3.3 3.3 4296 1
## log_a[5] 2.6 0.0 0.1 2.4 2.6 2.9 3082 1
## log_a[6] 2.8 0.0 0.1 2.6 2.7 3.0 2671 1
## log_a[7] 2.7 0.0 0.1 2.5 2.7 2.9 3556 1
## log_a[8] 2.6 0.0 0.1 2.4 2.6 2.8 4566 1
## log_a[9] 3.5 0.0 0.0 3.4 3.5 3.5 4736 1
## log_a[10] 2.9 0.0 0.1 2.8 2.9 3.0 3759 1
## log_a[11] 3.0 0.0 0.1 2.9 3.0 3.1 4019 1
## log_a[12] 2.8 0.0 0.1 2.7 2.8 2.9 3533 1
## log_a[13] 2.4 0.0 0.1 2.2 2.4 2.7 3427 1
## log_a[14] 3.2 0.0 0.1 3.1 3.2 3.3 3749 1
## log_a[15] 2.8 0.0 0.1 2.7 2.8 2.9 4944 1
## log_a[16] 3.5 0.0 0.0 3.4 3.5 3.5 5652 1
## log_b[1] -1.4 0.0 0.3 -2.1 -1.4 -0.8 2797 1
## log_b[2] -1.3 0.0 0.4 -2.1 -1.3 -0.6 3320 1
## log_b[3] -1.3 0.0 0.3 -1.9 -1.3 -0.7 2825 1
## log_b[4] -1.1 0.0 0.1 -1.4 -1.1 -0.8 4392 1
## log_b[5] -1.7 0.0 0.3 -2.2 -1.6 -1.1 2427 1
## log_b[6] -1.5 0.0 0.3 -2.1 -1.5 -0.9 2291 1
## log_b[7] -1.4 0.0 0.2 -1.9 -1.4 -0.9 3004 1
## log_b[8] -1.1 0.0 0.3 -1.6 -1.1 -0.5 4281 1
## log_b[9] -0.9 0.0 0.1 -1.2 -0.9 -0.7 4760 1
## log_b[10] -0.9 0.0 0.2 -1.3 -0.9 -0.5 4093 1
## log_b[11] -1.0 0.0 0.2 -1.4 -1.0 -0.5 4274 1
## log_b[12] -0.7 0.0 0.2 -1.1 -0.7 -0.2 3163 1
## log_b[13] -1.5 0.0 0.3 -2.1 -1.5 -0.9 2669 1
## log_b[14] -1.4 0.0 0.2 -1.7 -1.4 -1.1 3722 1
## log_b[15] -1.2 0.0 0.2 -1.6 -1.2 -0.8 4275 1
## log_b[16] -0.6 0.0 0.2 -1.0 -0.6 -0.3 2854 1
## s_a 0.4 0.0 0.1 0.3 0.4 0.6 3899 1
## s_b 0.4 0.0 0.1 0.2 0.4 0.7 1351 1
## s_Y 1.7 0.0 0.2 1.5 1.7 2.1 2396 1
## a[1] 13.6 0.0 2.1 10.2 13.3 18.5 2878 1
## a[2] 8.4 0.0 1.2 6.3 8.3 10.9 3854 1
## a[3] 26.9 0.1 6.0 18.1 25.9 41.5 2454 1
## a[4] 25.9 0.0 1.2 23.7 25.9 28.3 4260 1
## a[5] 14.2 0.0 1.6 11.2 14.0 17.7 3057 1
## a[6] 15.7 0.0 1.6 13.0 15.6 19.2 2638 1
## a[7] 14.9 0.0 1.4 12.4 14.8 17.8 3539 1
## a[8] 13.5 0.0 1.2 11.3 13.4 15.9 4534 1
## a[9] 31.9 0.0 1.1 29.7 31.9 34.0 4708 1
## a[10] 17.9 0.0 1.1 15.8 17.9 20.1 3741 1
## a[11] 19.6 0.0 1.1 17.4 19.6 22.1 4064 1
## a[12] 16.3 0.0 1.0 14.3 16.3 18.3 3531 1
## a[13] 11.6 0.0 1.4 9.1 11.6 14.7 3312 1
## a[14] 24.8 0.0 1.4 22.2 24.8 27.7 3778 1
## a[15] 16.7 0.0 1.1 14.6 16.7 19.0 4916 1
## a[16] 31.6 0.0 1.0 29.6 31.6 33.7 5652 1
## b[1] 0.3 0.0 0.1 0.1 0.2 0.5 3056 1
## b[2] 0.3 0.0 0.1 0.1 0.3 0.5 3477 1
## b[3] 0.3 0.0 0.1 0.1 0.3 0.5 3138 1
## b[4] 0.3 0.0 0.0 0.2 0.3 0.4 4370 1
## b[5] 0.2 0.0 0.1 0.1 0.2 0.3 2347 1
## b[6] 0.2 0.0 0.1 0.1 0.2 0.4 2402 1
## b[7] 0.2 0.0 0.1 0.1 0.2 0.4 3078 1
## b[8] 0.4 0.0 0.1 0.2 0.3 0.6 4163 1
## b[9] 0.4 0.0 0.1 0.3 0.4 0.5 4814 1
## b[10] 0.4 0.0 0.1 0.3 0.4 0.6 4061 1
## b[11] 0.4 0.0 0.1 0.3 0.4 0.6 4005 1
## b[12] 0.5 0.0 0.1 0.3 0.5 0.8 3297 1
## b[13] 0.2 0.0 0.1 0.1 0.2 0.4 2893 1
## b[14] 0.3 0.0 0.0 0.2 0.2 0.3 3436 1
## b[15] 0.3 0.0 0.1 0.2 0.3 0.5 4228 1
## b[16] 0.5 0.0 0.1 0.4 0.5 0.7 3148 1
## mu[1,1] 2.9 0.0 0.6 1.9 2.9 4.3 3556 1
## mu[1,2] 5.2 0.0 0.9 3.6 5.1 7.1 3621 1
## mu[1,3] 8.3 0.0 1.0 6.3 8.2 10.3 3859 1
## mu[1,4] 11.3 0.0 1.0 9.2 11.3 13.3 4426 1
## mu[1,5] 12.5 0.0 1.3 10.0 12.5 15.0 3831 1
## mu[1,6] 13.4 0.0 1.9 10.2 13.3 17.6 3022 1
## mu[2,1] 2.0 0.0 0.5 1.1 1.9 3.1 3798 1
## mu[2,2] 3.4 0.0 0.8 2.1 3.4 5.0 3865 1
## mu[2,3] 5.3 0.0 0.9 3.7 5.3 7.1 4086 1
## mu[2,4] 7.1 0.0 0.8 5.5 7.1 8.8 4938 1
## mu[2,5] 7.8 0.0 0.9 6.1 7.8 9.6 5084 1
## mu[2,6] 8.3 0.0 1.1 6.3 8.2 10.5 4239 1
## mu[3,1] 6.3 0.0 0.7 5.0 6.3 7.8 3899 1
## mu[3,2] 11.0 0.0 1.0 9.1 11.0 13.1 3966 1
## mu[3,3] 17.3 0.0 1.5 14.2 17.3 20.3 3713 1
## mu[3,4] 23.1 0.1 3.1 17.5 23.0 29.4 3095 1
## mu[3,5] 25.3 0.1 4.3 18.0 24.9 34.7 2811 1
## mu[3,6] 26.7 0.1 5.7 18.1 25.8 40.3 2536 1
## mu[4,1] 7.2 0.0 0.7 5.9 7.1 8.7 5264 1
## mu[4,2] 12.3 0.0 1.0 10.4 12.3 14.4 5561 1
## mu[4,3] 18.7 0.0 1.0 16.7 18.7 20.7 6446 1
## mu[4,4] 23.9 0.0 0.8 22.2 23.9 25.5 6404 1
## mu[4,5] 25.3 0.0 1.0 23.4 25.3 27.2 4846 1
## mu[4,6] 25.9 0.0 1.2 23.7 25.9 28.2 4290 1
## mu[5,1] 2.5 0.0 0.5 1.7 2.4 3.5 2658 1
## mu[5,2] 4.5 0.0 0.7 3.2 4.5 6.1 2764 1
## mu[5,3] 7.5 0.0 0.9 5.7 7.5 9.5 3097 1
## mu[5,4] 10.9 0.0 0.9 9.1 10.9 12.7 4432 1
## mu[5,5] 12.5 0.0 0.9 10.6 12.5 14.3 4480 1
## mu[5,6] 13.9 0.0 1.4 11.2 13.9 16.7 3343 1
## mu[6,1] 3.3 0.0 0.6 2.2 3.2 4.6 2728 1
## mu[6,2] 5.8 0.0 0.9 4.1 5.7 7.7 2832 1
## mu[6,3] 9.3 0.0 1.1 7.3 9.3 11.5 3134 1
## mu[6,4] 13.0 0.0 0.9 11.2 13.0 14.7 4095 1
## mu[6,5] 14.5 0.0 1.0 12.6 14.5 16.4 4119 1
## mu[6,6] 15.6 0.0 1.4 13.0 15.5 18.4 2826 1
## mu[7,1] 3.2 0.0 0.6 2.2 3.2 4.4 3621 1
## mu[7,2] 5.7 0.0 0.8 4.1 5.6 7.4 3744 1
## mu[7,3] 9.1 0.0 1.0 7.1 9.1 11.0 4156 1
## mu[7,4] 12.5 0.0 0.9 10.7 12.5 14.3 5241 1
## mu[7,5] 13.9 0.0 1.0 12.0 13.8 15.7 4765 1
## mu[7,6] 14.8 0.0 1.3 12.4 14.7 17.4 3687 1
## mu[8,1] 3.9 0.0 0.8 2.6 3.8 5.7 4878 1
## mu[8,2] 6.6 0.0 1.0 4.7 6.6 8.9 5120 1
## mu[8,3] 9.9 0.0 1.0 7.8 9.9 11.9 5615 1
## mu[8,4] 12.4 0.0 0.9 10.6 12.4 14.1 5650 1
## mu[8,5] 13.1 0.0 1.0 11.2 13.1 15.1 5007 1
## mu[8,6] 13.4 0.0 1.1 11.3 13.4 15.9 4571 1
## mu[9,1] 10.3 0.0 0.9 8.7 10.3 12.2 5441 1
## mu[9,2] 17.3 0.0 1.2 15.1 17.2 19.6 5688 1
## mu[9,3] 25.1 0.0 1.0 23.1 25.1 27.1 6577 1
## mu[9,4] 30.4 0.0 0.8 28.7 30.4 32.0 6018 1
## mu[9,5] 31.5 0.0 1.0 29.5 31.5 33.4 4996 1
## mu[9,6] 31.9 0.0 1.1 29.7 31.9 34.0 4714 1
## mu[10,1] 6.1 0.0 0.9 4.5 6.0 8.0 4980 1
## mu[10,2] 10.0 0.0 1.1 7.9 10.0 12.2 5325 1
## mu[10,3] 14.3 0.0 1.0 12.2 14.3 16.2 6088 1
## mu[10,4] 17.1 0.0 0.9 15.4 17.1 18.8 4975 1
## mu[10,5] 17.7 0.0 1.0 15.8 17.7 19.7 4021 1
## mu[10,6] 17.9 0.0 1.1 15.8 17.9 20.1 3750 1
## mu[11,1] 6.3 0.0 0.9 4.8 6.2 8.2 4942 1
## mu[11,2] 10.5 0.0 1.1 8.4 10.5 12.8 5287 1
## mu[11,3] 15.3 0.0 1.0 13.2 15.3 17.3 6300 1
## mu[11,4] 18.6 0.0 0.9 16.9 18.6 20.3 4784 1
## mu[11,5] 19.4 0.0 1.0 17.4 19.4 21.4 4093 1
## mu[11,6] 19.6 0.0 1.1 17.4 19.6 22.0 4056 1
## mu[12,1] 6.5 0.0 1.0 4.8 6.5 8.7 3932 1
## mu[12,2] 10.4 0.0 1.2 8.1 10.3 12.6 4231 1
## mu[12,3] 14.0 0.0 0.9 12.2 14.0 15.7 5355 1
## mu[12,4] 15.9 0.0 0.9 14.2 15.9 17.6 4511 1
## mu[12,5] 16.2 0.0 1.0 14.3 16.2 18.1 3765 1
## mu[12,6] 16.3 0.0 1.0 14.3 16.3 18.3 3537 1
## mu[13,1] 2.3 0.0 0.5 1.5 2.3 3.4 3380 1
## mu[13,2] 4.2 0.0 0.7 2.8 4.1 5.7 3507 1
## mu[13,3] 6.8 0.0 0.9 5.0 6.8 8.6 3922 1
## mu[13,4] 9.5 0.0 0.9 7.7 9.5 11.2 4992 1
## mu[13,5] 10.6 0.0 0.9 8.7 10.6 12.4 4815 1
## mu[13,6] 11.5 0.0 1.3 9.1 11.5 14.2 3630 1
## mu[14,1] 5.5 0.0 0.6 4.4 5.5 6.8 3744 1
## mu[14,2] 9.7 0.0 0.9 8.1 9.7 11.7 3845 1
## mu[14,3] 15.6 0.0 1.0 13.6 15.6 17.7 4069 1
## mu[14,4] 21.3 0.0 0.9 19.5 21.3 23.1 4233 1
## mu[14,5] 23.5 0.0 1.0 21.5 23.4 25.4 3889 1
## mu[14,6] 24.7 0.0 1.3 22.2 24.7 27.5 3761 1
## mu[15,1] 4.5 0.0 0.6 3.3 4.4 5.8 5020 1
## mu[15,2] 7.7 0.0 0.9 6.0 7.7 9.5 5231 1
## mu[15,3] 11.8 0.0 1.0 9.8 11.8 13.6 5750 1
## mu[15,4] 15.2 0.0 0.9 13.5 15.2 16.9 5997 1
## mu[15,5] 16.2 0.0 1.0 14.4 16.2 18.1 5388 1
## mu[15,6] 16.7 0.0 1.1 14.6 16.7 19.0 4967 1
## mu[16,1] 13.1 0.0 1.7 9.9 13.2 16.5 2934 1
## mu[16,2] 20.7 0.0 1.9 16.7 20.8 24.4 2826 1
## mu[16,3] 27.7 0.0 1.5 24.5 27.9 30.2 2935 1
## mu[16,4] 31.1 0.0 1.0 29.2 31.1 33.0 5933 1
## mu[16,5] 31.6 0.0 1.0 29.6 31.6 33.6 5951 1
## mu[16,6] 31.6 0.0 1.0 29.6 31.6 33.7 5657 1
## lp__ -79.7 0.1 5.3 -91.0 -79.3 -70.5 1288 1
##
## Samples were drawn using NUTS(diag_e) at Fri Dec 29 21:59:27 2023.
## 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).可視化用のデータを作成する。
ms <- rstan::extract(fit)
d$Time <- Time[d$TimeID]
d_est <- data.frame()
for (n in 1:N) {
qua <- apply(ms$mu[,n,], 2, quantile, prob = c(0.025, 0.5, 0.975))
d_est <- rbind(d_est, data.frame(PersonID = n, Time = Time, t(qua), check.names = FALSE))
}このデータを列方向にパーセンタイルを計算する。
## [1] 4000 6結果として、時点(TimeID)ごとのパーセンタイルが計算されているので、これをrbind()関数で縦に積み重ねる。
## 2.5% 50% 97.5%
## [1,] 1.926651 2.884747 4.286138
## [2,] 3.607455 5.135751 7.107914
## [3,] 6.289151 8.240309 10.332984
## [4,] 9.215038 11.322596 13.282221
## [5,] 9.991287 12.523979 15.033055
## [6,] 10.219319 13.256807 17.600268可視化した結果を見ると、欠損値が発生している患者では信用区間が広くなっているが、エラーなく推定できている。
ggplot(data = d_est, aes(x = Time, y = `50%`)) +
theme_bw(base_size = 18) +
facet_wrap(~PersonID) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = 'black', alpha = 1/5) +
geom_line(linewidth = 0.5) +
geom_point(data = d, aes(x = Time, y = Y), size = 3) +
labs(x = 'Time (hour)', y = 'Y') +
scale_x_continuous(breaks = Time, limit = c(0,24))