混合効果のモデル

混合効果のモデルの紹介

統計学における "効果" とは、特定の設定の予測子変数で応答変数の値に何らかの影響を与えるものを指します。効果はモデル パラメーターに変換されます。線形モデルでは、効果は係数になり、モデル項の比例の寄与を表します。非線形モデルでは、効果はしばしば特定の物理的な解釈をもち、より一般的な非線形の組み合わせで現れます。

"固定効果" は、データを収集するたびに同じ値になると見なされる母集団のパラメーターを表します。固定効果の推定は、回帰モデリングの従来のドメインです。これに対して、"変量効果" は標本によって異なる確率変数です。モデリングでは、変量効果は追加の誤差項のように機能します。その分布と共分散を指定しなければなりません。

たとえば、血流から薬物を除去するモデルを考えてみましょう。モデルは、予測子として時間 t を使用し、応答として薬物 C の集中を使用します。非線形モデル項 C₀e^–rt は、パラメーター C₀ と r を結合し、それぞれ初期の集中と除去レートを表します。複数の個体についてデータを収集する場合、除去レートは個々の i ごとに異なる確率変数 r_i であり、母集団平均 $$\overline{r}$$ の周辺で変動すると仮定することが妥当です。項 C₀e^–rt は、以下のようになります。

$$C_0{e}^{-[\overline{r}+(r_i-\overline{r})]t}=C_0{e}^{-(\beta +b_i)t}\text{,}$$

ここで、β = $$\overline{r}$$ は固定効果、b_i = $$r_i-\overline{r}$$ は変量効果です。

データが自然なグループになる場合に、変量効果が役立ちます。薬物除去モデルでは、グループは単に研究対象の個体です。より洗練されたモデルでは、個人の年齢、体重、食事などでデータをグループ化することが考えられます。グループは研究の焦点ではありませんが、変量効果をモデルに追加すると、特定の個人の標本を対象にする場合よりも推測の信頼性が高まります。

"混合効果モデル" では、固定効果と変量効果の両方を考慮します。すべての回帰モデルのように、その目的は、応答変数を予測子変数の関数として記述することです。しかし、混合効果のモデルは、標本サブグループ内の相関関係を認識します。この方法で、データ グループの完全な無視と別個のモデルによる各グループの近似の間の妥協案が提供されます。

混合効果のモデルの階層

非線形回帰モデルのデータが m 異なるグループ i = 1, ..., m の 1 つになったとしてみましょう モデルのグループを説明するには、以下のように、グループ i の応答 j を記述します。

$$y_{ij}=f(\phi ,x_{ij})+{\epsilon }_{ij}$$

y_ij は応答です。x_ij は予測子のベクトルです。φ モデル パラメーターのベクトルです。ε_ij は測定または処理誤差です。インデックス j の範囲は、1 ～ n_i です。ここで、n_i は、グループ i の観測数です。関数 f はモデルの形式を指定します。多くの場合、x_ij は単に観測時間 t_ij です。通常、誤差は独立していると仮定され、一定の分散で同一かつ正常に分布します。

その推定がすべてのグループに対して同じであると仮定して、φ におけるパラメーターの推定は母集団を説明します。しかし、推定がグループによって変化する場合は、モデルは以下のようになります。

$$y_{ij}=f({\phi }_i,x_{ij})+{\epsilon }_{ij}$$

混合効果のモデルでは、φ_i は固定効果と変量効果の組み合わせである場合があります。

$${\phi }_i=\beta +b_i$$

変量効果 b_i は、通常、平均 0 と共分散 Ψ により正常に分布された多変量として記述されます。固定効果 β と変量効果 Ψ の共分散を推定することにより、パラメーター φ_i がグループ全体で同じであると仮定しない母集団が記述されます。また、変量効果 b_i を推定すると、データ内の特定のグループも記述されます。

モデル パラメーターを個々の効果で特定する必要はありません。一般に、"計画行列" A および B が固定効果と変量効果の線形結合におけるパラメーターを指定するために使用されます。

$${\phi }_i=A\beta +B{b}_i$$

計画行列がグループ間で異なる場合は、モデルは以下のようになります。

$${\phi }_i=A_i\beta +B_i{b}_i$$

また、計画行列が観測間で異なる場合は、モデルは以下のようになります。

$$\begin{array}{l}{\phi }_{ij}=A_{ij}\beta +B_{ij}{b}_i\\ y_{ij}=f({\phi }_{ij},x_{ij})+{\epsilon }_{ij}\end{array}$$

x_ij 内のグループ固有の予測子によっては、観測 j と交換されないものがあります。それらの v_i を呼び出すと、モデルは以下のようになります。

$$y_{ij}=f({\phi }_{ij},x_{ij},v_i)+{\epsilon }_{ij}$$

混合効果のモデルを指定

非線形回帰モデルのデータが m 異なるグループ i = 1, ..., m の 1 つになったとしてみましょう (特に、グループが入れ子にされていないと仮定します)。このデータに対する一般的な非線形混合効果モデルを指定するには、以下のようにします。

1. グループ固有のモデル パラメーター φ_i を固定効果 β と変量効果 b_i の線形結合と定義します。

2. 応答値 y_i をパラメーターとグループ固有の予測子変数 X_i の非線形の関数 f として定義します。

モデルは以下のようになります。

$$\begin{array}{l}{\phi }_i=A_i\beta +B_i{b}_i\\ y_i=f({\phi }_i,X_i)+{\epsilon }_i\\ b_i\sim N(0,\Psi )\\ {\epsilon }_i\sim N(0,{\sigma }^2)\end{array}$$

この非線形混合効果モデルの定式化は以下の符号を使用します。

φi	グループ固有のモデル パラメーターのベクトル
β	固定効果のベクトル、母集団のパラメーターのモデリング
bi	通常、グループ固有の変量効果が分布される多変量のベクトル
Ai	固定効果を結合するためのグループ固有の計画行列
Bi	変量効果を結合するためのグループ固有の計画行列
Xi	グループ固有の予測子値のデータ行列
yi	グループ固有の応答値のデータ ベクトル
f	φi と Xi の一般的実数値の関数
εi	bi とは別に同一かつ正常に分布され、独立していると仮定されるグループ固有の誤差のベクトル
Ψ	変量効果に対する共分散行列
σ2	観測全体に一定であると仮定される誤差分散

たとえば、血流から薬物を除去するモデルを考えてみましょう。モデルは以下の 2 つの重なり段階を含んでいます。

• 薬物の血中濃度が周辺の組織と平衡になる初期段階 p

• 薬物が血流から除去される 2 番目の段階 q

複数の個体 i のデータの場合、モデルは以下のようになります。

$$y_{ij}=C_{pi}{e}^{-r_{pi}{t}_{ij}}+C_{qi}{e}^{-r_{qi}{t}_{ij}}+{\epsilon }_{ij},$$

ここで、y_ij は、t_ij 時の、個々の i の観測された集中です。モデルでは、個人によって異なるサンプリング時間と異なる観測数が許容されます。

除去レートの r_pi および r_qi は、物理的に有意義である陽性でなければなりません。log レートの R_pi = log(r_pi) および R_qi = log(r_qi) を導入し、モデルを再パラメーターすることにより、これを強制します。

$$y_{ij}=C_{pi}{e}^{-\exp (R_{pi})t_{ij}}+C_{qi}{e}^{-\exp (R_{qi})t_{ij}}+{\epsilon }_{ij}$$

変量効果でモデル化するパラメーターを選択することは、混合効果のモデルを構成する際に重要です。その手法の 1 つは、変量効果をすべてのパラメーターに追加し、それらの分散の推定を使用して、モデル内でその意義を決定することです。もう 1 つの方法は、変量効果なしにモデルを別々に各グループに当てはめ、パラメーター推定の変化を見ることです。推定がグループ全体に多様であるか、グループごとに信頼区間が最小の重なりをもっている場合は、パラメーターは変量効果に対して良い候補です。

すべてのモデル パラメーターのために固定効果 β と変量効果 b_i を導入するには、以下のモデルを再表現します。

$$\begin{array}{lll}{y}_{ij}\hfill & =\hfill & [{\overline{C}}_p+(C_{pi}-{\overline{C}}_p)]e^{-\exp [{\overline{R}}_p+(R_{pi}-{\overline{R}}_p)]t_{ij}}+\hfill \\ \hfill & \hfill & [{\overline{C}}_q+(C_{qi}-{\overline{C}}_q)]e^{-\exp [{\overline{R}}_q+(R_{qi}-{\overline{R}}_q)]t_{ij}}+{\epsilon }_{ij}\hfill \\ \hfill & =\hfill & ({\beta }_1+b_{1i})e^{-\exp ({\beta }_2+b_{2i})t_{ij}}+\hfill \\ \hfill & \hfill & ({\beta }_3+b_{3i})e^{-\exp ({\beta }_4+b_{4i})t_{ij}}+{\epsilon }_{ij}\hfill \end{array}$$

一般的モデルの表記法で、以下のようにします。

$$\beta =\left(\begin{array}{c}{\beta }_1\\ ⋮\\ {\beta }_4\end{array}\right),b_i=\left(\begin{array}{c}{b}_{i1}\\ ⋮\\ b_{i4}\end{array}\right),y_i=\left(\begin{array}{c}{y}_{i1}\\ ⋮\\ y_{i{n}_i}\end{array}\right),X_i=\left(\begin{array}{c}{t}_{i1}\\ ⋮\\ t_{i{n}_i}\end{array}\right),$$

ここで、n_i は、個体 i の観測数です。この例では、計画行列 A_i と B_i は、少なくとも最初は 4 行 4 列の単位行列です。必要に応じて、計画行列を変更し、個々の効果または時間従属性の重みを導入することができます。

モデルを近似し、共分散行列 Ψ を推定することにより、多くの場合さらに改善できます。変量効果の分散に対する比較的小さい推定は、モデルから除去できることを示します。同様に、一定の変量効果の間の共分散に対する比較的小さい推定は、完全な共分散行列が不要であることを示します。変量効果が観測されないため、Ψ を間接的に推定しなければなりません。Ψ に対する対角またはブロック対角共分散パターンを指定することにより、近似アルゴリズムの収束と効率を改善することができます。

Statistics and Machine Learning Toolbox™ の関数 nlmefit および nlmefitsa は、固定効果と変量効果を推定して、データに一般的非線形混合効果モデルを当てはめます。また、関数は変量効果に対する共分散行列 Ψ も推定します。追加の診断出力により、モデル パラメーターの数と適合度の間のトレードオフを評価することができます。

共変量モデルの指定

混合効果のモデルを指定のモデルで体重 (w) などのようにグループに依存する共変量を仮定すると、モデルは次のようになります。

$$\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\\ {\phi }_3\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& w_i\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right)\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\\ \begin{array}{l}{\beta }_3\\ {\beta }_4\end{array}\end{array}\right)+\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right)$$

したがって、i 番目のグループではすべての個体についてパラメーター φ_i が次のようになります。

$$\left(\begin{array}{c}{\phi }_{1_i}\\ {\phi }_{2_i}\\ {\phi }_{3_i}\end{array}\right)=\left(\begin{array}{c}{\beta }_1+{\beta }_4{w}_i\\ {\beta }_2\\ {\beta }_3\end{array}\right)+\left(\begin{array}{c}{b}_{1_i}\\ b_{2_i}\\ b_{3_i}\end{array}\right)$$

共変量モデルを指定するには、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)[1]という別の近似法を使用できます。

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

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

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

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 フィールドは、反復ごとにソルバーが呼び出す関数を指定します。通常、出力関数は各反復で点をプロットしたり、アルゴリズムからの最適化量を表示するために使用します。出力関数を設定するには、以下を行います。

1. 出力関数を MATLAB^® ファイル関数またはローカル関数として記述します。

2. statset を使用して Outputfcn の値を関数ハンドル、すなわち記号 @ が前に付いた関数名を設定します。たとえば、出力関数が outfun.m の場合、コマンド

    options = statset('OutputFcn', @outfun);

   は、OutputFcn を outfun へのハンドルになるように指定します。複数の出力関数を指定するには、次の構文を使用します。

    options = statset('OutputFcn',{@outfun, @outfun2});

3. 入力引数に options を使用して最適化関数を呼び出します。

出力関数の例については、出力関数の例を参照してください。

出力関数の構造

出力関数の関数定義行は、次の形式になります。

stop = outfun(beta,status,state)

ここで

• beta は現在の固定効果です。

• status は、現在の反復から得られるデータを含む構造体です。この構造体についての詳細は、status のフィールドを参照してください。

• state は、アルゴリズムの現在の状態です。値は、アルゴリズムの状態に示されています。

• stop は、最適化ルーチンを終了するか継続するかによって、true または false を表すフラグです。詳細は、stop フラグを参照してください。

ソルバーは、反復処理ごとに入力引数の値を outfun に渡します。

status のフィールド

次の表は status 構造体のフィールドを示します。

フィールド	説明
procedure	'ALT' ― 線形混合効果または制限された線形混合効果の近似の最適化用代替アルゴリズム 'LAP' ― 1 次または 1 次条件付き推定に対するラプラシアン近似の最適化
iteration	0 から始まる整数
inner	ALT と LAP 手順内で内部反復の状況を説明する構造体。次の各フィールドがあります。 procedure ― procedure が 'ALT' である場合 'PNLS' (ペナルティが課された非線形最小二乗法) 'LME' (線形混合効果の推定) 'none' procedure が 'LAP' である場合 'PNLS' (ペナルティが課された非線形最小二乗法) 'PLM' (プロファイルされた尤度最大化) 'none' state ― 以下のいずれかです。 'init' 'iter' 'done' 'none' iteration ― 0 から始まる整数、またはNaN。バーンイン反復のある nlmefitsa では、出力関数は STATUS.iteration で負の値をもつ各反復の後で呼び出されます。
fval	現在の対数尤度
Psi	現在の変量効果の共分散行列
theta	Psi の現在のパラメーター表現
mse	現在の誤差分散

アルゴリズムの状態

次の表は、state の値を示しています。

状態	説明
'init'	アルゴリズムは、最初の反復の前の初期状態にあります。
'iter'	アルゴリズムは、反復の最後にあります。
'done'	最後の反復後、アルゴリズムは最終状態になっています。

以下のコードは現在の反復で、どのタスクを実行するかを決めるために出力関数が 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, …)

生成された Figure を閉じると、nlmefitoutputfcn は nlmefit または nlmefitsa を停止します。

参考

nlmefit | nlmefitsa

トピック

• nlmefit と nlmefitsa を使用する混合効果のモデル
• モデルの検証のための残差の検査

参照

[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.