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共に体重と強い相関を持っていることがわかりました。
解析
上の結果に基いて、各パラメータへの体重の影響を定量的に評価します。
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の相関は、体重が共通の因子として影響している事による間接的な相関である、と考えられます。
今回のモデルを用いることで、各個人のCovariateに基いてパラメータが精度良く予測できるようになり、各個人の薬物動態プロファイルがより正確に予測できるようになります。