今回は点滴静注投与の場合を解析してみます。
コンパートメントモデルの部分を変更しただけなので、モデル式さえ書ければいろいろな場合を解析できますね。
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(070518);
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;
if Time <= &Tin then do;
ext = (1 - exp(-keli * Time));
end;
else if Time > &Tin then do;
ext = ((1 - exp(-keli * &Tin)) * exp(-keli * (Time - &Tin)));
end;
Pred = (Dose / &Tin) / (V1i * keli) * ext;
LnConc = log(Pred) + sqrt(&s2) * rand('Normal');
Conc = exp(LnConc);
if Conc < &MLimit then do;
Pred = .;
LnConc = .;
Conc = .;
end;
if Time = 0 then Conc = 0;
output;
i + 1;
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;
/* 非線形混合効果モデルによる血中濃度の */
/* 1コンパートメントモデルによる当てはめ */
/* 誤差分布に対数正規分布を仮定 */
/* 点滴静注投与の場合 */
/* Laplace approximation */
proc nlmixed data = nlmixed1c(drop = pred where=(Conc ^= 0))
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);
if Time <= &Tin then do;
ext = (1 - exp(-kel * Time));
end;
else if Time > &Tin then do;
ext = ((1 - exp(-kel * &Tin)) * exp(-kel * (Time - &Tin)));
end;
Pred = (Dose / &Tin) / (V1 * kel) * ext;
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;
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;
/* 参考(誤差分布が正規分布) */
proc nlmixed data = nlmixed1c(drop = pred where=(Conc ^= 0))
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);
if Time <= &Tin then do;
ext = (1 - exp(-kel * Time));
end;
else if Time > &Tin then do;
ext = ((1 - exp(-kel * &Tin)) * exp(-kel * (Time - &Tin)));
end;
Pred = (Dose / &Tin) / (V1 * kel) * ext;
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;
さすがにモデルがかなり複雑になってきたので、推定に結構時間がかかるようになったと思います。
参考に載せたように、誤差分布が間違っているとそれなりにバイアスの入った推定値が得られることが分かりました。
グラフをしっかり見て適切な誤差分布を仮定して解析しなければ危ないという事でしょう。
まあ、実データでうまく当てはまる分布をどう選択するか等の問題があるので、基準を検討できれば研究テーマにもなるかと思います。
あとは経口投与の場合のみ、2 コン以降はモデル式の誘導程度ですませて適当にプログラムを乗っけようと思います。
