このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

nlmefitsa

確率的な EM アルゴリズムを使用した非線形混合効果モデルの近似

構文

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0)
[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0,'Name',Value)


説明

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0) は、非線形混合効果の回帰モデルを当てはめ、BETA に固定効果の推定を返します。既定の設定では、nlmefitsa は、各モデル パラメーターが対応する固定効果と変量効果の和となり、変量効果の共分散行列が対角行列、すなわち、無相関の変量効果となるモデルを当てはめます。

この関数が返す BETA、PSI、およびその他の値は、パラメーターの最尤推定値に収束するように設計されているランダム (モンテカルロ) シミュレーションの結果です。結果はランダムであるため、結果に対するシミュレーションのプロットを調べて、シミュレーションが収束したことを確認してください。また、複数の開始値を使用して関数を複数回実行したり、'Replicates' パラメーターを使用して複数のシミュレーションを実行することもお勧めします。

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0,'Name',Value) は、コンマで区切られたパラメーターの名前と値のペアを 1 つまたは複数受け入れます。一重引用符で囲んで Name を指定します。

入力引数

定義

次に示す引数のリストでは、以下の変数定義が適用されます。

• n — 観測数

• h — 予測子変数の数

• m — グループ数

• g — グループ固有の予測子変数の数

• p — パラメーター数

• f — 固定効果の数

X	h 個の予測子変数上の n 個の観測の n 行 h 列の行列。
Y	n 行 1 列の応答ベクトル。
GROUP	m 個のグループのどれに各観測値が属するかを示すグループ化変数。GROUP は、カテゴリカル変数、数値ベクトル、各行がグループ名に対応する文字行列、string 配列、または文字ベクトルの cell 配列が可能です。
V	データ内の m 個のグループごとの g 個のグループ固有の予測子変数の m 行 g 列の行列。これらは、グループ内のすべての観測に対して同じ値を取る予測子値です。V の行の順序は、GRP2IDX(GROUP) に従います。グループ固有の予測子値のサイズがグループ全体で変化する場合は、V に m 行 g 列の cell 配列を使用します。グループの予測子がない場合は、V に [] を指定します。
MODELFUN	予測子値とモデル パラメーターを受け入れ、近似値を返す関数のハンドル。MODELFUN は、YFIT = MODELFUN(PHI,XFUN,VFUN) の形式で、入力引数は次のとおりです。 PHI — モデル パラメーターの 1 行 p 列のベクトル。 XFUN — 予測子変数の l 行 h 列の配列。ここで、次のようになります。 l は 1 (XFUN が X の単一行の場合) l は ni (XFUN がサイズ ni の単一グループに対する X の行を含む場合) l は n (XFUN が X のすべての行を含む場合) VFUN — 次のいずれかです。 V の単一行に対応する単一グループのグループ固有の予測子の 1 行 g 列のベクトル。 n 行 g 列の配列。ここで、次の場合、k 行目の VFUN は V (i,:) です。k 番目の観測がグループ i に含まれる場合。 V が空の場合、nlmefitsa は 2 つの入力で MODELFUN を呼び出します。 MODELFUN は、近似値 YFIT の l 行 1 列のベクトルを返します。PHI または VFUN のいずれかが単一行を含む場合、その行は他の 2 つの入力引数内のすべての行に対応します。MODELFUN が 1 回の呼び出しでモデル パラメーターの複数のベクトルに対して YFIT を計算することが可能な場合、パフォーマンスを向上するには、'Vectorization' パラメーターの名前と値のペア (後述) を使用してください。
BETA0	f の固定効果に対する初期推定をもつ f 行 1 列のベクトル。既定の設定では、f はモデル パラメーターの数 p に等しくなります。BETA0 は、f 行 REPS 列の行列にすることもできます。この場合、推定は REPS 回繰り返され、BETA0 の各列が一連の開始値として使用されます。

名前と値の引数

既定の設定では、nlmefitsa は、各モデル パラメーターが対応する固定効果と変量効果の和となるモデルを当てはめます。固定効果または変量効果の数が異なるモデル、またはこれらの効果に依存するモデルを当てはめるには、以下のパラメーターの名前と値のペアを使用します。最大で、'FE' 接頭辞をもつ 1 つのパラメーター名と 'RE' 接頭辞をもつ 1 つのパラメーター名を使用します。下記で詳しく説明するように、選択肢によっては、nlmefitsa による MODELFUN の呼び出し方法が変わることに注意してください。

FEParamsSelect	1:p 内で要素をもつ数値ベクトル、または 1 行 p 列の logical ベクトルとして、固定効果を含むモデル パラメーターのベクトル PHI の要素を指定するベクトル。モデルは f 個の固定効果を含みます。ここで、f は、指定した要素数です。
FEConstDesign	p 行 f 列の計画行列 ADESIGN。ここで、ADESIGN*BETA は PHI の p 要素の固定要素です。
FEGroupDesign	m 個のグループごとに異なる p 行 f 列の固定効果の計画行列を指定する p x f x m 配列。
REParamsSelect	1:p 内で要素をもつ数値ベクトル、または 1 行 p 列の logical ベクトルとして、変量効果を含むモデル パラメーターのベクトル PHI の要素を指定するベクトル。モデルは r 個の変量効果を含みます。ここで、r は、指定した要素数です。
REConstDesign	p 行 r 列の計画行列 BDESIGN。ここで、BDESIGN*B は PHI の p 要素のランダム要素です。この行列は 0 と 1 で構成され、各行の 1 は最大で 1 つでなければなりません。

既定のモデルは、FEConstDesign と REConstDesign の両方を eye(p) に設定するか、FEParamsSelect と REParamsSelect の両方を 1:p に設定することと同じです。

追加のオプション パラメーターの名前と値のペアは、尤度を最大にするために使用する反復アルゴリズムを制御します。

CovPattern	変量効果の共分散行列 PSI のパターンを定義する r 行 r 列の logical 行列または数値行列 PAT を指定します。nlmefitsa は、PAT の非対角要素の非ゼロに対応する共分散と同様に、PSI の対角に沿った分散に対する推定を計算します。nlmefitsa は、残りの共分散 (PAT 内の非対角要素 0 に対応する) を 0 に制約します。PAT は、ブロック対角行列の行と列の置換でなければなりません。そのため、nlmefitsa は、必要に応じて非ゼロ要素を PAT に追加してそのようなパターンを生成します。PAT の既定値は、無相関の変量効果に対応する eye(r) です。 または、1:r の値を含む 1 行 r 列のベクトルとして PAT を指定します。この場合、等しい値をもつ PAT の要素は変量効果のグループを定義します。nlmefitsa は、グループ内でのみ共分散を推定し、グループにまたがる共分散を 0 に制約します。
Cov0	共分散行列 PSI の初期値。r 行 r 列の正定値行列でなければなりません。空の場合、既定値は BETA0 の値に依存します。
ComputeStdErrors	係数推定値の標準誤差を計算し、出力として STATS 構造体に保存する場合は true、この計算を省略する場合は false (既定値)。
ErrorModel	誤差項の形式を指定する文字ベクトルまたは string スカラー。既定値は 'constant' です。各モデルは、標準正規 (ガウス) 変数 e、関数値 f、1 つまたは 2 つのパラメーター a と b を使用して、誤差を定義します。以下から選択できます。 'constant' — y = f + a*e 'proportional' — y = f + b*f*e 'combined' — y = f + (a+b*f)*e 'exponential' — y = f*exp(a*e)、または等価な log(y) = log(f) + a*e このパラメーターが指定された場合、出力される STATS.errorparam フィールドの値は以下のようになります。 a ('constant' および 'exponential' の場合) b (次の場合) 'proportional' [a b] (次の場合) 'combined'
ErrorParameters	誤差モデルのパラメーターの開始値を指定するスカラーまたは 2 要素ベクトル。ErrorModel パラメーターに応じて、a、b、または [a b] の値を指定します。
LogLikMethod	対数尤度を近似する方法を指定します。選択肢は、以下のとおりです。 'is' — 重要度のサンプリング 'gq' — ガウス求積法 'lin' — 線形化 'none' — 対数尤度の近似を省略 (既定の設定)
NBurnIn	初期バーンイン反復回数。この間、パラメーター推定は再計算されません。既定の設定は 5 です。
NChains	シミュレートされた "連鎖" の数 c。既定値は 1 です。c > 1 に設定すると、各反復で各グループについて c 個のシミュレートされた係数ベクトルが計算されます。既定値はデータによって異なり、連鎖全体で約 100 グループを与えるよう選択されます。
NIterations	反復回数。これは、スカラーまたは 3 要素ベクトルです。アルゴリズムの次の 3 つの各段階に対して実行する反復回数を制御します。 シミュレーテッド アニーリング 完全なステップ サイズ 縮小されたステップ サイズ 既定の設定は [150 150 100] です。スカラーは、既定値として同じ比率で 3 段階に分散されます。
NMCMCIterations	マルコフ連鎖モンテカルロ (MCMC) の反復回数。これは、スカラーまたは 3 要素ベクトルです。主な反復の各段階で、次の 3 つのタイプの MCMC 更新が実行される回数を制御します。 完全な多変量更新 単一の座標更新 複数の座標更新 既定の設定は [2 2 2] です。スカラー値は、すべての要素がスカラーに等しい 3 要素ベクトルとして扱われます。
OptimFun	尤度関数を最大化する推定プロセスの最適化関数。'fminsearch' (既定) または 'fminunc' として指定します。'fminunc' を使用するには、Optimization Toolbox™ が必要です。関数 fminsearch は、関数評価のみを使用する直接探索法を使用します。関数 fminunc (Optimization Toolbox) は勾配法を使用し、一般に、尤度関数を最大化する最適化問題ではこの関数の方が効率的です。
Options	statset への呼び出しで作成される構造体。nlmefitsa では、以下の statset パラメーターを使用します。 DerivStep — 有限差分の勾配の計算に使われる相対微分。スカラーまたは長さがモデル パラメーター p 数であるベクトル。既定の設定は eps^(1/3) です。 Display — 推定時の表示のレベル。 'off' (既定) — 何の情報も表示しません。 'final' — 推定アルゴリズムの最後の反復の後に情報を表示します。 'iter' — 各反復で情報を表示します。 FunValCheck 'on' (既定値) — MODELFUN の無効な値 (NaN や Inf など) をチェックします。 'off' — このチェックをスキップします。 OutputFcn — @ を使用して指定された関数ハンドル、関数ハンドルをもつ cell 配列、または空配列。nlmefitsa は、各反復の後、すべての出力関数を呼び出します。出力関数の例は、nlmefitoutputfcn.m (nlmefitsa の既定の出力関数) を参照してください。 Streams — 関数 nlmefitsa 内で乱数の生成に使用するストリーム。RandStream オブジェクトとして指定します。既定の設定は、現在のグローバル ストリームです。 TolX — 推定された固定効果と変量効果の終了許容誤差。既定の設定は 1e-4 です。
ParamTransform	p パラメーターごとに変換関数 f() を指定する p 値のベクトル。 XB = ADESIGN*BETA + BDESIGN*B PHI = f(XB) ベクトルの各要素は、対応する PHI 値に変換を指定する次のいずれかの整数コードでなければなりません。 0: PHI = XB (すべてのパラメーターの既定値) 1: log(PHI) = XB 2: probit(PHI) = XB 3: logit(PHI) = XB
Replicates	ベクトル BETA0 の開始値から開始する推定の数 REPS。BETA0 が行列である場合は、REPS は、BETA0 の列数と一致しなければなりません。既定の設定は、BETA0 の列数です。
Vectorization	MODELFUN への PHI、XFUN、VFUN の入力引数の可能なサイズを定義します。利用可能な値は以下のとおりです。 'SinglePhi' — MODELFUN は一度に 1 つのモデル パラメーター セットに対してのみ YFIT を計算できる関数 (ODE ソルバーなど) です。つまり、PHI は、呼び出しごとに 1 つの行ベクトルでなければなりません。nlmefitsa は、必要に応じて、1 つの PHI ベクトルと、一度に 1 つの観測またはグループに対する行を含む XFUN を使用して MODELFUN をループ内で呼び出します。VFUN は、XFUN のすべての行に適用する 1 つの行、または XFUN 内の行に対応する行をもつ行列です。 'SingleGroup' — MODELFUN は、データ内の 1 つのグループに対応する入力のみを受け入れるため、XFUN は呼び出しごとに 1 つのグループの X の行を含んでいなければなりません。モデルに応じて、PHI は、グループ全体に適用する 1 つの行、または観測ごとに 1 つの行をもつ行列です。VFUN は、1 つの行です。 'Full' — MODELFUN は、データ内の複数のパラメーター ベクトルと複数のグループの入力を受け入れます。PHI または VFUN は、XFUN のすべての行に適用する 1 つの行、または XFUN 内の行に対応する行をもつ行列です。このオプションを使用すると、MODELFUN への呼び出し回数が減ることでパフォーマンスが向上しますが、PHI または V で単項拡張を実行するために MODELFUN が必要です。 'Vectorization' の既定値は 'SinglePhi' です。いずれの場合も、V が空の場合、nlmefitsa は 2 つの入力で MODELFUN を呼び出します。

出力引数

BETA	固定効果の推定。
PSI	変量効果に対する r 行 r 列の推定された共分散行列。既定の設定では、r はモデル パラメーター数 p と等しくなります。
STATS	構造体には次のフィールドがあります。 logl — 近似したモデルに対する最大対数尤度。LogLikMethod パラメーターに 'none' の既定値がある場合は空です。 rmse — 推定誤差分散の平方根 (exponential 誤差モデルの対数スケールで計算) errorparam — 誤差分散モデルの推定パラメーター aic — 赤池情報量基準 (logl が空の場合は空)。aic = –2 * logl + 2 * numParam のように計算します。ここで次のようになります。 logl は、最大対数尤度です。 numParam は、変量効果に対する共分散行列の自由度、固定効果の数、および誤差モデルのパラメーターの数を含む近似パラメーターの数です。 bic — ベイズ情報量基準 (logl が空の場合は空)。bic = -2*logl + log(M) * numParam のように計算します。 M は、グループの数です。 logl と numParam は、aic で定義されます。 一部の文献では、bic の計算を bic = -2*logl + log(N) * numParam としています。N は観測数です。出力の値を調整するために、bic を以下のように再定義できます。 bic = bic - numel(unique(group)) + numel(Y) sebeta — BETA に対する標準誤差 (ComputeStdErrors パラメーターが既定値の false である場合は空)。 covb — パラメーター推定の推定された共分散 (ComputeStdErrors が false の場合は空)。 dfe — 誤差自由度。 pres — 母集団の残差 (y-y_population)。ここで y_population は、母集団の予測値です。 ires — 母集団の残差 (y-y_population)。ここで y_population は、個別の予測値です。 pwres — 母集団の重み付きの残差 cwres — 条件付き重み付きの残差 iwres — 個別の重み付きの残差

例

確率的な EM アルゴリズムを使用した非線形混合効果モデル

標本データを読み込みます。

load indomethacin

8 時間にわたり 6 人の被験者の血流内における薬物インドメタシンの濃度に関するデータにモデルを当てはめます。

model = @(phi,t)(phi(:,1).*exp(-phi(:,2).*t)+phi(:,3).*exp(-phi(:,4).*t));
phi0 = [1 1 1 1];
xform = [0 1 0 1]; % log transform for 2nd and 4th parameters
[beta,PSI,stats,br] = nlmefitsa(time,concentration, ...
   subject,[],model,phi0,'ParamTransform',xform)

beta = 4×1

    0.8630
   -0.7897
    2.7762
    1.0785

PSI = 4×4

    0.0585         0         0         0
         0    0.0248         0         0
         0         0    0.5068         0
         0         0         0    0.0139

stats = struct with fields:
          logl: []
           aic: []
           bic: []
        sebeta: []
           dfe: 57
          covb: []
    errorparam: 0.0811
          rmse: 0.0772
          ires: [66x1 double]
          pres: [66x1 double]
         iwres: [66x1 double]
         pwres: [66x1 double]
         cwres: [66x1 double]

br = 4×6

   -0.2302   -0.0033    0.1625    0.1774   -0.3334    0.1129
    0.0363   -0.1502    0.0071    0.0471    0.0068   -0.0481
   -0.7631   -0.0553    0.8780   -0.8120    0.5429    0.1695
   -0.0030   -0.0223    0.0192   -0.0830    0.0505   -0.0066

データを母集団全体の当てはめに沿ってプロットします。

figure
phi = [beta(1),exp(beta(2)),beta(3),exp(beta(4))]; 
h = gscatter(time,concentration,subject);
xlabel('Time (hours)')
ylabel('Concentration (mcg/ml)')
title('{\bf Indomethacin Elimination}')
xx = linspace(0,8);
line(xx,model(phi,xx),'linewidth',2,'color','k')

変量効果推定値に基づいて個別曲線をプロットします。

for j=1:6
    phir = [beta(1)+br(1,j), exp(beta(2)+br(2,j)), ...
            beta(3)+br(3,j), exp(beta(4)+br(4,j))];
    line(xx,model(phir,xx),'color',get(h(j),'color'))
end

アルゴリズム

非線形混合効果モデルのパラメーターを推定するために、尤度関数を最大にするパラメーター値を選択する必要があります。これらの値は、最尤推定値と呼ばれます。尤度関数は、次の形式で記述できます。

$$p\left(y|\beta ,{\sigma }^2,\sum \right)={\displaystyle \int p\left(y|\beta ,b,{\sigma }^2\right)p\left(b|\sum \right)db},$$

ここで

• y は、応答データです。

• β は、母集団係数のベクトルです。

• σ² は、残差分散です。

• ∑ は、変量効果に対する共分散行列です。

• b は、観測されない変量効果のセットです。

右側の各関数 p() は、共変量因子に応じた正規 (ガウス) 尤度関数です。

積分は閉形式でないため、最大になるパラメーターを見つけることは困難です。Delyon、Lavielle、および Moulines [1] は、E ステップを確率的な手順で置き換える期待値最大化 (EM) アルゴリズムを使用して、最尤推定値を見つけることを提案しています。彼らはこの方法をアルゴリズム SAEM (確率的近似 EM) と名付けました。また、実際条件における収束や尤度関数のローカル最大値への収束など、このアルゴリズムには望ましい理論特性があることを示しました。彼らの提案には、以下の 3 つのステップがあります。

1. シミュレーション: 現在のパラメーター推定値を前提として事後密度 p(b|Σ) から変量効果 b の値をシミュレートして生成します。

2. 確率近似: 前のステップの値を使用し、シミュレートした変量効果から計算した対数尤度の平均値に近づけて、対数尤度関数の期待値を更新します。

3. 最大化ステップ: シミュレートした変量効果の値を使用して対数尤度関数が最大になるように新しいパラメーター推定値を選択します。

参考文献

[1] Delyon, B., M. Lavielle, and E. Moulines, "Convergence of a stochastic approximation version of the EM algorithm." Annals of Statistics, 27, 94-128, 1999.

[2] Mentré, F., and M. Lavielle, "Stochastic EM Algorithms in Population PKPD analyses." American Conference on Pharmacometrics, 2008.

バージョン履歴

R2010a で導入

参考

nlinfit | nlpredci | nlmefit

トピック

• 混合効果のモデル
• グループ化変数