同じネタで 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 アルゴリズムという方法が流行っているようです。
