Stan & POPPK (1): 一被験者の薬物動態解析

StanでPopulation pharmacokinetic解析 - yoshidk6’s blog の続きです。 最初のステップとして、一被験者のデータを用いた解析を行います。
(2017/5/27 に、データ・コード含めて更新を行いました。)
(2017/10/26 にコードの再更新を行いました。)

使用する仮想データ

作成された仮想データは、各被験者がTIME=0で薬物を経口投与された後、各時間ごとに血中濃度を測定(CONC)されたという状況を想定しています。 今回はこのうち、ID==1の被験者のデータを解析します。
simple_pk_stan/sim_pk_20170521.csv at master · yoshidk6/simple_pk_stan · GitHub

図示してみると、薬物が吸収されるに従って血中濃度CONCが上昇し、次第に体内から排出されるに従って減少していくことが分かります。

library(tidyverse)

data.pk <- read_csv("https://raw.githubusercontent.com/yoshidk6/simple_pk_stan/master/data/sim_pk_20170521.csv") 

data.pk %>%
  filter(ID==1) %>% 
  ggplot(aes(TIME, CONC)) +
  geom_line() +
  geom_point() 

Stanモデル

transformed parameters

このブロックで、1-コンパートメントモデルの解析解に基づいて濃度の算出を行っています(算出式に興味が有る方は一つ前の記事を御覧ください)。 薬物動態パラメータから計算された各時点での血中濃度の平均値がmuに代入されます。

model

Y(測定値)は標準偏差s_Yの対数正規分布に従うと設定しています。

各薬物動態パラメータについては、推定を安定させるため弱情報事前分布を設定しています。また、初期値もある程度plausibleな値を中心にランダムに生成するように設定しています*1

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

Rコード

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

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

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

data.pk   <- 
  read_csv("https://raw.githubusercontent.com/yoshidk6/simple_pk_stan/master/data/sim_pk_20170521.csv") 

## Prepare data for Stan
data.pk.id1 <- filter(data.pk, ID==1)

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.id1),
       TIME = data.pk.id1$TIME,
       DOSE = 10,  # ID==1の被験者への投与量は10と設定しています
       Y    = data.pk.id1$CONC)

## Run Stan
fit.stan <-
  stan(file = "mod_00_single_subj.stan", data = data)

結果の評価

安定した推定結果が得られ、推定値のmean (CL=0.33, VD=5.06, KA=0.93, s_y=0.17)も、もともと設定した値(CL=0.30, VD=4.89, KA=1.0, s_y=0.15)と非常に近い値となりました。推定されたパラメータを元に算出された平均値と90%予測区間も、実測値と比較的良く合っています。

上図を生成するためのコードは以下の通りです。

mcmc.sample <- rstan::extract(fit.stan)

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.id1,
          y.pred.interval) %>% 
  ggplot(aes(TIME, `50%`)) +
  geom_line() +
  geom_ribbon(aes(ymin=`5%`, ymax=`95%`), alpha=0.1) +
  geom_point(data=data.pk.id1, aes(TIME,CONC))

次の記事では、複数被験者のデータを用いたパラメータ推定を行います。

*1:今回用いているモデルには2つの解が存在してしまうことが知られているので(flip-flop)、ラベルスイッチングを何らかの方法で避ける必要がありました。一被験者の場合は値に上下関係をつけることも可能なのですが、後にあるように複数被験者に対して各個人のパラメータを推定すると制限をつけるのが難しくなってしまうため、今回は初期値をうまく設定することにしました。

StanでPopulation pharmacokinetic解析

簡単な薬物動態モデルを使ったPopulation pharmacokinetic (母集団薬物動態) 解析を、練習を兼ねてStanで実装してみました。

項目(予定)

1. 一被験者の薬物動態解析
2. 階層モデルによる複数被験者の薬物動態解析
3. 推定パラメータ間の相関と多変量正規分布
4. パラメータ個人差を説明する共変量の導入

モチベーション

• 薬物動態・臨床薬理界隈の人にStanを触ってもらうきっかけにする
• 薬物動態に馴染みはないが興味があるStan/Rユーザーへの紹介

モデル

仮想データの作成と実際のデータ解析共に、以下のシンプルな1-コンパートメントモデルを用いています。 消化管に投与された薬物が一次の吸収速度 kaに従って体内に移行し、体内では薬物が分布容積Vdに従って分布、最終的にクリアランスCLによって体外へと排出されるというモデルです(Web上にあまり良い図がなかったので、模式図を簡単に作成しました)。

このモデルからは簡単に解析解を導出することができます。興味の有る方は↓の記事を御覧ください（上の図と対応させるには、A₀を投与量/Vd、k₀をCL/Vdと読み替えてください）。
differential equations - Analytic solution to the one-compartment model - Mathematics Stack Exchange

解析コード・データ

仮想データ、その作成コード、及びに解析コードはGitHubで公開しています。
GitHub - yoshidk6/simple_pk_stan: Practice Stan with simple POPPK models and virtual data

ダウンロードされる場合はリポジトリごとダウンロードするのが一番楽かと思います。↓を参考にしてください。
GitHubで公開されているファイルをダウンロードするときの注意点（右クリックじゃできないよ） | 非IT企業に勤める中年サラリーマンのIT日記

仮想データはmrgsolveを用いて作成しました。
GitHub - metrumresearchgroup/mrgsolve: Simulate from ODE-based population PK/PD and QSP models in R

Acknowledgement

Stanの使用を始めるにあたり、 StanとRでベイズ統計モデリング / 松浦 健太郎 著 石田 基広 監修 市川 太祐 高橋 康介 高柳 慎一 福島 真太朗 編集 | 共立出版 を頭から読んで勉強しました。実践と背景知識のバランスがとれた素晴らしい本です。このブログの書式も作者の方のブログ StatModeling Memorandum に大いに影響されています。