Stan & POPPK (3): 多変量正規分布を用いたパラメータ間の相関まで含めた解析

Stan & POPPK (2): 階層モデルを用いた複数被験者の薬物動態解析 - yoshidk6’s blog の続きです。前回使用したモデルを拡張して、パラメータ間の相関も含めた解析を行います。

使用する仮想データ
解析

使用する仮想データ

作成された仮想データは前回と同じなので、説明を省略します。

解析

Stanモデル

モデル構造は前回とほぼ同じですが、多変量正規分布を扱うための変更を行っています。
Stan超初心者入門の114ページ目から詳しい説明が載っています。
簡単にまとめると、個人間変動を考慮している２つのパラメータに対して相関行列rhoと各々の標準偏差sigmaを定義し、そこから分散共分散行列Omegaを算出するようになっています。

GitHubリンク: https://github.com/yoshidk6/simple_pk_stan/blob/master/code/mod_03_multinorm.stan

Rコード

上記のモデルを動かすRのコードは以下の通りです。

library(tidyverse)
library(rstan)
library(ggmcmc)

rstan_options(auto_write=T)
options(mc.cores=parallel::detectCores())

data.pk   <- read_csv("../data/sim_pk_20170521.csv") 
data.subj <- read_csv("../data/subj_dose_20170521.csv")

init <- function(){
    list(KA = exp(rnorm(1, log(0.5), 0.2)),
         CL = exp(rnorm(1, log(0.5), 0.2)),
         VD = exp(rnorm(1, log(5), 0.2)),
         sigma = exp(rnorm(2, log(0.2), 0.5)),
         rho = diag(2),
         s_Y = runif(1, 0.5, 2),
         CLVDiLog=matrix(rep(log(c(0.5,5)), ea = nrow(data.subj)), 
                      nrow = nrow(data.subj)))
}

data <- 
  list(N    = nrow(data.pk),
       N_ID = nrow(data.subj),
       ID   = data.pk$ID,
       TIME = data.pk$TIME,
       DOSE = data.subj$DOSE,
       Y    = data.pk$CONC)

fit.stan <-
  stan(file = "mod_03_multinorm.stan", 
       data = data, init=init)