2017-06-02

Metrum Research GroupによるSTANの入門コース

Metrum Research Groupから昨日、STANの入門コースの教材がアップロードされました*1。Metrumは、STANをPharmacometricsに適用するため積極的に開発に携わっています。今回公開されたものは、昨年度のPAGE/ACoP meetingで開催された Getting Started with Bayesian PKPD Modeling Using Stan: Practical use of Stan & R for PKPD applications というワークショップが元になっており、基礎的なSTANの使用法に重点が置かれています。コードや解説PDFに加え、実際にコードを動かすチュートリアルもYouTubeにアップロードされています。英語に抵抗の無い方はこのブログより百倍ためになると思うので、ご覧になることをおすすめします。

metrumrg.com

以下に簡単にコース内容をまとめてみました。番号はアップロードされている動画の番号に(おおよそ)対応しています。特に既にStanに触ったことのある方は、興味がある場所から見てみて下さい。

ベイズ統計・統計モデリングの基本
Stanの紹介・インストール・linear regressionを用いた簡単なデモ
単純なPK-PDモデル(Emaxモデル)を用いた非線形回帰*2と階層モデル（試験間差）への拡張
個人間差を考慮したPK-PDモデル(Emaxモデル)
User-defined functionとそれを用いたPopulation PK解析
5.の続きとCensored dataに対する尤度計算法(解説のみ)

Population PK解析で用いているコンパートメントモデルには、このチュートリアルでも解析解を用いており、数値的にODE(常微分方程式)は解いていません。とはいえ、僕の書いた様な簡易型のものではなく、しっかりとNONMEM形式のデータに対応することを意識して書かれています。ODEの数値解法*3についてカバーされていなかったのは残念ですが、今年のPAGE/ACoPでMetrumが開催するワークショップでカバーするようなので楽しみです。

*1:この分野にいる人はご存知だと思いますが、Metrumはこれに限らず幅広いトピックに関する教材をCreative Commonsライセンスの元で公開しています。

*2:いきなりODEを解く必要があるPOPPKから入るのではなく、Emaxモデルから入るのは上手いですね。

*3:Torstenという名前のサブモジュールが開発中のようです

2017-05-29

Stan & POPPK (2): 階層モデルを用いた複数被験者の薬物動態解析

http://yoshidk6.hatenablog.com/entry/2017/05/24/134934 の続きです。前回使用したモデルを拡張して、複数の被験者からのデータを解析します。解析は以下の手順で行っていきます。

使用する仮想データ
全被験者が同じPKパラメータを持つと仮定
- Stanモデル
  - transformed parameters
- Rコード
- 結果の評価
階層モデルを導入した解析
- Stanモデル
  - transformed parameters
  - model
- Rコード
- 結果の評価
- パラメータの相関
次のステップ

使用する仮想データ

作成された仮想データは前回と同じく、各被験者が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

2017-05-28

Google Cloud PlatformのRにRStanをインストールする

RStudio ServerをGoogle Computing Engineで動かす - yoshidk6’s blog に引き続き、VMインスタンスにRStanをインストールします。
以下の記事を参考にしてインストールを進めていきます。
Installing RStan on Mac or Linux · stan-dev/rstan Wiki · GitHub

Toolchain

build-essential, g++, libssl-dev をインストールします。

sudo apt-get update
sudo apt-get install build-essential g++ libssl-dev

Configuration

dpkg -s g++でg++のバージョンが4.9以降であることを確認します。 RStudio serverにログインし、以下のコードを実行してMakevarsファイルを作成します。

dotR <- file.path(Sys.getenv("HOME"), ".R")
if (!file.exists(dotR)) dir.create(dotR)
M <- file.path(dotR, "Makevars")
if (!file.exists(M)) file.create(M)
cat("\nCXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function", 
    file = M, sep = "\n", append = TRUE)

# Run only if g++ version is 4.9 or higher
cat("\nCXXFLAGS+=-flto -ffat-lto-objects  -Wno-unused-local-typedefs", 
    file = M, sep = "\n", append = TRUE)

Installing RStan

引き続きRStudio server上で以下のコードを実行し、RStanをインストールします。

install.packages("rstan", repos = "https://cloud.r-project.org/", dependencies=TRUE)

インストールに成功したら、Session -> Restart R でRを再起動します。
以下のコードを実行し、10が返されることを確認すればRStanのインストールは完了です。

fx <- inline::cxxfunction( signature(x = "integer", y = "numeric" ) , '
    return ScalarReal( INTEGER(x)[0] * REAL(y)[0] ) ;
' )
fx( 2L, 5 ) # should be 10

2017-05-27

RStudio ServerをGoogle Computing Engineで動かす

モチベーション

家で使っているコンピュータがMacbook Air 2011モデル(1.6GHz Core i5, メモリ2GB)と非常に貧弱で、Stanを動かすのに支障が出ている一方で、特段買いたいパソコンも見当たらないので、クラウドのサーバーを試してみることにしました。AWSについては公式含めかなり多くの解説記事がありますが、今回はGoogleのCloudサービスであるGoogle Computing Engineを試してみます。

モチベーション
Virtual Machineをセットアップ
R/RStudioのインストール
tidyverseのインストール
ファイル転送
最後に
References

Virtual Machineをセットアップ

Quickstart Using a Linux VM | Compute Engine Documentation | Google Cloud Platform を参考にしています。

[プロジェクト]ページに移動し、"Computing Engineを試す"のチュートリアルに従うと、仮想マシンインスタンスの作成が体験できます。途中で無料トライアルを有効にするためクレジットカード情報などの登録が必要になります。

僕は一旦このトライアルで作成したインスタンスを削除した後、"r-studio"という名前のインスタンスを全く同じ設定で新たに作成しました。

R/RStudioのインストール

まず、VMインスタンスページからSSHボタンをクリックしてインスタンスのターミナルに接続し、以下の作業を行います。RとRStudio serverのインストールは、ローカルのマシンにインストールするのと同様に行います(Download RStudio Server – RStudio)。既に手元のサーバーマシンなどにRStudio serverをインストールしたことがある方には馴染み深い手順になっていると思います。

Rのインストール

まず、そのままapt-getを実行すると古いバージョン(2017/5/27現在 3.1.1 (2014-07-10))がインストールされてしまうため、apt-getのインストール元にCRANを追加する必要があります。エディタで /etc/apt/source.list を開き(sudo vi /etc/apt/sources.list)、以下の接続先を追加してください(CRAN mirror は CRAN - Mirrors から適当に選んで下さい。httpsではエラーとなってうまく行きませんでした)。

deb http://cran.cnr.berkeley.edu/bin/linux/debian jessie-cran34/

その後apt-get updateを実行すると、恐らくPublic keyが足りないというエラーで怒られます。

sudo apt-get update
# W: GPG error: http://cran.cnr.berkeley.edu jessie-cran34/ Release: The following signatures couldn't be verified be cause the public key is not available: NO_PUBKEY FCAE2A0E115C3D8A

Debian 7にRの実行環境(3.1.1)をapt-getでインストールする - Symfoware を参考に、末尾にある FCAE2A0E115C3D8A (選んだCRAN mirrorによって違うかも?) を使用して以下のように公開鍵を登録後、改めてapt-get update と apt-get install r-base を行います。

apt-key adv --keyserver keyserver.ubuntu.com --recv-keys FCAE2A0E115C3D8A
sudo apt-get update
sudo apt-get install r-base

コマンドライン版のRをRコマンドで起動し、最近のバージョンがインストールされていることを確認します。(2017/5/27現在 3.4.0 (2017-04-21))

RStudio serverのインストール

再び Download RStudio Server – RStudio を参考に、RStudioをインストールします

sudo apt-get install gdebi-core
wget https://download2.rstudio.org/rstudio-server-1.0.143-amd64.deb
sudo gdebi rstudio-server-1.0.143-amd64.deb

RStudio severにアクセス

RStudioにアクセスできるようにするため、ファイアーウォールの設定を追加する必要があります。

Google Cloud Platformのトップページのメニューから、"ネットワーキング>ファイアーウォールルール"に移動し、以下のルールを追加します( RStudio Server on a Google Compute Engine instance · Joe Roe を参考にしました)。
- 名前: allow-rstudio
- トラフィックの方向: 上り
- ソースIPの範囲: 0.0.0.0/0 (allow from any source)
- 指定したプロトコルとポート: tcp:8787
- ターゲットタグ: allow-rstudio
ターゲットタグであるallow-rstudioを、作成したVMインスタンスの"編集"の中にある"ネットワークタグ"に追加します。
http://<your-instance-ip>:8787 からRStudio serverに接続します。

RStudio serverにログイン

RStudio Server on a Google Compute Engine instance · Joe Roe
上のサイトの"Users and permissions"に詳細が書いてありますが、VMインスタンス上のデフォルトアカウントがkey-basedなログイン形式を用いている一方、RStudio serverはパスワードによるログインにしか対応していないため、RStudio server用のアカウントをLinux上に別に作成する必要があるようです。指摘されているように、既にデフォルトのアカウントでいろいろとファイルを作成してしまっていたら少し面倒そうです。

ともあれ、ブラウザ上のterminalからsudo adduser <username>によって新規ユーザーを作成します。その後、http://<your-instance-ip>:8787 にあるRStudio serverを開き、Linuxアカウントを用いてログインすると無事RStudio serverが使用できます。

`tidyverse`のインストール

ライブラリーのインストール中に幾つかのLinux用ライブラリが足りないと叱られました。以下のコマンドで足りないライブラリをインストールしたところ無事tidyverseがインストールできました。

sudo apt-get install libcurl4-openssl-dev libssl-dev libxml2-dev

元々のモチベーションであったRStanのインストールについては、以下に記載しています。

yoshidk6.hatenablog.com

ファイル転送

まだ試していませんが、VMインスタンスへのファイルの転送は通常のSCPによって行うことができるようです。
Transferring Files to Instances | Compute Engine Documentation | Google Cloud Platform

最後に

終わったらインスタンスを停止することを忘れないようにしましょう。

References

googleComputeEngineRという、ローカルのRから(?)GCEにアクセスしてサーバーをセットアップするパッケージ。Shiny serverを作成する等色々な用途に使うようなら便利かもしれません。今回は完全にブラウザから作成することを目指したので除外。 Launch RStudio Server in the Google Cloud with two lines of R · Mark Edmondson
Googleが提供しているCloud SDKを用いたRStudioのインストール法。これもローカルPCにソフトウェアのインストールが必要になるので除外しましたが、記事で使ったコマンドの多くはこのサイトのものを参考にしました。 Deploying R Studio on Compute Engine
もう一つ(日本語の)記事を見かけました。2. 同様にGoogleが提供しているコマンドラインツールを使用しています。 Google Container Engineの無限のリソースでRStudioを動かす - Technically, technophobic.

2017-05-24

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

Stan

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

f:id:yoshidk6:20170528160031p:plain

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%予測区間も、実測値と比較的良く合っています。 f:id:yoshidk6:20170528160259p:plain

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

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)、ラベルスイッチングを何らかの方法で避ける必要がありました。一被験者の場合は値に上下関係をつけることも可能なのですが、後にあるように複数被験者に対して各個人のパラメータを推定すると制限をつけるのが難しくなってしまうため、今回は初期値をうまく設定することにしました。

2017-05-22

StanでPopulation pharmacokinetic解析

Stan

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

項目(予定)

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

モチベーション

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

モデル

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

f:id:yoshidk6:20170522122333p:plain:w500

このモデルからは簡単に解析解を導出することができます。興味の有る方は↓の記事を御覧ください（上の図と対応させるには、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 に大いに影響されています。

使用する仮想データ

Toolchain

Configuration

Installing RStan

モチベーション

Virtual Machineをセットアップ

R/RStudioのインストール

Rのインストール

RStudio serverのインストール

RStudio severにアクセス

RStudio serverにログイン

tidyverseのインストール

ファイル転送

最後に

References

使用する仮想データ

Stanモデル

transformed parameters

model

Rコード

結果の評価

項目(予定)

モチベーション

モデル

解析コード・データ

Acknowledgement

`tidyverse`のインストール