読者です 読者をやめる 読者になる 読者になる

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

SAS

同じネタで 5 回も引っ張ってしまった。


何はともあれ、シミュレーション用のプログラムができました。

dm 'log; clear; output; clear';
proc datasets kill; run;
option linesize = 130 pagesize = 500;

%let execpath = " ";
%let Path = " ";
%macro setexecpath;
  %let execpath = %sysfunc(getoption(sysin));
  %if %length(&execpath) = 0 %then
    %let execpath = %sysget(sas_execfilepath);
  data _null_;
    do i = length("&execpath") to 1 by -1;
      if substr("&execpath", i, 1) = "\" then do;
        call symput("Path", substr("&execpath", 1, i));
        stop;
      end;
    end;
  run;
%mend setexecpath;
%setexecpath;

libname Out "&Path";

/* パラメータ等の設定 */
%let Dose = 200;
%let N = 50;
%let MTime = 0, 0.1, 0.15, 0.25, 0.5, 1, 1.5, 2, 
             4, 6, 8, 10, 12, 14, 16, 20, 24;
%let Tin = 0.5;

%let kel = 0.3; * log(kel) = -1.204;
%let V1 = 60; * log(V1) = 4.094;

/* 個体間変動と測定誤差 */
%let e2_kel = 0.01; %let e2_V1 = 0.01;
%let s2 = 0.05;

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

/* パラメータの対数変換 */
data _null_;
  call symput ('Lnkel', log(&kel));
  call symput ('LnV1', log(&V1)); 
run;

%macro sim(start, rep);

proc datasets lib = work;
  delete Out;
run;

ods exclude all;
%do i = %eval(&start) %to %eval(&start + &rep - 1);

proc datasets lib = work;
  delete Est;
run;

/* 前時点濃度を上回らない様にしたダミーデータ */
data nlmixed1c;
  call streaminit(&i);
  retain BConc 1E10;
  i = 1;
  do Dose = &Dose;
    do ID = 1 to &N;
      Lnkeli = &Lnkel + sqrt(&e2_kel) * rand('Normal');
      LnV1i  = &LnV1 + sqrt(&e2_V1) * rand('Normal');
      keli = exp(Lnkeli);
      V1i = exp(LnV1i);
      do Time = &MTime;
        do until (Conc < BConc);
          Pred = Dose / V1i * exp(-keli * Time);
          LnConc = log(Pred) + sqrt(&s2) * rand('Normal');
          Conc = exp(LnConc);
        end;
        BConc = Conc;
        if Conc < &MLimit then do;
          Pred = .;
          LnConc = .;
          Conc = .;
        end;
        output;
        i + 1;
      end;
      BConc = 1E10;
    end;
  end;
  drop BConc;
run;

/* 非線形混合効果モデルによる血中濃度の */
/* 1コンパートメントモデルによる当てはめ */
/* 誤差分布に対数正規分布を仮定 */
/* Laplace approximation */
proc nlmixed data = nlmixed1c(drop=pred)
  method = gauss qpoints = 1 cov hess;
  parms Lnkel -1 LnV1 4
        e2_kel 0.01 e2_V1 0.01
        s2 0.05;
  kel = exp(Lnkel + eta_kel);
  V1 = exp(LnV1 + eta_V1);
  Pred = Dose / V1 * exp(-kel * Time);
  likelihood = 1 / (sqrt(2 * constant('PI') * s2) * Conc) *
               exp(-1 * (log(Conc) - log(Pred)) ** 2 / (2 * s2));
  logL = log(likelihood);
  model i ~ general(logL);
  random eta_kel eta_V1 ~ normal(
    [0, 0],    [
      e2_kel,
      0     , e2_V1
    ]
  )
  subject = ID;
  ods output ParameterEstimates = Est;
run;

data Est;
  set Est;
  seed = &i;
run;

proc append base = Out data = Est; run;

%end;
ods select all;

/* ME = (Predicted - True) / N */
/* MSE = (Predicted - True) ^ 2 / N */
data Out.Out2;
  set Out;
  if Parameter = 'LnV1' then ME = Estimate - log(&V1);
  if Parameter = 'Lnkel' then ME = Estimate - log(&kel);
  if Parameter = 'e2_V1' then ME = Estimate - &e2_V1;
  if Parameter = 'e2_kel' then ME = Estimate - &e2_kel;
  if Parameter = 's2' then ME = Estimate - &s2;
  MSE = ME ** 2;
run;

proc sort data = Out.Out2; by Parameter; run;
proc means data = Out.Out2;
  var ME MSE;
  by Parameter;
run;

%mend sim;

%sim(070530, 1000);

とりあえず動作確認で 2 回だけ回しました。
まだ何とも言えないですね、2 回だし。
1,000 回やってみて、いい結果出てたら 10,000 回ぐらいやってみようかな。

追記

2011-06-02
結局このやり方ではいい結果になりませんでした。
しかし、NLMIXED procedure を使って Population PK の解析をする方法をまとめることはできたかなと思います。


非線形混合効果モデルの推測方法は、

  • 周辺尤度に基づく推測
  • EM アルゴリズムによる全尤度に基づく推測

があり、固定効果部分に関しては基本的にルート n 一致性のある推定量が得られます。
近似手法の場合はルート n 一致性が無いものがあり、FOCE 系の手法は一番収束が遅く (n と p (個人ごとの測定回数) が両方大きくないといけない)、次いでラプラス近似が収束が早いはずです (これも n と p が両方大きくないといけない)。
尤度を周辺化する積分をまじめに計算するとクソ重いので、SAEM アルゴリズムという方法が流行っているようです。