NLMIXEDプロシジャの紹介

岸本 淳司(慶応義塾大学SFC)

非線形混合モデル
M 非線形モデル– パラメタが非線形に入っているモデル

M 混合効果モデル– 固定（母数）効果とランダム（変量）効果とが同時に 入っているモデル

M 非線形混合モデル– 固定効果とランダム効果とが非線形に入っているモデル– たいへん複雑なモデルだ！
– なんか新しそうだ！（発想自体は昔からある）
2

1

応用例１：薬物動態
M 濃度の変化は非線形（コンパートメントモデル）
M 個人の効果をランダム効果とする
M 個人ごとにモデルを当てはめるのではなく、集団としてのモデルを推定する
M Population PK
M ブリッジングに使う！

3

応用例２：Overdispersion のある二項分布
M 二項分布では平均が決まると分散も決まる
M 実際には平均から期待されるより広い分散をもつデータが観測される（Overdispersion)
M 二項分布では、すべての個体は同一のパラメタを持つと仮定している
M 個体ごとにパラメタの違いを認めて、それが分布していると仮定すると、Overdispersionがうまく説明できる
42

NLMIXED指定の基本方針
M 説明変数の値を与えたときのデータの分布（条件付分布）を記述する– DATAステップに類似したプログラミングステートメン
トを使う

M 条件付分布には、正規分布はもちろん、二項分布やポアソン分布が、さらには対数尤度関数を指定することによって任意の分布が指定できる
M ランダム効果について分布を記述する– ランダム効果に指定できるのは正規分布のみ
5NLMIXED推定の基本方針
M ランダム効果について積分した最尤推定を行う– 非線形モデルではこの積分計算が厄介– ここに Adaptive Gaussian Quadrature を採用– NONMEM の1次法もサポート（正規分布のみ）
次法もサポート（正規分布のみ）

M 最適化計算のテクニックはマニアックに揃っている– SAS/OR の NLP プロシジャからコードを借用

M 柔軟に指定できることを重視– 定型的使用にはテンプレートの開発が必要

63

MIXEDプロシジャとの比較
M 線形混合モデルのREMLに相当する方法はな
い
M データの分布として正規分布以外の二項分布やポアソン分布も利用できる
M 誤差の分布形の指定（AR1とかCSなど）はできない
M 線形混合モデルの「ランダム係数モデル」の非線形への拡張である
7正規分布の例（オレンジの木）
yij = [(b1 + ui1)/(1 + exp[-(dij - b2 )/ b3 ])] + eij
M ロジスティック成長曲線の混合モデル– ｙij ： 測定値– ｄij ：日付
– b1,b2,b3 : 固定効果
– ui1 ： ランダム効果 N(0,σu 2)84

NLMIXEDの起動とパラメタ指定
proc nlmixed data=tree;parms b1=190 b2=700 b3=350
s2u=1000 s2e=60;
M NLINプロシジャと同様
M 初期値はグリッド指定可
9モデルの指定
num = b1+u1;
ex = exp(‑(day‑b2)/b3);
den = 1 + ex;
M 複数ステートメントに分けて記述できる
M u1はランダム効果を示すパラメタ
105

分布の指定model ynormal(num/den,s2e);
random u1normal(0,s2u) subject=tree;run;
M 条件付分布をmodelステートメントで指定– ~ 記号を使うSASプロシジャは初めて

M ランダム効果の分布をrandomステートメントで指定– ランダム効果が複数ある場合は共分散を指定
11

出力（問題のサイズ）
Dimensions
Observations Used
Observations Not Used
Total Observations
Subjects
Max Obs Per Subject
Parameters
Quadrature Points

35
035
575
1126

コメント
M Subjects の項目について– NLMIXEDにはCLASSステートメントがない– Subjectsでソートしてないといけない– ここで正しい数になっているかチェックすべし

M Quadrature Point について
– ここでは u は線形関数なので、積分点は１でよい

13

出力（反復過程）
IterationHistory
Iter Calls NegLogLike Diff
MaxGradSlope
14 131.686742 0.805045 0.010269 ‑0.633
26 131.64466 0.042082 0.014783 ‑0.0182
38 131.614077 0.030583 0.009809 ‑0.02796
410 131.572522 0.041555 0.001186 ‑0.01344
511 131.571895 0.000627
0.0002 ‑0.00121
613 131.571889 5.549E‑6 0.000092 ‑7.68E‑6715 131.571888 1.096E‑6 6.097E‑6 ‑1.29E‑6NOTE: GCONV convergence criterion satisfied.
147

コメント
M slopeがゼロに近くなったら収束
M １０-6 がデフォルトの基準

15

出力（当てはまりの情報）
Fit Statistics‑2 Log Likelihood
263.1
AIC (smaller is better)
273.1
BIC (smaller is better)
271.2
Log Likelihood‑131.6
AIC (larger is better)‑136.6
BIC (larger is better)‑135.6168

コメント
M この情報は、同じデータに対して複数のモデルを当てはめてみたとき、そのモデル間の比較に有用である

17

出力（パラメタ推定値）
Parameter Estimates
Standardt
Parameter Estimate
Error DF Value Pr > ￨t￨ Alpha Lower Upper Gradient
B1
B2
B3
S2u
S2e192.05 15.6473
727.90 35.2472
348.07 27.0790
999.88 647.4461.5139 15.8831

4 12.27 0.0003
4 20.65 <.0001
4 12.85 0.0002
4 1.54 0.1974
4 3.87 0.0179

0.05 148.61 235.50 1.154E‑60.05 630.04 825.76 5.289E‑60.05 272.88 423.25 ‑6.1E‑60.05 ‑797.70 2797.45 ‑3.84E‑60.05 17.4153 105.61 2.892E‑6189コメント
M 標準誤差は、2次導関数を使ったWald流
M 自由度は４– オレンジの木5本 － ランダム効果の次元１– t 分布とみなす– 正規分布とみなすより保守的

M 分散推定値の下限が負になっている– データの対数をとるとか平方根をとるともっとあてはまりがよくなるはず

19ロジスティックモデル・二項分布
・Overdispersionのモデル

xij ui ~ Binomial ( nij , pij )ηij = logpij
1 − pij= β 0 + β1t j + ui

M 固定効果：治療法 ｊM ランダム効果：施設 ｉ
20

10

ＮＬＭＩＸＥＤの起動と初期パラメタの指定
proc nlmixed data=infection;parms beta0=‑1 beta1=1 s2u=2;
M 二項分布だから残差のパラメタはない

21

モデルの指定eta = beta0 + beta1*t + u;
expeta = exp(eta);
p = expeta/(1+expeta);

2211分布の指定model xrandom ubinomial(n,p);normal(0,s2u)
subject=clinic;

M データの分布は二項分布を指定
M ランダム効果については、正規分布しか用意されていない
23

オプショナルな計算
predict eta out=eta;
estimate '1/beta1' 1/beta1;
run;
M 予測値の計算は Empirical Bayes によるM パラメタの関数については、あらゆる表現が可能
2412出力（問題のサイズ）
Dimensions
Observations Used
Observations Not Used
Total Observations
Subjects
Max Obs Per Subject
Parameters
Quadrature Points

16016823525

出力（パラメタ推定値）
Parameter Estimates
Standardt
Parameter Estimate Error DF Value Pr > ￨t￨ Alpha Lower Upper Gradient
Beta0‑1.1974 0.5561 7 ‑2.15 0.0683 0.05 ‑2.5123 0.1175 ‑3.1E‑7
Beta10.7385 0.3004 7
2.46 0.0436 0.05 0.02806 1.4488 ‑2.08E‑6S2u1.9591 1.1903 7
1.65 0.1438 0.05 ‑0.8554 4.7736 ‑2.48E‑7

M ランダム効果S2uの標準誤差はOverdispersionを考慮している– 普通のLogistic回帰分析ならもっと狭いはず

2613

ESTIMATEステートメントの効果
Label1/beta

Additional Estimates
Standardt
Estimate Error DF Value Pr > ￨t￨ Alpha Lower Upper
11.3542 0.5509 7
2.46 0.0436 0.05 0.05146 2.6569

M この推定値の標準誤差はδ法で計算されている

27

コンパートメントモデルの例
Cit = [(D k ei kai)/(Cli(kai - kei))] [exp(-k ei t) - exp(-kai t)]+ e itCli = exp(β1+bi1)
kai = exp(β2 +bi2)
kei = exp(β3 )C:濃度 D:用量 k:定数 Cl:クリアランス
β：固定効果 b：ランダム効果
28

14

プログラム例
proc nlmixed data=theoph;parms beta1=‑3.22 beta2=0.47 beta3=‑2.45
s2b1=0.03 cb12=0 s2b2=0.4 s2=0.5;
cl = exp(beta1 + b1);
ka = exp(beta2 + b2);
ke = exp(beta3);
pred = dose*ke*ka*(exp(‑ke*time)‑exp(‑ka*time))/cl/(ka‑ke);model conc
normal(pred,s2);
normal([0,0],[s2b1,cb12,s2b2])
random b1 b2
subject=subject;
run;
29

結論
M 時代の要求に応える新しい道具が用意された
M モデルの扱いが柔軟で、応用範囲が広い– そのため、定型的作業のためにはテンプレートの開発が必要

M 使いこなすためにはかなり勉強する必要がある

30

15

16