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()

f:id:yoshidk6:20170529072239p:plain

使用するデータは以下からダウンロード出来ます。 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_NTIMEがベクトル化されているため、途中の掛け算がベクトルの要素同士の演算を行う.*となっています。

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%予測区間を見てみると、全体の傾向はおおよそ再現できています。 f:id:yoshidk6:20170529074032p:plain

一方、実測値と予測区間を個人ごとにプロットしてみると、個人内変動に比べて予測区間がかなり広くなっています。 f:id:yoshidk6:20170529074014p:plain

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モデル

大きな変更点として、被験者集団のパラメータ典型値であるCLVDに加え、個人ごとのパラメータ値であるCLiVDiが定義されています。これら個人ごとのパラメータは、平均log(CL)log(VD)、標準偏差s_CLs_VDの対数正規分布に従って生成されると仮定しています*1

transformed parameters

前節のDOSE_Nと同様に、CL_NVD_Nによって各被験者のパラメータが各測定時点に割り振られています。 また、パラメータもベクトルで表されるようになったため、多くの掛け算・割り算がベクトルの要素同士の演算を行う.*./になっています。

model

このブロック内で、CLiVDiが、平均がlog(CL)log(VD)・標準偏差がs_CLs_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と小さくなり、正しく個人内変動を表していると考えられます。

f:id:yoshidk6:20170529080239p:plain f:id:yoshidk6:20170529080213p:plain

図の作成に用いたコードは前節と同じです。

パラメータの相関

以上で見てきたように、階層モデルを用いて個人間差を正しく表現したパラメータ推定を行うことができると分かりました。一方で、推定された各個人のパラメータを詳しく見てみると、CLiVDiが相関していることが見て取れます。 f:id:yoshidk6:20171027123546p:plain

今回推定されたパラメータの分布を用いて新たな試験の結果予測を行おうとすると、この相関を考慮していないことが問題になります。例えば、上のグラフでCLi~0.2, VDi~8というエリアには殆ど被験者がいませんが、単純に2つの対数正規分布を使用するとこのエリアにも仮想被験者が生成されてしまいます。

次のステップ

続いての記事では、上記の問題点を解決するため以下の評価を行っていきます。

  • 多変量正規分布を用いたパラメータ間の相関まで含めた解析
  • これら両方のパラメータ分布に関わっている内在的な因子(covariate)の評価

*1:薬物動態の世界では、多くの場合パラメータの個人間変動は対数正規分布に従うと考えられています。