Pharmacometrics解析のためのRStanへのTorstenライブラリの導入

Metrum Research Groupが開発しているTorstenライブラリは、Stan上でPharmacometrics解析を行うための拡張ライブラリです。 BaseのStanパッケージに比べ、以下の拡張がなされています。

NONMEMフォーマットのデータセットの読み込みとイベントスケジュール(薬剤投与、測定)の処理
標準薬物動態ライブラリの搭載 (NONMEMのPREDPPのようなもの)
- 線形1-/2-コンパートメントモデル+一次吸収の解析解
- 線形コンパートメントモデル用の汎用モデル (Padé approximantionによる解)
- 汎用ODEモデル (Runge KuttaもしくはBDF)

実用上では一番目のポイントがかなり重要で、わざわざStan用にデータセットを編集/Stanモデル内でイベントスケジュールの管理をする必要がなくなるため、Pharmacometrics解析にStanを使う心理的・労力的ハードルが低くなります。解析解などを用いている線形コンパートメントモデルはもちろん高速ですが、汎用モデルのODE solverもC++上で動いている(はず)ので、計算時間は比較的早いと思われます。

この記事では、TorstenライブラリをRStanに導入する方法を解説し、一例として2-コンパートメントモデルを使ってみます。

インストール
- 必要なもの
- 手順
使用例
参照したリンク先のまとめ

インストール

必要なもの

R (とRStudio)
devtoolsライブラリ(CRANからインストール可能)
RStanライブラリ

手順

最近までインストールが複雑でしたが、Metrumが専用ライブラリを作成してくれたので、RStanへのインストールが格段に楽になりました。

R上で以下のコマンドを入力するとインストールが行われます。RStanライブラリのヘッダーを置き換えてコンパイルをし直すようなので、インストールに比較的時間がかかります。

devtools::install_github('metrumresearchgroup/TorstenHeaders')
torstenHeaders::install_torsten()

使用例

こちらもMetrumによるGitHubリポジトリに、Torstenを使用した幾つかのモデルの例が載っています。今回はこのうちの2-コンパートメントモデルの例を使用してみます。
https://github.com/metrumresearchgroup/example-models
Torstenのマニュアルもこのリポジトリに含まれています(PDFへのリンク)。

はじめに、今回使用するライブラリを読み込んでおき、RStanの設定もしてしまいます。

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

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

データセット

解析に使用するデータセットをRに読み込みます。下記のRスクリプトをダウンロードした上で、read_rdump関数を使用して読み込みを行います。 https://github.com/metrumresearchgroup/example-models/blob/torsten-0.83/TwoCptModel/TwoCptModel.data.R

data <- rstan::read_rdump("TwoCptModel.data.R")

元のデータセットとしては、これらすべての要素がカラムとして格納されたデータフレームを想像してもらえるとイメージが掴めると思います。各行が一つのイベントに対応し、それぞれ薬物投与やサンプリングに対応しています。各データが何に対応しているかを簡単に解説すると、

amt: 薬物の投与量
addl: 薬物投与を計何回行うか
ii: 投与間隔
cmt: イベントが発生するコンパートメントのインデックス (1は吸収コンパートメント、2は中心コンパートメント)
cObs: 測定値
evid: 測定の場合は0、投与の場合は1に対応するインデックス
iObs: 全イベントのうち測定に対応する行のインデックス
nObs: 測定データの個数
nt: 全イベントの個数
rate: 定速静注の速度(今回は使用なし)
ss: 定常状態から計算の開始(今回は使用なし)
time: 各イベントが発生した時間

となっています。

実際にCSVファイルなどからデータを読み込む場合は、下記のRコードの57~75行目の様な処理を行います。基本的には、with関数を使って各カラムを独立したベクトルとしてリストに格納するだけです。この際、DVカラム(観測データが含まれる)に関しては、投与イベント等、観測データが欠損しているレコードを除く必要があります。インデックスiObsを使用することで、DVのうち観測データのみがcObsに格納されています。
https://github.com/metrumresearchgroup/example-models/blob/torsten-0.83/TwoCptModel/TwoCptModelSimulation.R

Stanモデル

以下に2-コンパートメントモデルを使用したモデルの例があげられています。1被験者のパラメータを推定するだけのシンプルなモデルです。
https://github.com/metrumresearchgroup/example-models/blob/torsten-0.83/TwoCptModel/TwoCptModel.stan

鍵となるのは65~69行目で、Torstenに含まれるPKModelTwoCpt関数を使用してODEの(解析)解を求めています。

  // PKModelTwoCpt takes in the NONMEM data, followed by the parameter
  // arrays and returns a matrix with the predicted amount in each 
  // compartment at each event.
  x = PKModelTwoCpt(time, amt, rate, ii, evid, cmt, addl, ss,
                   theta, biovar, tlag);

timeからssまではデータセットの中身をそのまま渡しています。 thetaがパラメータ(CL, Q, V1, V2, ka)を含む配列、biovarとtlagは追加のパラメータです。それぞれ、

CL: 中心コンパートメントからの消失クリアランス
Q: 中心コンパートメントと末梢コンパートメント間のクリアランス
V1: 中心コンパートメントの分布容積
V2: 末梢コンパートメントの分布容積
ka: 吸収速度定数
biovar: 各コンパートメントへのバイオアベイラビリティ
tlag: 各コンパートメントへの投与ラグタイム

を表します。
現在、PKModelOneCptとPKModelTwoCptの引数としてはこれらすべてが必要であり、例えばラグタイムを設定する必要がない場合でも0を渡す必要があるようです*1。

この関数を走らせると、渡したtimeにおける全コンパートメントの薬物量がmatrixとして返されます。

Rスクリプト

上記のStanコードを動かすためのRスクリプトは以下にあります。 https://raw.githubusercontent.com/metrumresearchgroup/example-models/torsten-0.83/R/TwoCpt.R

ただ、色々と不可欠で無いものが入っているので、最小限の要素のみ残したものを以下に記載します。ここでは、Stanモデルファイルがworking directlyに保存されていると仮定しています。また、データファイルは上記のコードを使って読み込み済みと仮定しています。

## create initial estimates
init <- function() {
  list(CL = exp(rnorm(1, log(10), 0.2)),
       Q = exp(rnorm(1, log(15), 0.2)),
       V1 = exp(rnorm(1, log(35), 0.2)),
       V2 =  exp(rnorm(1, log(105), 0.2)),
       ka = exp(rnorm(1, log(2), 0.2)),
       sigma = runif(1, 0.5, 2))
}

## run Stan
fit <- stan(file = "TwoCptModel.stan",
            data = data,
            init = init)

推定が終わったら、ggmcmcライブラリなどを用いて結果を評価します。

fit.param <- ggs(fit)
list.param <- c("CL", "Q", "V1", "V2", "ka")

fit.param.plot <- 
  fit.param %>% filter(Parameter %in% list.param)

ggs_density(fit.param.plot)
ggs_traceplot(fit.param.plot)
ggs_autocorrelation(fit.param.plot)
ggs_Rhat(fit.param.plot)

推定値に基づいたシミュレーション結果と測定値を比較します。

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

obs.pred.interval <- 
  mcmc.sample$cObsPred %>% 
  apply(MARGIN=2,quantile, prob=c(0.05, 0.5, 0.95)) %>% 
  t() %>% 
  tbl_df()

obs.pred.interval %>% 
  mutate(Time = data$time[data$iObs],
         cObs = data$cObs) %>% 
  ggplot(aes(Time, `50%`)) +
  geom_line() +
  geom_ribbon(aes(ymin=`5%`, ymax=`95%`),alpha=0.1) +
  geom_point(aes(y=cObs)) +
  scale_x_continuous(breaks=seq(0, 200, by=24)) +
  xlab("Time (hour)") +
  ylab("Concentration")

全体の傾向・ばらつき共によく再現されています。図では複数回投与後の薬物動態プロファイルに見えませんが、これはday 2-7の間で投与前濃度しか測定を行っていないため、それ以外の時点のシミュレーション結果が出力されていないことに依ります*2。 f:id:yoshidk6:20171030065007p:plain