Triad sou.

NLMIXED Procedureによる非線形混合効果モデルの当てはめ (1)

SAS

静注の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による変量効果の推定も可能、そのまま使ったらエラーでた、怖い怖い。

ダミーデータ

今回発生させているダミーデータが少々気になる。
個別のグラフプロットを見ると濃度が上がったりしていて、実際には濃度が上がることはほとんど無いはずなので、何かしらの対策を考えなくては。

R

R だと lmer という関数があるようだ。
実際さわっていないからどんな解析が出来るのかは不明だが、PQL とか ML とか MCMC とかができるみたい。

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.