SUGI25 at Indianapolis
参加報告
Advanced Tutorialのトピックスの１つ
PROC NLMIXEDの紹介
国立がんセンター研究所
がん情報研究部
石塚直樹

NLMIXEDプロシジャとは？
• 何のためのプロシジャ？
「非線型混合モデル：nonlinear mixed models」
をあてはめる• 必要とされる状況は？
薬物動態、
Over dispersionのある二項分布のデータ

1

今までのプロシジャとは?
• 非線型回帰モデルの解析
NLINプロシジャがあったけど、• 混合モデルの解析
MIXEDプロシジャ→これを一般化するとNLMIXED• GLIMMIX、NLINMIXマクロが近い→ 推定方法が違う• 非線型の最適化→ SAS/ORのNLPプロシジャが近い

統計的な特徴は?
• 変量効果について積分した尤度の近似を最大化する(最尤法)• 積分の近似法
adaptive Gaussian quadratureが使える
first-order Taylor series approximation• 最大化の方法デフォルト dual quasi-Newton algorithm.

2統計的な特徴

つづき• 標準誤差
対数尤度関数の2次導関数の行列に基づいた
Waldタイプの近似、t 値は推定値／標準誤差• 予測値の計算
変量効果に経験ベイズ推定値を用いる方法が可能• 固定効果の任意の関数の推定
デルタ法により利用可能

数値例

その１

• オレンジの木：Draper and Smith(1981)• ロジスティック成長曲線の混合モデル：Lindstrom and Bates (1990), Pinheiro and Bates (1995)

b1 + ui1
yij =+ eij
1 + exp[− ( d ij − b2 ) / b3 ]
yij：i番目の木のj番目の測定値
dij：i番目の日

b1, b2, b3：固定効果
ui1：変量効果 N(0,σu2)
eij：誤差 N(0,σe2)

3

SASプログラム
未知パラメータの初期値の指定

proc nlmixed data=tree;parms b1=190 b2=700 b3=350 s2u=1000 s2e=60;
num = b1+u1;
モデルの左辺の分子
ex = exp(-(day-b2)/b3);
den = 1 + ex;
モデルの左辺の分母model y ~ normal(num/den,s2e);
random u1 ~ normal(0,s2u) subject=tree;run;
変量効果の分散成分

誤差の分散成分

SASアウトプット(1)• Specification
入力データセット、従属変数、変量効果の変数、それぞれの分布、個体(被験者)を示す変数、最適化の方法• Dimension
オブザベーション数、個体数、パラメータ数• Parameters
未知パラメータの初期値• Iteration Histroy-log尤度 の最適化の経過

4

SASアウトプット(2)• Fit Statistics
あてはめの良さ（尤度）に関する情報• Parameter Estimates
未知パラメータの最尤推定値
標準誤差：ヘシアン行列に基づく
自由度：個体数−変量効果のパラメータ数
分散成分(s2u, s2e)には解釈に注意が必要
この例では変量効果が線型なので積分の近似法は1次のテーラー展開でも同じ結果になる

数値例

その２

• 感染症治療薬の多施設共同試験(8施設、2剤、応答は2値)：Beitler and Landis(1985)• ロジスティック混合モデル

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

nij：施設iで治療法jの患者数、xij：治癒した患者数
tj：治療法j，ui：施設の変量効果N(0, σu2)5SASプログラム
proc nlmixed data=infection; 変量効果の分散成分
parms beta0=-1 beta1=1 s2u=2;
eta = beta0 + beta1*t + u;
expeta = exp(eta);
p = expeta/(1+expeta);
二項分布のmodel x ~ binomial(n,p);モデルの指定random u ~ normal(0,s2u) subject=clinic;
predict eta out=eta;
estimate '1/beta1' 1/beta1;run;
ηijの経験ベイズ推定値を1/β1 の推定含むデータセットetaの作成

SASアウトプット? Additional Estimate
ESTIMATEステートメントで指定した結果
Parameter
beta0
beta1

Estimate-1.1974
0.7385

StandardError
0.55610.3004

DF77t Value
-2.15
2.46

Pr > |t|0.06830.0436

推定値1/0.7385=1.3542, 標準誤差0.3004/0.73852=0.5509
Label1/beta1

Estimate
1.3542StandardError
0.5509DF7t Value
2.46

Pr > |t|0.04366

SASアウトプット つづき
信頼区間
Parameter

Alpha

Lower

Upper

beta0
beta1

0.05
0.05

-2.51230.02806

0.11751.4488線型の変換ではない場合に、信頼限界から直接変換したほうが良いかもしれない

デルタ法では、
1/0.02806=12.4069≠2.6569
1/1.4488=0.6923≠0.05146

LabelAlpha

Lower

Upper

1/beta1

0.05

0.051462.6569数値例

その３• 薬物動態の解析(テオフィリン経口投与)：Pinheiro and Bates(1995)• １コンパートメント・モデル

Cit =Dkei kai
Cl i[exp(− k t ) − exp(− k t )]+ eeiaiit

Cit：被験者iの時間tの血中濃度
D ：用量kei：被験者iの消失速度定数
kai：被験者iの吸収速度定数
Cli：被験者iのクリアランス
eit ：被験者iの時間tの正規分布に従う誤差

7

数値例

その３ つづき• 個体差を表す変量効果bi1, bi2の導入

Cit = exp(β1 + bi1 )kai = exp(β 2 + bi 2 )kei = exp(β 3 )  0  s 2b1 cb12   bi1   ~ N   ,  b0cb12
s2b2
 i2   

SASプログラム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);
random b1 b2 ~ normal([0,0],[s2b1,cb12,s2b2])
subject=subject;run;

8MODELステートメントで指定可能な分布
normal(m,v)：平均m、分散vの正規分布
binary(p)：確率pの二値(ベルヌーイ)分布
binomial(n,p)：カウントn、確率pの二項分布
gamma(a,b)：形状パラメータa、尺度パラメータbのガンマ分布
negbin(n,p)：カウントn、確率pの負の二項分布
poisson(m)：平均mのポアソン分布
general(ll)：PROGRAMMINGステートメントで指定した対数尤度関数llに従う分布開発担当者Dr.Wolfingerに質問してきました• NONMEMと同じことができるか？
積分の近似方法をMETHOD=FIROとすることで、1次のテーラー展開を用いれば可能• GENMOD(GEE)と同じ結果が得られるか？
Subject-specific(NLMIXED)に対して、Population-avaraged(GEE)という違いがある• MIXEDプロシジャと同様にrepeatedステートメントは利用可能か？
使えない：非線型混合モデルでは未だ一般的ではない

9

質問のつづき• 初期値の設定方法は？
指定のない時はすべて初期値は１→ gridサーチも可能 ex. parms -10 to 10 by 5;
注：Dr.Wolfingerは初期値によらず収束すると言うが、パラメータ推定値が異なることが多々ある!!
→ 収束条件を変えて確認？ (GCOV, XCOVオプション)
最適化の方法にオプションを加えるか？
(LINESEARCH=オプション, UPD=オプション)• 分散成分の自由度の計算は?
個体(被験者)数−変量効果のパラメータ数、あるいはオブザベーション数で固定効果と同じ
ただし、任意の自由度をDF=d として指定可能

質問のつづき• パラメータに非負の制約があるときは？
BOUNDSステートメントで範囲の指定が可能
ex. bounds 0 <= a1-a9 x <=1, b1-b10 y>=0;• ３つ以上の変量効果の指定方法は？
３つの変量効果 b1,b2,b3の場合
  0   s11 s12 s13  b1     b2  ~ N   0 ,  s 21 s 22 s 23  0   s31 s32 s33 b  3  

ただし、 s 21 = s12, s 23 = s 32, s31 = s13, s 32 = s23
random b1 b2 b3 ~ normal([0,0,0],[s11,s21,s22,s31,s32,s33]) subject=person

10質問のつづき• 施設と被験者のような二重に枝別れした指定は？
指定できない。ただし、新しいNLINMIXマクロ(nlmm801.sas)で可能• ２群のパラメータの比較方法は？
PROGRAMMINGステートメントで、群毎に別の変数を指定し(CLASSステートメントは使えない!) 、ESTIMATEステートメントで差を指定する
例えば、
if (group=1) then eta = alpha1 + beta1*logtstd + e;else eta = alpha2 + beta2*logtstd + e;
estimate ‘alpha1-alpha2’ alpha1-alpha2;
estimate ‘beta1-beta2’ beta1-beta2;

質問のつづき• 変量効果の正規性のチェック方法は？
PREDICTステートメントで経験ベイズ推定値を出力して確認することが可能• パラメータが多すぎないかチェックする方法は？
極端な場合にはヘシアン行列がランク落ちするパラメータ間の共分散、相関係数が大きくなる→パラメータ推定値の分散共分散行列を出力(COVオプション)→パラメータ推定値の相関行列を出力(CORRオプション)→ヘシアン行列の出力(HESSオプション)

11最後に教科書は？
• Davidian, M. and Giltinan, D.M. (1995)
Nonlinear Model for repeated Measurement Data,
New York: Chapman & Hall.
• Vonesh, E.F. and Chinchihill, V.M. (1997)
Linear and Nonlinear Models for the Analysis of
Repeated Measurements,
New York: Marcel Dekker.
• Agresti, Booth, Hobert, Caffoのカテゴリカルデータについての論文
http://www.stat.ufl.edu/~aa/

やっぱりSASは凄かった!!
• NLMIXEDプロシジャの開発者のDr.
WolfingerはMULTTESTプロシジャ、MIXEDプロシジャの開発も担当• Nonlinear Mixed ModelについてBiometrikaなどに彼自身も論文を書いている• アカデミアでの実績だけではなく、研究成果を広く一般に、多くの人が使えるプロシジャにして提供するってただ者ではない

12