今回は誤差分布に対数正規分布を仮定したモデルによる解析を行ってみようと思います。
依然シミュレーションデータの問題は時間の都合上修正していませんが、アイデアは浮かんだので近いうちに直す予定です。
シミュレーションデータは誤差分布に対数正規分布を仮定して発生させているので、解析も仮定通りにやってみます。
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 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;
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;
/* グラフにより確認 */
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)
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;
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;
誤差分布に対数正規分布を仮定する場合、測定誤差の分散よりも変量効果の分散が大きい場合に収束しにくい感じがしました。
このあたりは数をこなして経験を積むか、時間があるときにちゃんと検討しておくべきなんだろうと思います。
シミュレーションデータ修正の指針は、1 ショット静注ならば濃度が一つ前の時点より高くならなければいいので、濃度が一つ前の時点を超えてしまったら生成し直す形で発生させるようにしようと思います。
こうしてしまうと、きちんと対数正規分布に従わず、誤差の高い値が切断されてしまうので、そのあたりはシミュレーションで確認してまとめてみようかと思います、適当にねｗ
