Stan & POPPK (2): 階層モデルを用いた複数被験者の薬物動態解析
http://yoshidk6.hatenablog.com/entry/2017/05/24/134934 の続きです。前回使用したモデルを拡張して、複数の被験者からのデータを解析します。解析は以下の手順で行っていきます。
使用する仮想データ
作成された仮想データは前回と同じく、各被験者がTIME=0で薬物を経口投与された後、各時間ごとに血中濃度(CONC)を測定されたという状況を想定しています。 被験者は各10人ずつ4群に割り当てられており、それぞれ10, 20, 30, 40の投与量がT与えられています。
simple_pk_stan/sim_pk_20170521.csv at master · yoshidk6/simple_pk_stan · GitHub
図示してみると、投与量に応じて血中濃度が変動する一方で、血中濃度推移の形は変動せず、投与量に比例して濃度が上昇していることがわかります。
library(tidyverse) data.pk <- read_csv("https://raw.githubusercontent.com/yoshidk6/simple_pk_stan/master/data/sim_pk_20170521.csv") data.pk.plot <- mutate(data.pk, ID=factor(ID), DOSE_LEVEL=factor(DOSE_LEVEL)) data.pk.plot %>% ggplot(aes(TIME, CONC, group=ID, color=DOSE_LEVEL)) + geom_line() + geom_point() + facet_wrap(~DOSE_LEVEL) + scale_y_log10()
使用するデータは以下からダウンロード出来ます。 https://raw.githubusercontent.com/yoshidk6/simple_pk_stan/master/data/sim_pk_20170521.csv https://raw.githubusercontent.com/yoshidk6/simple_pk_stan/master/data/subj_dose_20170521.csv
全被験者が同じPKパラメータを持つと仮定
まずは、全ての被験者が同じPKパラメータを持つと仮定した、単純な解析を行ってみます。
Stanモデル
全被験者が同じパラメータを持つと仮定しているため、parameters
ブロックの設定や大まかな構成は前回と全く同じです。一方、今回のデータでは被験者ごとに異なる投与量が与えられているため、data
ブロックとモデル式にそれを反映する必要があります。主な変更点を以下で説明しています。
transformed parameters
今回のモデルでは、投与量DOSE
が各被験者ごとに与えられています。mu
の計算でベクトル演算を可能にするために、それぞれのデータ点に対応した投与量を設定し直す必要があります。このブロック内で、各被験者の投与量DOSE
と各データ点に対応する被験者のID
を組み合わせて、長さN
のベクトルDOSE_N
を生成しています。
muの計算式では、DOSE_N
とTIME
がベクトル化されているため、途中の掛け算がベクトルの要素同士の演算を行う.*
となっています。
GitHubリンク: https://github.com/yoshidk6/simple_pk_stan/blob/master/code/mod_01_noIIV.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)), s_Y = runif(1, 0.5, 2)) } data <- list(N = nrow(data.pk), N_ID = nrow(data.subj), ID = data.pk$ID, TIME = data.pk$TIME, DOSE = data.subj$DOSE, # 各被験者の投与量(長さN_ID) Y = data.pk$CONC) fit.stan <- stan(file = "mod_01_noIIV.stan", data = data, init = init)
結果の評価
推定されたパラメータを元に算出された平均値と90%予測区間を見てみると、全体の傾向はおおよそ再現できています。
一方、実測値と予測区間を個人ごとにプロットしてみると、個人内変動に比べて予測区間がかなり広くなっています。
PKパラメータは前回と近い推定値となりましたが、Yの標準偏差s_Y
は0.29と、前回の0.17に比べてかなり大きく推定されています。これは、各個人の生理学的な差(=PKパラメータの差)を無視して個人間変動・個人内変動の全てをs_Y
に負わせていることが原因と考えられます。
上図を生成するためのコードは以下の通りです。
mcmc.sample <- rstan::extract(fit.stan) # Plot for each subject (ID 1~9) y.pred.interval <- mcmc.sample$y_new %>% apply(MARGIN=2,quantile, prob=c(0.05, 0.5, 0.95)) %>% t() %>% tbl_df() bind_cols(data.pk, y.pred.interval) %>% filter(ID<=9) %>% mutate(ID=factor(ID)) %>% ggplot(aes(TIME, `50%`, color=ID)) + facet_wrap(~ID) + geom_line() + geom_ribbon(aes(ymin=`5%`, ymax=`95%`, fill=ID),alpha=0.1) + geom_point(data=data.pk %>% filter(ID<=9) %>% mutate(ID=factor(ID)), aes(TIME,CONC)) + ylab("CONC") # Plot for each dose level dat.fit.plot.each.dose <- ## Convert to tibble mcmc.sample$y_new %>% t() %>% tbl_df() %>% ## Re-groupe predicted values according to TIME and DOSE_LEVEL bind_cols(data.pk,.) %>% gather(key,prediction,-(ID:DOSE_LEVEL)) %>% group_by(TIME,DOSE_LEVEL) %>% ## Calc quantiles summarize(`5%` =quantile(prediction, probs=0.05), `50%`=quantile(prediction, probs=0.5), `95%`=quantile(prediction, probs=0.95)) %>% ungroup() %>% mutate(DOSE_LEVEL=factor(DOSE_LEVEL)) ggplot(dat.fit.plot.each.dose, aes(TIME, `50%`, color=DOSE_LEVEL)) + facet_wrap(~DOSE_LEVEL, scale="free") + geom_line() + geom_ribbon(aes(ymin=`5%`, ymax=`95%`),alpha=0.1) + geom_line( data=data.pk.plot, aes(TIME,CONC,group=ID), linetype="dotted") + geom_point(data=data.pk.plot, aes(TIME,CONC,group=ID)) + ylab("CONC")
階層モデルを導入した解析
各個人が異なるPKパラメータを持っていると考えた場合、各々のデータを用いて各個人のパラメータ推定を行うという方針もあります。実際、今回のデータセットならばその方針を取ることも可能です。しかし実世界では、各個人ごとに得られるデータの数はこれほど多くなく(例えば1被験者に対して1~2時点のデータのみなど)、推定計算を個別に安定して行うことが難しいというシチュエーションは珍しくありません。また、個人間のパラメータにはなんらかの類似性があると考えるのが自然であり、各個人ごとに推定を行うと他の被験者から得られる貴重な情報を無視してしまうことになります。
これらの問題点を解決するため、各個人ごとのパラメータは類似しており、ある特定の分布から確率的に生成されると見なすことで、得られるすべてのデータを用いて安定的かつ理にかなった推定を行えるようにするというのが階層モデルの考え方です(より一般的な単語として混合効果モデル・mixed effect modelと呼ばれることもあり、薬物動態の世界ではそのように呼ばれる事が殆どです)。
Stanモデル
大きな変更点として、被験者集団のパラメータ典型値であるCL
とVD
に加え、個人ごとのパラメータ値であるCLi
とVDi
が定義されています。これら個人ごとのパラメータは、平均log(CL)
とlog(VD)
、標準偏差s_CL
とs_VD
の対数正規分布に従って生成されると仮定しています*1。
transformed parameters
前節のDOSE_N
と同様に、CL_N
とVD_N
によって各被験者のパラメータが各測定時点に割り振られています。
また、パラメータもベクトルで表されるようになったため、多くの掛け算・割り算がベクトルの要素同士の演算を行う.*
と./
になっています。
model
このブロック内で、CLi
とVDi
が、平均がlog(CL)
とlog(VD)
・標準偏差がs_CL
とs_VD
の対数正規分布に従って生成されると設定しています。
GitHubリンク: https://github.com/yoshidk6/simple_pk_stan/blob/master/code/mod_02_withIIV.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)), s_CL= exp(rnorm(1, log(0.2), 0.5)), s_VD= exp(rnorm(1, log(0.2), 0.5)), s_Y = runif(1, 0.5, 2), CLi = rep(0.5,nrow(data.subj)), VDi = rep(5, 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_02_withIIV.stan", data = data, init = init)
結果の評価
推定されたパラメータを元に算出された平均値と90%予測区間を見てみると、全体の傾向は前節と変わりませんが、個人ごとのプロットでも観測された変動と予測区間がよく一致していることが分かります。Yの標準偏差s_Y
は0.15と小さくなり、正しく個人内変動を表していると考えられます。
図の作成に用いたコードは前節と同じです。
パラメータの相関
以上で見てきたように、階層モデルを用いて個人間差を正しく表現したパラメータ推定を行うことができると分かりました。一方で、推定された各個人のパラメータを詳しく見てみると、CLi
とVDi
が相関していることが見て取れます。
今回推定されたパラメータの分布を用いて新たな試験の結果予測を行おうとすると、この相関を考慮していないことが問題になります。例えば、上のグラフでCLi~0.2, VDi~8
というエリアには殆ど被験者がいませんが、単純に2つの対数正規分布を使用するとこのエリアにも仮想被験者が生成されてしまいます。
次のステップ
続いての記事では、上記の問題点を解決するため以下の評価を行っていきます。
- 多変量正規分布を用いたパラメータ間の相関まで含めた解析
- これら両方のパラメータ分布に関わっている内在的な因子(covariate)の評価
*1:薬物動態の世界では、多くの場合パラメータの個人間変動は対数正規分布に従うと考えられています。