血中濃度のシミュレーションデータの問題について

問題点

• 血中濃度のデータを想定しているくせに、前の時点よりも濃度が高くなってしまうシミュレーションデータが生成されてしまう

根本的な解決になるかどうかは分からないが、シミュレーションデータを作ってみた。

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

/* パラメータ等の設定 */
%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;

/* 前時点濃度を上回らない様にしたダミーデータ */
data nlmixed1c;
  call streaminit(070525);
  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;

goptions reset = all;
goptions ftext = 'Times New Roman' ftitle = 'Times New Roman'
         htitle = 1.6 htext = 1.6; * hsize = 15cm vsize = 15cm;
options linesize=120 pagesize = 300;

title 'Truncated';
proc gplot data = nlmixed1c;
  plot Conc * Time = ID / nolegend haxis = axis1 vaxis = axis2;
  plot LnConc * Time = ID / nolegend haxis = axis1 vaxis = axis2;
  symbol i = j v = none c = black r = &N;
  axis1 order = (0 to 24 by 2);
  axis2 label = (a = 90) minor = none;
run; quit;

/* 何も処理しないダミーデータ */
data nlmixed1cn;
  call streaminit(070525);
  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;
        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;
        i + 1;
      end;
    end;
  end;
run;
/* グラフにより確認 */
title 'Not truncated';
proc gplot data = nlmixed1cn;
  plot Conc * Time = ID / nolegend haxis = axis1 vaxis = axis2;
  plot LnConc * Time = ID / nolegend haxis = axis1 vaxis = axis2;
  symbol i = j v = none c = black r = &N;
  axis1 order = (0 to 24 by 2);
  axis2 label = (a = 90) minor = none;
run; quit;

グラフを描いてみると、前の時点の濃度を超えてしまうようなデータは生成されていない事が確認できた。

グラフを見ると濃度が過小になりそうですが、果たして推定精度とか妥当性とかはどうなるんでしょうか？

今日は眠いのでまた今度。