今日は indicator function を導入して、測定濃度限界を考慮したモデルによる解析をしてみよう。
indicator function については文献[1]に書いてあります。
/* パラメータ等の設定 */
%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;
%let V1 = 60;
/* 個体間変動と測定誤差 */
%let e2_kel = 0.1;%let e2_V1 = 0.1;
%let s2 = 0.01;
%let MLimit = 0.001; /* 濃度の測定限界 */
/* パラメータの対数変換 */
data _null_;
call symput ('Lnkel', log(&kel));
call symput ('LnV1', log(&V1));
run;
/* ダミーデータの作成 */
data nlmixed1c;
call streaminit(070428);
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);
Lim = 0;
if Conc < &MLimit then do;
Lim = 1;
Pred = .;
LnConc = .;
Conc = 0; /* MissingはNLMIXEDが弾くため */
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.1 e2_V1 0.1
s2 0.01;
kel = exp(Lnkel + eta_kel);
V1 = exp(LnV1 + eta_V1);
Pred = Dose / V1 * exp(-kel * Time);
if Lim = 0 then do;
likelihood = 1 / (sqrt(2 * constant('PI') * s2)) *
exp(-1 * (Conc - Pred) ** 2 / (2 * s2));
logL = log(likelihood);
end;
else do;
logL = logcdf('normal', &Mlimit, Pred, sqrt(s2));
end;
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)
method = gauss qpoints = 1 cov hess;
where Lim = 0;
parms Lnkel -1 LnV1 4
e2_kel 0.1 e2_V1 0.1
s2 0.01;
kel = exp(Lnkel + eta_kel);
V1 = exp(LnV1 + eta_V1);
Pred = Dose / V1 * exp(-kel * Time);
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;
ダミーデータの測定値が高い、要は測定誤差の分散が異常に高いと収束しないで吹っ飛ぶので乱数のシードを固定しています。
CDF関数を入れると
NOTE: A finite difference approximation is used ....
という出力がログに出る、うーむ数値解析を勉強しなきゃ。
ダミーデータの作り方もそのままだ、GW中に考えよう。
参考文献
[1] 矢船 明史, 石黒 真木夫, 赤池 弘次. 母集団薬物データの解析. 朝倉書店 2004.
履歴
ミス発見。
でも修正したら収束しなくなった・・・
とりあえず pdf, cdf 系の関数を NLMIXED Procedure 中につっこむとどうなるかを検討してみた。
indicator function を考慮しないモデルで、logpdf 関数を使った場合と正規分布の対数尤度関数をベタ打ちした場合を比較。
/* logpdf関数 */
logL = logpdf('normal', Conc, Pred, sqrt(s2));
model i ~ general(logL)
/* 対数尤度関数ベタ打ち */
logL2 = log(1 / (sqrt(2 * constant('PI') * s2)) *
exp(-1 * (Conc - Pred) ** 2 / (2 * s2)));
model i ~ general(logL2);
logL を使うと反復計算すら始まらない、パラメータの初期値をセットした段階でエラー。
logL2 を使うと model statement の normal を使用した場合と全く同じ結果が得られる。
一応 Data Step で両方の対数尤度関数の値を計算してみると、1.0E-14程度の差しかなかったため、関数の引数が間違っているわけではなさそうだった。
思いつきでやってみたものの、cdf 関数を使って indicator function を入れるのは危なそうだなぁ。
