NLIN Procedureで1-コンパートメントモデルのパラメータ推定

NLINで非線形モデルの推定をしてみよう。
測定誤差は正規分布を仮定、薬物動態パラメータの自然対数をとったもの推定します。
NLINを使ったのは、Marquardt法を使ってみたかったから。NLMIXEDでやると最尤法。

/* 静注1-コンパートメントモデルの当てはめ */
%let kel = 0.3;
%let V1 = 60;

/* 個体間変動と測定誤差 */
%let eta_kel = 0.1;
%let eta_V1 = 0.1;
%let s2 = 0.01;

/* 濃度の測定限界 */
%let MLimit = 0.001; 

/* ダミーデータ */
data nlin1c;
  do Dose = 200;
    do ID = 1 to 50;
      Lnkeli = log(&kel) + sqrt(&eta_kel) * rand('Normal');
      LnV1i  = log(&V1) + sqrt(&eta_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 = nlin1c;
  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;

/* 各個体の血中濃度の1コンパートメントモデルによる当てはめ */
proc nlin data = nlin1c method = marquardt noitprint;
  parms Lnkel -5 to 5 by 1
        LnV1 -5 to 5 by 1;
  kel = exp(Lnkel);
  V1 = exp(LnV1);
  Pred = Dose / V1 * exp(-kel * Time);
  model Conc = Pred;
  by ID;
  output out = Estimate p = Predict residual = Residual
         SSE = SSE Parms = Lnkel LnV1;
run;

PK関係はダミーデータの作り方がちょっとめんどくさい。