静注の1-コンパートメントモデル。
周辺尤度を求める積分の近似にはラプラス近似を使用、[1]を見ると Adaptive Gaussian quadrature (適応的ガウス求積法) が平均に対するパラメータ推定に関しては誤差が少ないように見える。
分散に対するパラメータの推定については、FO (Taylor展開一次近似) だとロバスト分散を使っても悪そうですね (ロジスティックモデルもコンパートメントモデルも)。
収束の成功割合を考えるとラプラス近似辺りがリーズナブルなのかも。
[1]は1995年の論文なので、新しめな手法も含めた比較はされていないみたい。
/* パラメータ等の設定 */ %let kel = 0.3; %let V1 = 60; /* 個体間変動と測定誤差 */ %let e2_kel = 0.1; %let e2_V1 = 0.1; %let s2 = 0.01; %let MLimit = 0.001; /* 濃度の測定限界 */ /* ダミーデータの作成 */ data nlmixed1c; do Dose = 200; do ID = 1 to 50; Lnkeli = log(&kel) + sqrt(&e2_kel) * rand('Normal'); LnV1i = log(&V1) + sqrt(&e2_V1) * rand('Normal'); keli = exp(Lnkeli); V1i = exp(LnV1i); do Time = 0, 0.1, 0.15, 0.25, 0.5, 1, 1.5, 2, 4, 6, 8, 10, 12, 14, 16, 20, 24; Pred = Dose / V1i * exp(-keli * Time); LnConc = log(Pred) + sqrt(&s2) * rand('Normal'); Conc = exp(LnConc); if Conc < &MLimit then do; Pred = .; LnConc = .; Conc = .; end; output; end; end; end; run; /* グラフにより確認 */ proc gplot data = nlmixed1c; plot Conc * Time = ID / nolegend haxis = axis1; plot LnConc * Time = ID / nolegend haxis = axis1; symbol i = j v = none c = black r = &N; axis1 order = (0 to 24 by 1); run; quit; /* Laplace approximation */ proc nlmixed data = nlmixed1c(drop=pred) method = gauss qpoints = 1 cov hess; parms Lnkel -1 LnV1 4 e2_kel 0.1 e2_V1 0.1 s2 0.01; kel = exp(Lnkel + eta_kel); V1 = exp(LnV1 + eta_V1); Pred = Dose / V1 * exp(-kel * Time); model Conc ~ normal(Pred, s2); random eta_kel eta_V1 ~ normal( [0, 0], [ e2_kel, 0 , e2_V1 ] ) subject = ID; predict kel out = Pred_kel; predict V1 out = Pred_V1; predict e2_kel out = Pred_e2_kel; predict e2_V1 out = Pred_e2_V1; predict Pred out = Pred; run;
積分の近似オプション
method = value
FIRO: Taylor展開の一時近似
GAUSS: 適応的ガウス求積法(qpoints=1を加えるとラプラス近似、数理は[1])
ISAMP: Importance Sampling、モンテカルロ積分の一種。尤度関数が期待値の形になる性質を利用して、効率よく乱数を発生させて積分を近似する。異常に時間がかかって悩ましい。積分→パラメータ更新→積分→・・・・
random statement
SAS 9.1.3は変量効果が従う確率分布に多変量正規分布以外を指定できない。
random statementは1回しか指定できない、ネストしたモデルは無理。
RANDOM random-effects ~ distribution SUBJECT=variable <options>;
random-effects: プログラム上で指定した変量効果の変数名
distribution: normal([平均ベクトル], [分散共分散行列])
subject: 変量効果のグループ変数(個体とか施設とか)
optionsはHelpを参照。
分散共分散行列の指定は、下三角行列っぽく書く。
きっとSASがバージョンアップするといろんな分布が使えるようになるんだろう。
predict statement
変量効果の経験ベイズ推定ができます、[4]とSAS Helpを参考に。
EBLUP(Empirical Best Linear Unbiased Predictor)のMSEには条件付きMSEで何かをやるみたいな事が書いてあるため、SASの経験ベイズ推定は条件付きMSEを用いているみたい、そのうち[4]を読もう。
predict expression OUT=SAS-data-set <options>;
expressionにはSASの簡単な関数を書ける、exp(x1)とかx1+x2とか。
optionはHelpを参照。
FOCE
NONMEMのFOCEに対応する方法は[3]にある%NLINMIX macroを使ってできる模様?
EBLUPによる変量効果の推定も可能、そのまま使ったらエラーでた、怖い怖い。
ダミーデータ
今回発生させているダミーデータが少々気になる。
個別のグラフプロットを見ると濃度が上がったりしていて、実際には濃度が上がることはほとんど無いはずなので、何かしらの対策を考えなくては。
BUGS
MCMC をやる場合は BUGS がいい。
無料だけど登録が必要なので、ここから登録して、メールで送られてくる key を読み込む。
Population PK/PD については PKBUGS がある。ちょっとだけ触った。
[2]を読みたいのでメモ。
参考文献
[1] Pinheiro, J.C. Bates, D.M. Approximations to the Log-likelihood Function in the Nonlinear Mixed-effect Model. Journal of Computational and Graphical Statistics 1995; 4: 12–35.
[2] J. Wang. EM algorithms for nonlinear mixed effects models. Computational Statistics & Data Analysis 2007; 51(6): 3244–3256.
[3] SAS Customer Support Center. Sample 539: %NLINMIX macro to fit nonlinear mixed models [cited Apr. 21, 2007]. http://support.sas.com/ctx/samples/index.jsp?sid=539&tab=details
[4] Booth, J.G. and Hobert, J.P. Standard Errors of Prediction in Generalized Linear Mixed Models. Journal of the American Statistical Association 1998; 93: 262–272.