混合効果のモデル
混合効果のモデルの紹介
統計学における "効果" とは、特定の設定の予測子変数で応答変数の値に何らかの影響を与えるものを指します。効果はモデル パラメーターに変換されます。線形モデルでは、効果は係数になり、モデル項の比例の寄与を表します。非線形モデルでは、効果はしばしば特定の物理的な解釈をもち、より一般的な非線形の組み合わせで現れます。
"固定効果" は、データを収集するたびに同じ値になると見なされる母集団のパラメーターを表します。固定効果の推定は、回帰モデリングの従来のドメインです。これに対して、"変量効果" は標本によって異なる確率変数です。モデリングでは、変量効果は追加の誤差項のように機能します。その分布と共分散を指定しなければなりません。
たとえば、血流から薬物を除去するモデルを考えてみましょう。モデルは、予測子として時間 t を使用し、応答として薬物 C の集中を使用します。非線形モデル項 C0e–rt は、パラメーター C0 と r を結合し、それぞれ初期の集中と除去レートを表します。複数の個体についてデータを収集する場合、除去レートは個々の i ごとに異なる確率変数 ri であり、母集団平均 の周辺で変動すると仮定することが妥当です。項 C0e–rt は、以下のようになります。
ここで、β = は固定効果、bi = は変量効果です。
データが自然なグループになる場合に、変量効果が役立ちます。薬物除去モデルでは、グループは単に研究対象の個体です。より洗練されたモデルでは、個人の年齢、体重、食事などでデータをグループ化することが考えられます。グループは研究の焦点ではありませんが、変量効果をモデルに追加すると、特定の個人の標本を対象にする場合よりも推測の信頼性が高まります。
"混合効果モデル" では、固定効果と変量効果の両方を考慮します。すべての回帰モデルのように、その目的は、応答変数を予測子変数の関数として記述することです。しかし、混合効果のモデルは、標本サブグループ内の相関関係を認識します。この方法で、データ グループの完全な無視と別個のモデルによる各グループの近似の間の妥協案が提供されます。
混合効果のモデルの階層
非線形回帰モデルのデータが m 異なるグループ i = 1, ..., m の 1 つになったとしてみましょう モデルのグループを説明するには、以下のように、グループ i の応答 j を記述します。
yij は応答です。xij は予測子のベクトルです。φ モデル パラメーターのベクトルです。εij は測定または処理誤差です。インデックス j の範囲は、1 ~ ni です。ここで、ni は、グループ i の観測数です。関数 f はモデルの形式を指定します。多くの場合、xij は単に観測時間 tij です。通常、誤差は独立していると仮定され、一定の分散で同一かつ正常に分布します。
その推定がすべてのグループに対して同じであると仮定して、φ におけるパラメーターの推定は母集団を説明します。しかし、推定がグループによって変化する場合は、モデルは以下のようになります。
混合効果のモデルでは、φi は固定効果と変量効果の組み合わせである場合があります。
変量効果 bi は、通常、平均 0 と共分散 Ψ により正常に分布された多変量として記述されます。固定効果 β と変量効果 Ψ の共分散を推定することにより、パラメーター φi がグループ全体で同じであると仮定しない母集団が記述されます。また、変量効果 bi を推定すると、データ内の特定のグループも記述されます。
モデル パラメーターを個々の効果で特定する必要はありません。一般に、"計画行列" A および B が固定効果と変量効果の線形結合におけるパラメーターを指定するために使用されます。
計画行列がグループ間で異なる場合は、モデルは以下のようになります。
また、計画行列が観測間で異なる場合は、モデルは以下のようになります。
xij 内のグループ固有の予測子によっては、観測 j と交換されないものがあります。それらの vi を呼び出すと、モデルは以下のようになります。
混合効果のモデルを指定
非線形回帰モデルのデータが m 異なるグループ i = 1, ..., m の 1 つになったとしてみましょう (特に、グループが入れ子にされていないと仮定します)。このデータに対する一般的な非線形混合効果モデルを指定するには、以下のようにします。
グループ固有のモデル パラメーター φi を固定効果 β と変量効果 bi の線形結合と定義します。
応答値 yi をパラメーターとグループ固有の予測子変数 Xi の非線形の関数 f として定義します。
モデルは以下のようになります。
この非線形混合効果モデルの定式化は以下の符号を使用します。
φi | グループ固有のモデル パラメーターのベクトル |
β | 固定効果のベクトル、母集団のパラメーターのモデリング |
bi | 通常、グループ固有の変量効果が分布される多変量のベクトル |
Ai | 固定効果を結合するためのグループ固有の計画行列 |
Bi | 変量効果を結合するためのグループ固有の計画行列 |
Xi | グループ固有の予測子値のデータ行列 |
yi | グループ固有の応答値のデータ ベクトル |
f | φi と Xi の一般的実数値の関数 |
εi | bi とは別に同一かつ正常に分布され、独立していると仮定されるグループ固有の誤差のベクトル |
Ψ | 変量効果に対する共分散行列 |
σ2 | 観測全体に一定であると仮定される誤差分散 |
たとえば、血流から薬物を除去するモデルを考えてみましょう。モデルは以下の 2 つの重なり段階を含んでいます。
薬物の血中濃度が周辺の組織と平衡になる初期段階 p
薬物が血流から除去される 2 番目の段階 q
複数の個体 i のデータの場合、モデルは以下のようになります。
ここで、yij は、tij 時の、個々の i の観測された集中です。モデルでは、個人によって異なるサンプリング時間と異なる観測数が許容されます。
除去レートの rpi および rqi は、物理的に有意義である陽性でなければなりません。log レートの Rpi = log(rpi) および Rqi = log(rqi) を導入し、モデルを再パラメーターすることにより、これを強制します。
変量効果でモデル化するパラメーターを選択することは、混合効果のモデルを構成する際に重要です。その手法の 1 つは、変量効果をすべてのパラメーターに追加し、それらの分散の推定を使用して、モデル内でその意義を決定することです。もう 1 つの方法は、変量効果なしにモデルを別々に各グループに当てはめ、パラメーター推定の変化を見ることです。推定がグループ全体に多様であるか、グループごとに信頼区間が最小の重なりをもっている場合は、パラメーターは変量効果に対して良い候補です。
すべてのモデル パラメーターのために固定効果 β と変量効果 bi を導入するには、以下のモデルを再表現します。
一般的モデルの表記法で、以下のようにします。
ここで、ni は、個体 i の観測数です。この例では、計画行列 Ai と Bi は、少なくとも最初は 4 行 4 列の単位行列です。必要に応じて、計画行列を変更し、個々の効果または時間従属性の重みを導入することができます。
モデルを近似し、共分散行列 Ψ を推定することにより、多くの場合さらに改善できます。変量効果の分散に対する比較的小さい推定は、モデルから除去できることを示します。同様に、一定の変量効果の間の共分散に対する比較的小さい推定は、フルの共分散行列が不要であることを示します。変量効果が観測されないため、Ψ を間接的に推定しなければなりません。Ψ に対する対角またはブロック対角共分散パターンを指定することにより、近似アルゴリズムの収束と効率を改善することができます。
Statistics and Machine Learning Toolbox™ の関数 nlmefit
および nlmefitsa
は、固定効果と変量効果を推定して、データに一般的非線形混合効果モデルを当てはめます。また、関数は変量効果に対する共分散行列 Ψ も推定します。追加の診断出力により、モデル パラメーターの数と適合度の間のトレードオフを評価することができます。
共変量モデルの指定
混合効果のモデルを指定のモデルで体重 (w) などのようにグループに依存する共変量を仮定すると、モデルは次のようになります。
したがって、i
番目のグループではすべての個体についてパラメーター φi が次のようになります。
共変量モデルを指定するには、nlmefit
または nlmefitsa
の名前と値の引数 'FEGroupDesign'
を使用します。
'FEGroupDesign'
は、m
グループごとに異なる p-by-q
の固定効果の計画行列を指定する p-by-q-by-m
の配列です。上記の例を使用する配列は以下のようになります。
% Number of parameters in the model (Phi) num_params = 3; % Number of covariates num_cov = 1; % Assuming number of groups in the data set is 7 num_groups = 7; % Array of covariate values covariates = [75; 52; 66; 55; 70; 58; 62 ]; A = repmat(eye(num_params, num_params+num_cov),... [1,1,num_groups]); A(1,num_params+1,1:num_groups) = covariates(:,1)
計画行列 A
を使用するため、'FEGroupDesign',A
を指定します。
nlmefit
または nlmefitsa
の選択
Statistics and Machine Learning Toolbox には、非線形混合効果モデルの近似を行う nlmefit
と nlmefitsa
の 2 つの関数があります。各関数には個別の機能があり、その機能に応じて使用する関数を決めます。
近似法
nlmefit
には、非線形混合効果モデルの近似を行う、以下の 4 つの近似法があります。
'LME'
—beta
とB
の現在の条件推定で線形混合効果モデルに対する尤度を使用します。これは既定の設定です。'RELME'
—beta
とB
の現在の条件推定で線形混合効果モデルに対する制限付き尤度を使用します。'FO'
— 変量効果のない 1 次ラプラシアン近似。'FOCE'
—B
の条件推定での 1 次ラプラシアン近似。
nlmefitsa
では、次の 3 つのステップがある確率近似期待値最大化 (Stochastic Approximation Expectation-Maximization: SAEM)[25]という別の近似法を使用できます。
シミュレーション: 現在のパラメーター推定値を前提として事後密度 p(b|Σ) から変量効果 b の値をシミュレートして生成します。
確率近似: 前のステップの値を使用し、シミュレートした変量効果から計算した対数尤度の平均値に近づけて、対数尤度関数の期待値を更新します。
最大化ステップ: シミュレートした変量効果の値を使用して対数尤度関数が最大になるように新しいパラメーター推定値を選択します。
nlmefit
と nlmefitsa
はどちらも尤度関数を最大化するパラメーター推定を見つけようとしますが、この計算は困難です。nlmefit
は、尤度関数をさまざまな方法で近似し、近似関数を最大化することにより、この問題に対応します。この関数は、収束基準や反復制限などを使用する従来の最適化手法を使用します。
一方、nlmefitsa
は、正確な尤度関数を最大化する値に結果的に収束する方法で、パラメーターのランダム値をシミュレートします。結果はランダムであり、従来の収束テストは適用されません。そのため、nlmefitsa
には、シミュレーションが進行すると結果をプロットし、シミュレーションを複数回再開始するオプションがあります。これらの機能を使用して、必要な精度に結果が収束したかどうかを判断できます。
nlmefitsa
専用のパラメーター
以下のパラメーターは nlmefitsa
専用です。ほとんどのパラメーターは確率アルゴリズムを制御します。
Cov0
― 共分散行列PSI
の初期値です。r 行 r 列の正定値行列でなければなりません。空の場合、既定値はBETA0
の値に依存します。ComputeStdErrors
― 係数の推定の標準誤差を計算し、出力のSTATS
構造体に格納する場合はtrue
、この計算を省略する場合はfalse
(既定値) です。LogLikMethod
― 対数尤度を近似する方法を指定します。NBurnIn
― パラメーター推定が再計算されない、最初のバーンイン反復の数です。既定の設定は 5 です。NIterations
— このアルゴリズムの 3 つの段階それぞれについて反復回数を制御します。NMCMCIterations
— マルコフ連鎖モンテカルロ (MCMC) 反復回数です。
モデルおよびデータの要件
nlmefit
と nlmefitsa
の機能にはいくつかの相違点があります。そのため、どちらの関数でも使用できるデータおよびモデルもありますが、一部のデータおよびモデルでは関数を選択する必要がある場合があります。
"誤差モデル" ―
nlmefitsa
はさまざまな誤差モデルをサポートします。たとえば、応答の標準偏差は一定で、関数の値または 2 つの値の組み合わせに比例します。nlmefit
は、応答の標準偏差が一定であるという仮定の下、モデルを近似します。誤差モデルの 1 つ、'exponential'
は、応答の対数に一定の標準偏差があると指定します。入力として対数の応答を指定したり、モデル関数を再記述して非線形関数値の対数を生成したりすることによって、nlmefit
を使用してこのようなモデルを当てはめることができます。"変量効果" ― どちらの関数もデータをパラメーターをもつ非線形関数に近似し、そのパラメーターが単純なスカラー値または共変量の線形関数である可能性があります。
nlmefit
を使用すると、線形関数の係数は固定効果と変量効果の両方をもつことができます。nlmefitsa
は、線形関数の定数 (切片) 係数のみに対する変量効果をサポートし、勾配係数に対する変量効果はサポートしません。このため、共変量モデルの指定の例では、nlmefitsa
は最初の 3 つのベータ値のみを変量効果として扱います。"モデル形式" ―
nlmefit
は、固定係数とモデル パラメーターへの変量効果を関連付ける計画行列に対する制限がほとんどない、一般的なモデルの設定をサポートします。nlmefitsa
には、さらに以下のような制限があります。固定効果の設計は、各グループ (各個体) で一定でなければならないため、観測に依存する設計はサポートされません。
変量効果の設計はデータ セット全体で一定でなければならないため、観測に依存する設計もグループに依存する設計もサポートされません。
"変量効果" で説明したとおり、変量効果の設計で勾配係数の変量効果を指定してはなりません。つまり、設計が 0 と 1 で構成されなければなりません。
変量効果の設計は複数の係数に対して同じランダムな影響を使用してはなりません。また、1 つの係数に対して複数の変量効果を使用できません。
固定効果の設計は、複数のパラメーターに対して同じ係数を使用してはなりません。つまり、各列に 0 以外の値を最大 1 つしかもてません。
共変量効果が無作為であるデータに
nlmefitsa
を使用するには、共変量を非線形モデル式に直接含めます。固定効果または変量効果の計画行列に共変量を含めないでください。"収束" ― 「モデル形式」で説明したとおり、
nlmefit
とnlmefitsa
は異なる手法で収束を測定します。nlmefit
は従来の最適化測定を使用し、nlmefitsa
はランダムなシミュレーションの収束の判定に役立つ診断を行います。
実際には、nlmefitsa
はロバストで難しい問題でも失敗することはあまりありません。ただし、収束する場合は nlmefit
は問題に対しすばやく収束できます。問題によっては、これらの手法を組み合わせると有効な場合もあります。たとえば、nlmefitsa
をしばらく実行して適当なパラメーター推定を取得し、nlmefit
を使用した次の反復の開始点としてそれらを使用できます。
混合効果のモデルによる出力関数の使い方
options
構造体の Outputfcn
フィールドは、反復ごとにソルバーが呼び出す関数を指定します。通常、出力関数は各反復で点をプロットしたり、アルゴリズムからの最適化量を表示するために使用します。出力関数を設定するには、以下を行います。
出力関数を MATLAB® ファイル関数またはローカル関数として記述します。
statset
を使用してOutputfcn
の値を関数ハンドル、すなわち記号@
が前に付いた関数名を設定します。たとえば、出力関数がoutfun.m
の場合、コマンドoptions = statset('OutputFcn', @outfun);
は、
OutputFcn
をoutfun
へのハンドルになるように指定します。複数の出力関数を指定するには、次の構文を使用します。options = statset('OutputFcn',{@outfun, @outfun2});
入力引数に
options
を使用して最適化関数を呼び出します。
出力関数の例については、出力関数の例を参照してください。
出力関数の構造
出力関数の関数定義行は、次の形式になります。
stop = outfun(beta,status,state)
beta は現在の固定効果です。
status は、現在の反復から得られるデータを含む構造体です。この構造体についての詳細は、status のフィールドを参照してください。
state は、アルゴリズムの現在の状態です。値は、アルゴリズムの状態に示されています。
stop は、最適化ルーチンを終了するか継続するかによって、
true
またはfalse
を表すフラグです。詳細は、stop フラグを参照してください。
ソルバーは、反復処理ごとに入力引数の値を outfun
に渡します。
status のフィールド
次の表は status
構造体のフィールドを示します。
フィールド | 説明 |
---|---|
procedure |
|
iteration | 0 から始まる整数 |
inner |
|
fval | 現在の対数尤度 |
Psi | 現在の変量効果の共分散行列 |
theta | Psi の現在のパラメーター表現 |
mse | 現在の誤差分散 |
アルゴリズムの状態
次の表は、state
の値を示しています。
状態 | 説明 |
---|---|
| アルゴリズムは、最初の反復の前の初期状態にあります。 |
| アルゴリズムは、反復の最後にあります。 |
| 最後の反復後、アルゴリズムは最終状態になっています。 |
以下のコードは現在の反復で、どのタスクを実行するかを決めるために出力関数が state
の値を使用する方法を示します。
switch state case 'iter' % Make updates to plot or guis as needed case 'init' % Setup for plots or guis case 'done' % Cleanup of plots, guis, or final plot otherwise end
stop フラグ
出力引数 stop
は、true
またはfalse
を示すフラグです。このフラグは、停止または連続しなければならないかどうかをソルバーに示します。次の例は、stop
フラグの一般的な使用方法を示したものです。
中間結果に基づき最適化を停止. 出力関数は、渡された引数の値に基づいて、任意の反復で推定を停止できます。たとえば、以下のコード セットは、ステータス構造体の 'fval'
フィールドに格納された対数尤度の値に基づいて、true で停止します。
stop = outfun(beta,status,state) stop = false; % Check if loglikelihood is more than 132. if status.fval > -132 stop = true; end
GUI の入力に基づき反復を停止. nlmefit
の反復を実行する GUI を設計する場合、GUI の [停止] ボタンをクリックすると出力関数が停止するように設計できます。たとえば、次のコードは計算をキャンセルするためのダイアログを実装します。
function retval = stop_outfcn(beta,str,status) persistent h stop; if isequal(str.inner.state,'none') switch(status) case 'init' % Initialize dialog stop = false; h = msgbox('Press STOP to cancel calculations.',... 'NLMEFIT: Iteration 0 '); button = findobj(h,'type','uicontrol'); set(button,'String','STOP','Callback',@stopper) pos = get(h,'Position'); pos(3) = 1.1 * pos(3); set(h,'Position',pos) drawnow case 'iter' % Display iteration number in the dialog title set(h,'Name',sprintf('NLMEFIT: Iteration %d',... str.iteration)) drawnow; case 'done' % Delete dialog delete(h); end end if stop % Stop if the dialog button has been pressed delete(h) end retval = stop; function stopper(varargin) % Set flag to stop when button is pressed stop = true; disp('Calculation stopped.') end end
出力関数の例
nmlefitoutputfcn
は、nlmefit
および nlmefitsa
の Statistics and Machine Learning Toolbox 出力関数の例です。これにより、固定効果 (BETA
) および変量効果 (diag
(STATUS.Psi
)) の分散をもつプロットを初期化および更新します。nlmefit
では、プロットには対数尤度 (STATUS.fval
) も含まれています。
nlmefitoutputfcn
は nlmefitsa
の既定の設定の出力関数です。これを nlmefit
で使用するには、オプション構造体でその関数ハンドルを指定します。
opt = statset('OutputFcn', @nlmefitoutputfcn, …) beta = nlmefit(…, 'Options', opt, …)
nlmefitsa
がこの関数を使用しないようにするには、出力関数に空の値を指定します。opt = statset('OutputFcn', [], …) beta = nlmefitsa(…, 'Options', opt, …)
nlmefitoutputfcn
は nlmefit
または nlmefitsa
を停止します。