Using non-vectorized function in dplyr::mutate and how to associate column names for purrr::pmap

map family functions of the purrr package are very useful for using non-vectorized functions in dplyr::mutate chain (see GitHub - jennybc/row-oriented-workflows: Row-oriented workflows in R with the tidyverse or Beware of Vectorize · Jim Hester's blog).
I encounter the needs for this especially when dealing with nested data frames.

One of the drawbacks is that name/input argument assignments become confusing when you want to use more than two columns of your data frames (and using pmap family) for the function of interest. This post first briefly review how mutate works in combination with map or map2, then provide two approaches to avoid confusions around name assignments when using pmap.

How mutate works with vectorized functions

In most cases, the processes you want to do in mutate is vectorized and there is no need to use map family function. This works because the output from the function of interest (c in the example below) has the same length as the original data frame, and mutate only need to append one column to the data frame.

library(tidyverse)

df0 <- tibble(a = 1:3, b = 4:6)

df0 %>% mutate(c = a + b)

Non-vectorized function with one or two input arguments (map or map2)

Imagine that we want to create a new column containing arithmetic progressions in each row [ref (in Japanese)]. Since seq function is not vectorized, we cannot directly use this in mutate chain.

df1 <- tibble(a = c(1, 2), b = c(3, 6), c = c(8, 10))

df1 %>% mutate(d = seq(a, b))
# Error in mutate_impl(.data, dots) : Evaluation error: 'from' must be of length 1.

Instead, we can use map family function here. map family function take list(s) as input arguments and apply the function of interest using each element of the given lists. Because each column of data frames in R is a list, map works very well in combination.

In this example, we want to provide two input arguments to the seq function, from and to. map2 is the appropriate function for this.

df2 <- df1 %>% mutate(d = map2(a, b, seq))

as.data.frame(df2)
#  a b  c             d
#1 1 3  8       1, 2, 3
#2 2 6 10 2, 3, 4, 5, 6

The figure below shows how map function handles this process in mutate chain.

f:id:yoshidk6:20180806153334p:plain

Like you do with map function outside mutate, we can use map_dbl or map_chr to create columns with double or character types.

If we want to explicitly specify names of the argument, .x and .y can be used. See what happens with this:

df2 <- mutate(df1, d = map2(a, b, ~seq(.y, .x)))

as.data.frame(df2)

Non-vectorized function with three or more input arguments (pmap)

Assignment of column names become confusing when using three or more columns, because we don't have shorthand like .x or .y any more. Let's take a look at the following example using rnorm function.

Case example

Generate a list of random numbers for each row with rnorm function. Each row of the original data frame contain different value of mean, sd, n.

We will first prepare a data frame with columns corresponding to mean, sd, n, and apply rnorm function for each row using pmap. Each element of the new column data contains a vector of random samples *1. This type of structure is called as "nested data frames" and there are many resources on this, such as R for Data Science.

A simple case

If your data frame has the exact same names and numbers of columns to the input arguments of the function of interest, a simple syntax like the one below works *2.

df4 <- 
  tribble(~mean, ~sd, ~n,
          1,  0.03, 2,
          10, 0.1,  4,
          5,  0.1,  4)

df4.2 <- 
  df4 %>% 
  mutate(data = pmap(., rnorm)) 

as.data.frame(df4.2)

One caution is that the syntax like the one below doesn't work. pmap thinks that you are calling rnorm(df4$n, df4$mean, df4$sd) for each row, and each element of the new column contain three random samples from the same list of mean and sd. *3

df4 %>% mutate(data = pmap(., ~rnorm(n, mean, sd))) %>% as.data.frame() # Wrong answer

Number of columns > Number of input arguments

In most cases, however, you will have more columns than the input arguments. pmap complains in this case, saying that you have unused argument.

df5 <- 
  tribble(~mean, ~sd, ~dummy, ~n,
          1,  0.03, "a", 2,
          10, 0.1,  "b", 4,
          5,  0.1,  "c", 4)

df5 %>% mutate(data = pmap(., rnorm))  # Error

There are two ways to avoid this error.

Make a small list on the fly

The first method is to create a small list that only contains the necessary columns (Ref: Dplyr: Alternatives to rowwise - tidyverse - RStudio Community )

df5.2 <- 
  df5 %>% 
  mutate(data = pmap(list(n=n, mean=mean, sd=sd), rnorm))

as.data.frame(df5.2)

Here, list(n=n, mean=mean, sd=sd) create a new list with three vectors named n, mean, and sd, which serves the same purpose as the df4 data frame in the above example.

Mind that if you don't give names to the elements of the new list, the order of the list items will be used to associate with input arguments of rnorm. My recommendation is to always assign names to the list elements.

df5 %>% mutate(data = pmap(list(n, mean, sd), rnorm)) # Correct but not recommended
df5 %>% mutate(data = pmap(list(mean, sd, n), rnorm)) # Wrong answer

Use ... to ignore unused columns

The second method is to absorb unused columns with ... (Ref: Map over multiple inputs simultaneously. — map2 • purrr). A syntax like the one below works because pmap automatically associate names of the input list and names in function(). In other word, columns names of the data frame must match the variable names in the function().

df5.3 <- 
  df5 %>% 
  mutate(data = pmap(., function(n, mean, sd, ...) rnorm(n=n, mean=mean, sd=sd))) 

as.data.frame(df5.3)

Input arguments of function() and rnorm() are not automatically associated with names. It is recommended to explicitly associate input argument name for the function of interest (rnorm in this case).

df5 %>% mutate(data = pmap(., function(n, mean, sd, ...) rnorm(n, mean, sd))) # Correct but not recommended
df5 %>% mutate(data = pmap(., function(n, mean, sd, ...) rnorm(mean, sd, n))) # Wrong answer
df5 %>% mutate(data = pmap(., function(mean, sd, n, ...) rnorm(mean, sd, n))) # Wrong answer

A syntax like the one below gives unexpected outputs, as you saw in the df4 example.

df5 %>% mutate(data = pmap(., function(...) rnorm(n=n, mean=mean, sd=sd))) # Wrong answer

Column names are different from the input argument names

You can use either of the two approaches above.

df6 <- 
  tribble(~mean1, ~sd1, ~dummy, ~n1,
          1,  0.03, "a", 2,
          10, 0.1,  "b", 4,
          5,  0.1,  "c", 4)

df6.2 <-
  df6 %>% mutate(data = pmap(list(mean=mean1, sd=sd1, n=n1), rnorm)) 
as.data.frame(df6.2)

df6.3 <- 
  df6 %>% mutate(data = pmap(., function(n1, mean1, sd1, ...) rnorm(n=n1, mean=mean1, sd=sd1)))
as.data.frame(df6.3)

*1:In the examples below (and above), we further use as.data.frame function to exposure actual numbers of vectors

*2:This works even if the order of the columns is different from the order of input arguments

*3:This happens because rnorm is actually vectorized. See ?rnorm: If length(n) > 1, the length is taken to be the number required.

dplyr::mutate内でベクトル化されていない関数を使う - purrr::pmapでのカラム名の対応付け

Rのデータフレームの各行に対して処理を行いたいが、関数をベクトル化することが難しい場合*1には、下で簡単に書くようにmutatemapを組み合わせると非常に便利です。しかし、map内で定義する関数の入力値が3つ以上になると話が少しややこしくなります。納得するまで結構混乱したので、自分の備忘録も兼ねて記事にしておくことにしました。

通常のmutateを使った処理

殆どの場合、mutate内で行う処理(例えば+*ifelse等)はベクトル化されているので、出力値(下のケースではc)が他の列と同じ長さになり、元のデータフレームに追加されます。

library(tidyverse)

df0 <- tibble(a = 1:3, b = 4:6)

df0 %>% mutate(c = a + b)

各行について別々に処理、入力値が2つまでの場合

こちらにあるように、データフレームの各行に対して等差数列を作成する例を考えてみます。 このとき、seq関数はベクトル化されていないため、以下のように書くとエラーになります。

df1 <- tibble(a = c(1, 2), b = c(3, 6), c = c(8, 10))

df %>% mutate(d = seq(x, y))
# Error in mutate_impl(.data, dots) : Evaluation error: 'from' must be of length 1.

ここでmap系関数の出番です。
map系関数は、リストを引数とし、リスト内の各要素に対して与えられた関数を適用した結果を返すという働きを持っています。 Rのデータフレームは各列がそれぞれリストとなっているので、map系巻数と相性が良くなっています。

今回の例ではseq関数にfromtoの2つの引数を与えたいので、map2関数を使います。

df2 <- df1 %>% mutate(d = map2(a, b, seq))

print(df2$d)

下にこの処理がどのように行われているかのイメージ図を貼ってみました。

f:id:yoshidk6:20180806153334p:plain

もちろん、使用したい処理が一つの引数のみを取る場合はmapを使えばいいですし、作成する列をdoublecharacter型にしたければmap_dblmap_chrなどを使います。 また、mapmap2では、関数内で明示的に引数を指定したい場合、.x.yを使います。例えば、上の例でdf2 <- mutate(df1, d = map2(a, b, ~seq(.y, .x)))を実行してみてください。

各行について別々に処理、入力値が3つ以上の場合

この様にmap系関数を使うと、ベクトル化されていない処理をmutate内で簡潔に書くことが出来ます。 3つ以上のカラムを使用したい場合、同様にpmapを使うことになるのですが、カラム名と関数の引数名の対応付けが少し複雑になります。

以下の例を用いて、どのように対応付けを行うことができるのか見ていきます。

使用する例

複数の平均値・標準偏差の組み合わせについて、それぞれ異なった数の乱数をrnorm関数を用いて算出したい。

ここではまず、meansdnに対応する3つのカラムを持ったデータフレームを作成し、各行ごとにrnorm関数を適用して乱数を発生させています。新しく生成されたカラムdataには、各行ごとに乱数値を保持したベクトルが格納されます。
結果を見やすくするため、この例では更にunnnest()関数を用いて縦長のデータフレームに変換していますが、変換せずにそのままmap系関数を使って処理を続けていくことも可能です。
ネスト化周りに関してはR4DSなど色々と解説があると思うので、馴染みのない方はそちらも参考にしてみてください。

単純な例

データフレーム内のカラムが、pmap内で使用する関数の引数と個数・名前共に一致する場合にのみ、以下のような簡潔な文法が通じます。

df4 <- 
  tribble(~mean, ~sd, ~n,
          1,  0.03, 2,
          10, 0.1,  4,
          5,  0.1,  4)

df4 %>% mutate(data = pmap(., rnorm)) %>% unnest()

ちなみに、下のような書き方だと、rnorm内の入力はdf4$n df4$mean df4$sdであると見做されてしまい、想定していない結果となってしまいます。

df4 %>% mutate(data = pmap(., ~rnorm(n, mean, sd))) %>% unnest() # Wrong answer

余分なカラムが含まれている場合

実際には、カラム数が関数の引数の数と一致するという状況は非常に稀だと思います。 余計なカラムが含まれていると、unused argumentがあると言って怒られます。

df5 <- 
  tribble(~mean, ~sd, ~dummy, ~n,
          1,  0.03, "a", 2,
          10, 0.1,  "b", 4,
          5,  0.1,  "c", 4)

df5 %>% mutate(data = pmap(., rnorm))  # Error

これを回避するには、2つの方法があります。

小さいリストを作ってしまう

1つ目は、pmapの第一引数をdf4のようなリストにその場で変えてしまう方法です。(参照: Dplyr: Alternatives to rowwise - tidyverse - RStudio Community )

df5 %>% mutate(data = pmap(list(n=n, mean=mean, sd=sd), rnorm)) %>% unnest()

注意点として、下のようにリスト内要素に名前をつけないと、rnormへの受け渡しが順序のみによって行われてしまいます。少し面倒ですが、上のようにきちんと名前を割り振る書き方を推奨します。

df5 %>% mutate(data = pmap(list(n, mean, sd), rnorm)) %>% unnest() # Correct but not recommended
df5 %>% mutate(data = pmap(list(mean, sd, n), rnorm)) %>% unnest() # Wrong answer

... で未使用変数を受け流す

2つ目は、pmapへの入力ではデータフレームをそのまま渡す一方で、関数宣言の所で未使用の変数を受け流す方法です(参照: Map over multiple inputs simultaneously. — map2 • purrr, "Use ... to absorb unused components of input list .l")。
pmapが、入力したリスト名と関数の引数名を自動で照合してくれる性質を利用しているので、function()内の変数名とデータフレームのカラム名が一致する必要があります。
引数を二回書く必要が出てきますが、リスト内要素に名前をつけ忘れる心配はありません。

df5 %>% mutate(data = pmap(., function(n, mean, sd, ...) rnorm(n=n, mean=mean, sd=sd))) %>% unnest()

ただしもちろん、普段関数を使う際と同様に、引数を名前で紐付けないと誤った結果となります。

df5 %>% mutate(data = pmap(., function(n, mean, sd, ...) rnorm(n, mean, sd))) %>% unnest() # Correct but not recommended
df5 %>% mutate(data = pmap(., function(n, mean, sd, ...) rnorm(mean, sd, n))) %>% unnest() # Wrong answer
df5 %>% mutate(data = pmap(., function(mean, sd, n, ...) rnorm(mean, sd, n))) %>% unnest() # Wrong answer

先にdf4の例で見たように、以下のような書き方も誤った結果となります。

df5 %>% mutate(data = pmap(., function(...) rnorm(n=n, mean=mean, sd=sd))) %>% unnest() # Wrong answer

カラム名がpmap内で使用する関数の引数名と一致しない場合

上2つのいずれの方法を使っても対応できます。

df6 <- 
  tribble(~mean1, ~sd1, ~dummy, ~n1,
          1,  0.03, "a", 2,
          10, 0.1,  "b", 4,
          5,  0.1,  "c", 4)

df6 %>% mutate(data = pmap(list(mean=mean1, sd=sd1, n=n1), rnorm)) %>% unnest()
df6 %>% mutate(data = pmap(., function(n1, mean1, sd1, ...) rnorm(n=n1, mean=mean1, sd=sd1))) %>% unnest()

map系関数の使い方について参考になるサイト

*1:特にnested dafa frameを扱っているときに頻発します

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-コンパートメントモデルを使ってみます。

インストール

必要なもの

手順

最近までインストールが複雑でしたが、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の間で投与前濃度しか測定を行っていないため、それ以外の時点のシミュレーション結果が出力されていないことに依ります*2f:id:yoshidk6:20171030065007p:plain

その他

今回の例ではまずTorstenライブラリを使ってみることに主眼をおいたため、単一被験者のデータ解析という非常にシンプルなものになっています。 ポピュレーション解析をに興味がある方は、同じexample-models内にあるTwoCptModelPopulationを見てみることをおすすめします。

参照したリンク先のまとめ

*1:省略できるよう改良を計画しているようです。

*2:RとStanコードを改良することで、測定点以外の時間についてもシミュレーション結果を表示することが可能です。そのうちそれについて記事を書くかもしれません。

Stan & POPPK (4): Covariateの影響を考慮した解析

Stan & POPPK (3): 多変量正規分布を用いたパラメータ間の相関まで含めた解析 - yoshidk6’s blog の続きです。パラメータの個人間変動を生む要因を分析してモデルに組み込みます。

パラメータの個人間変動を生じる因子の探索

前回までの解析で、薬物動態パラメータであるCLとVDに個人間変動が見られました。 個人間変動の要因を見つけ出すことができれば、新しい被験者に薬物を投与する際に、その被験者がどのような薬物動態特性を示すかをより高い精度で予測できるようになります。

今回のデータセットでは、↓にあるように各個人の体重と性別が与えられているので、推定された各個人のパラメータ値との関連を見てみます。 https://github.com/yoshidk6/simple_pk_stan/blob/master/data/subj_dose_cov_20170521.csv

# 推定パラメータ値の抽出
fit.CLi <- 
  summary(fit.stan, pars = c("CLi"))$summary %>% 
  tbl_df() %>% 
  select(mean)
fit.VDi <- 
  summary(fit.stan, pars = c("VDi"))$summary %>% 
  tbl_df() %>% 
  select(mean)

fit.indiv.params <-
  bind_cols(tibble(ID=1:nrow(fit.CLi)),
            fit.CLi %>% rename(CL=mean),
            fit.VDi %>% rename(VD=mean))

fit.indiv.params %>% 
  ggplot(aes(CL,VD)) +
  geom_point() +
  geom_smooth(method="lm")


# 体重・性別データと結合し、パラメータ値との相関をプロット
data.cov <- read_csv("../data/subj_dose_cov_20170521.csv") 

fit.indiv.params.cov <-
  full_join(fit.indiv.params, data.cov)

fit.indiv.params.cov %>% 
  gather(Parameter, Value, CL, VD) %>% 
  ggplot(aes(WT,Value)) +
  geom_point() +
  facet_wrap(~Parameter, scales="free") +
  geom_smooth(method="lm")

fit.indiv.params.cov %>% 
  gather(Parameter, Value, CL, VD) %>% 
  ggplot(aes(factor(SEX),Value)) +
  geom_boxplot() +
  facet_wrap(~Parameter, scales="free")

下図より明らかに、CLとVD共に体重と強い相関を持っていることがわかりました。 f:id:yoshidk6:20171028120238p:plain f:id:yoshidk6:20171028120240p:plain

解析

上の結果に基いて、各パラメータへの体重の影響を定量的に評価します。

Stanモデル

本モデルでは、各個人のパラメータ値を算出する際(57-60行)に体重の影響を組み込んでいます。影響度合いを記述するためには、薬物動態界隈で頻用されているexponential modelを使用しました。各個人のパラメータは、[個人の体重 (WT[k])/体重の中央値 (WTMED)]の指数によって影響されると仮定し、その係数(WTCLとWTVD)を推定しています。

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

Rコード

上記のモデルを動かすRのコードは以下の通りです。 サンプリングが一部不安定になってしまったため、thinningを導入しています。

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")
data.cov  <- read_csv("../data/subj_dose_cov_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)),
         sigma = exp(rnorm(2, log(0.2), 0.5)),
         rho = diag(2),
         s_Y = runif(1, 0.5, 2),
         WTCL= rnorm(1, 0.5, 1),
         WTVD= rnorm(1, 0.5, 1),
         CLVDiLog=matrix(rep(log(c(0.5,5)), ea = nrow(data.subj)), 
                      nrow = 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,
       WT   = data.cov$WT,
       WTMED= median(data.cov$WT),
       Y    = data.pk$CONC)


nChains <- 4
nPost <- 1000 ## Number of post-burn-in samples per chain after thinning
nBurn <- 1000 ## Number of burn-in samples per chain after thinning
nThin <- 10

nIter <- (nPost + nBurn) * nThin
nBurnin <- nBurn * nThin

fit.stan <-
  stan(file = "mod_04_cov.stan", 
       data = data, init=init,
            iter = nIter,
            warmup = nBurnin,
            thin = nThin, 
            chains = nChains,
            control = list(adapt_delta = 0.8))

結果の評価

推定の結果、観測値は正しく再現されました(図は繰り返しになるので省略します)。

WTCLはおおよそ0.75、WTVDはおおよそ1.07と推定されました。体重が1.5倍になるとCLとVDはそれぞれ1.35倍、1.5倍になるということになり、比較的大きな影響を持っていると言えます。

体重の影響を除いた後のパラメータの個人間変動を見てみると、前回のモデルと比べて標準偏差が非常に小さくなり(CLとVDでそれぞれ0.2, 0.3から0.06, 0.1に減少)、CLとVDの相関も0.9程度からほぼ0にまで低下しています。パラメータの個人間変動の推定値(赤)は、体重を考慮した予測値(青)によって正しく再現されています。更に、体重の影響を除いた上での予測値(緑)を見てみると、ばらつきが非常に小さいことが見て取れます。以上より、(1)観測された個人間変動は体重によって殆どが説明され得る、(2)前回見られたCLとVDの相関は、体重が共通の因子として影響している事による間接的な相関である、と考えられます。 f:id:yoshidk6:20171028122921p:plain

今回のモデルを用いることで、各個人のCovariateに基いてパラメータが精度良く予測できるようになり、各個人の薬物動態プロファイルがより正確に予測できるようになります。

Stan & POPPK (3): 多変量正規分布を用いたパラメータ間の相関まで含めた解析

Stan & POPPK (2): 階層モデルを用いた複数被験者の薬物動態解析 - yoshidk6’s blog の続きです。前回使用したモデルを拡張して、パラメータ間の相関も含めた解析を行います。

使用する仮想データ

作成された仮想データは前回と同じなので、説明を省略します。

解析

Stanモデル

モデル構造は前回とほぼ同じですが、多変量正規分布を扱うための変更を行っています。
Stan超初心者入門 の114ページ目から詳しい説明が載っています。
簡単にまとめると、個人間変動を考慮している2つのパラメータに対して相関行列rhoと各々の標準偏差sigmaを定義し、そこから分散共分散行列Omegaを算出するようになっています。

GitHubリンク: https://github.com/yoshidk6/simple_pk_stan/blob/master/code/mod_03_multinorm.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)),
         sigma = exp(rnorm(2, log(0.2), 0.5)),
         rho = diag(2),
         s_Y = runif(1, 0.5, 2),
         CLVDiLog=matrix(rep(log(c(0.5,5)), ea = nrow(data.subj)), 
                      nrow = 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_03_multinorm.stan", 
       data = data, init=init)

結果の評価

推定されたパラメータを元に算出された平均値と90%予測区間を見てみると、前回と同様に、全体・個人ごとのプロット共に観測値とよく合っています。

f:id:yoshidk6:20171027125833p:plain f:id:yoshidk6:20171027125837p:plain

一方で、推定された各個人のパラメータと、モデルに基づいた予測値を比較すると、CLiVDiの相関が正しく表現されるようになったことがわかります。(一つ目がパラメータ間の相関を考慮していない前回の結果、2つめが今回の結果)

f:id:yoshidk6:20171027130526p:plain f:id:yoshidk6:20171027130025p:plain

今回のモデルを用いることで、正しく個人間のパラメータ変動を反映した予測結果が算出されることが期待できます。

O-1ビザ取得顛末記

自分がO-1ビザの申請をした時、研究者・科学者の申請に関する情報がインターネット上にほとんど無く難儀したので、備忘録的に体験談を書いておきます。*1

O-1を申請することになった経緯

外国人がアメリカの企業で働くには、就労可能なビザが必要になります。 アメリカ企業による正規雇用の場合*2、最もメジャーなのはH1Bビザを取ることです。しかし、H1Bビザに申し込めるチャンスは年に一回しかなく、申請から発給までにもしばらく時間がかかります。年間発給数も決まっているため、抽選に外れると次年まで待たなければなりません。多くの人は、アメリカの大学・大学院を卒業するともらえるOPTというシステムを使って働き始め、有効期間の1~2年のうちにH1B取得を目指すようです。

僕の場合、日本の大学院を出た後にポスドクを経由して就職したので、これには該当しません。その上、ポスドク時代に取得していたJ-1ビザに"Two-Year Home-Country Physical Presence Requirement"*3 という制限がかかっていたため、そもそもH1Bを取得することができませんでした*4

幸いなことに、O-1というビザには

  1. いつでも申し込め、審査プロセスも早い
  2. Two-year ruleの対象外であり、Waiver手続きをしなくても取得できる

という利点があったため、僕の場合はこのカテゴリーでのビザ取得を目指すことになりました。このビザはIndividuals with Extraordinary Ability or Achievementの為とされており、"卓越した能力"を示すためにかなりの量の書類を用意する必要があります。比較的特殊なケースになるので弁護士費用なども高くなり、このビザのサポートをしてくれる企業も少なくなります。ただ、申請のハードルは思ったより高くなく、移民局のページに書いてあるようなノーベル賞を持っているようなレベルの科学者である必要は全く無いようです。

書類の用意

O-1を申請するにあたって何よりも重要なのは、科学者として卓越していることを示すための書類集めです。 USCIS(移民局)のページによると、以下のカテゴリーから少なくとも3つを満たす必要があるとされています。

  1. Receipt of nationally or internationally recognized prizes or awards for excellence in the field of endeavor
  2. Membership in associations in the field for which classification is sought which require outstanding achievements, as judged by recognized national or international experts in the field
  3. Published material in professional or major trade publications, newspapers or other major media about the beneficiary and the beneficiary’s work in the field for which classification is sought
  4. Original scientific, scholarly, or business-related contributions of major significance in the field
  5. Authorship of scholarly articles in professional journals or other major media in the field for which classification is sought
  6. A high salary or other remuneration for services as evidenced by contracts or other reliable evidence
  7. Participation on a panel, or individually, as a judge of the work of others in the same or in a field of specialization allied to that field for which classification is sought
  8. Employment in a critical or essential capacity for organizations and establishments that have a distinguished reputation

弁護士の人と話した感触だと、ある程度の論文数+引用数(4と5に該当?)をベースとした上で、その他をできれば2つ以上満たしたいという所でした。引用数はどれくらいがminimumなのかは分かりませんが、2桁でも大丈夫そうです(二桁前半・後半など細かくは分かりません)。特許も考慮されるようですが、どのカテゴリーに該当するのかはいまいち分かりません。僕の場合は結局、論文・引用・学会発表に加え、海外学振&国内学会のポスター賞を1に、論文誌のPeer-reviewerの経験を7に適用しました。O-1の取得を目指す可能性のある方は、カテゴリー数を増やすことも念頭に置いて、積極的にReviewerの役割を受けることをおすすめします。

併せて、申請には専門領域で確立された研究者からの複数の推薦状が必須です。僕は五通、(1)複数の国の人から and/or (2)自分の論文の論文を引用してくれている人から集めてくれと言われました。更に、そのうち少なくとも2通は今まで直接仕事・研究をしたことがない人からである必要がありました。最終的に、2通ずつをアメリカと日本の方に、1通をイギリスの方にお願いしました。

その他、僕が経験した・聞いたことのうち、役に立ちそうな事を列挙しておきます。

  • 6に関しては、本当に飛び抜けた給料である必要があるようで、通常の会社勤めの人が満たすのは難しそうです。
  • Citationのカウントやリストの作成にはGoogle Scholar Citationsを使いました(そのまま印刷して提出しました)。
  • 僕は使っていないのでどうなのか分かりませんが、Web上で検索すると、O-1が通る確率を無料で評価してくれるようなサービスが有るようです。
  • 貢献を数字として示しやすいIFやjournalの分野内でのランキング(Web of Scienceで調べられる)は、自分の論文・引用元の論文・Reviewerの経験のすべてにおいて重視されていました。

経過

ビザが手に入るまでの時間はおおよそ三ヶ月となりましたが、人によってかなりばらつきがあるようです(半年近くかかったという話も聞いたことがあります)。

  • 11月中旬: 弁護士と連絡を取り始め、申請に用いる書類のカテゴリーを模索(およそ二週間)
  • 11月下旬~12月上旬: Reference lettersを除く必要書類の準備の完了
  • 12月中旬~1月上旬: Reference lettersも含めた必要書類の準備の完了*5
  • 1月下旬: USCISへのPetition申請
  • (2月上旬: 日本に一時帰国)
  • 2月上旬: Petitionの承認(申請から約一週間強*6 )
  • 2月上旬: ↑の数日後、I-797等面接に必要な書類の原本を弁護士からFEDEXで受け取る
  • 2月中旬: 日本の米国大使館で面接 *7
  • 2月中旬: ビザの貼付されたパスポートを取得(面接から一週間強)

結びに

上にも書きましたが、O-1取得のハードルは思ったほど高くないので、アメリカ就職を考えている科学者の方(特にポスドクからアメリカに来られている方)には一考の価値があると思います。

*1:もちろん必要な書類や手続の詳細などは個人・申請年によって違うので、詳細はケースを扱ってくれる弁護士の人にご相談下さい。

*2:インターンシップやポスドク等の場合はJ-1というビザを使うことが多いと思います。

*3:J-1プログラムの終了後、二年間母国に物理的に住まないと、多くのアメリカビザ・グリーンカードが取得できないという制限。政府系機関から給料を貰っていたりした場合に該当します。

*4:一応Waiver手続きをすることは可能なのですが、それについてはWeb上に色々と情報があるので見てみて下さい。

*5:ホリデーシーズンにかかってしまったため、少し手間取りました。

*6:Premium processingという、追加料金を払う代わりに15 calendar days以内に返事をもらえるという制度を使用しました。

*7:J-1の二年ルールに該当している人はアメリカ国内でビザをO-1に切り替えることが出来ないらしく、一時帰国して申請することになりました。第三国(カナダなど)でも一応は大丈夫なことになっているようですが、母国のほうが安全なようです。

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に触ったことのある方は、興味がある場所から見てみて下さい。

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

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

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

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

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