ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

混合効果のモデル

混合効果のモデルの紹介

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

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

たとえば、血流から薬物を除去するモデルを考えてみましょう。モデルは、予測子として時間 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 つになったとしてみましょう (特に、グループが入れ子にされていないと仮定します)。このデータに対する一般的な非線形混合効果モデルを指定するには、以下のようにします。

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

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

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

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

φiグループ固有のモデル パラメーターのベクトル
β固定効果のベクトル、母集団のパラメーターのモデリング
bi通常、グループ固有の変量効果が分布される多変量のベクトル
Ai固定効果を結合するためのグループ固有の計画行列
Bi変量効果を結合するためのグループ固有の計画行列
Xiグループ固有の予測子値のデータ行列
yiグループ固有の応答値のデータ ベクトル
fφi と Xi の一般的実数値の関数
εibi とは別に同一かつ正常に分布され、独立していると仮定されるグループ固有の誤差のベクトル
Ψ変量効果に対する共分散行列
σ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 Toolbox™ の関数 nlmefit および nlmefitsa は、固定効果と変量効果を推定して、データに一般的非線形混合効果モデルを近似します。また、関数は変量効果に対する共分散行列 Ψ も推定します。追加の診断出力により、モデル パラメーターの数と適合度の間のトレードオフを評価することができます。

共分散モデルの指定

混合効果のモデルを指定 のモデルが、重み (w) などのグループ依存の共分散を仮定する場合、モデルは次のようになります。

そのため、i 番目のグループ内の任意の固体に対する φi パラメーターは

になります。共分散モデルを指定するには、'FEGroupDesign' オプションを使用します。

'FEGroupDesign' は、m グループごとに異なる p-by-q の固定効果の計画行列を指定する p-by-q-by-m の配列です。上記の例を使用する配列は以下のようになります。

  1. 配列を作成します。

    % 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)
  2. 指定の計画行列で struct を作成します。

    options.FEGroupDesign = A; 
    
  3. nlmefit と nlmefitsa を使用した混合効果のモデル」に示すように、nlmefit (または nlmefitsa) の引数を指定します。

nlmefit または nlmefitsa の選択

Statistics Toolbox には、非線形混合効果モデルの近似を行う nlmefitnlmefitsa の 2 つの関数があります。各関数には個別の機能があり、その機能に応じて使用する関数を決めます。

近似法

nlmefit には、非線形混合効果モデルの近似を行う、以下の 4 つの近似法があります。

  • 'LME'beta および B の現在の条件付き推定での、線形混合効果モデルの尤度を使用します。これは既定の設定です。

  • 'RELME'beta および B の現在の条件付き推定での、線形混合効果モデルの制限付き尤度を使用します。

  • 'FO' ― 変量効果のない 1 次ラプラシアン近似です。

  • 'FOCE'B の条件付き推定での 1 次ラプラシアン近似です。

nlmefitsa には、これ以外にも、以下の 3 つのステップを行う確率近似期待値最大化 (Stochastic Approximation Expectation-Maximization (SAEM) [24] という近似法があります。

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

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

  3. 最大化ステップ: 新しいパラメーター推定を選択し、変量効果のシミュレーションから得られる値を指定して対数尤度関数を最大化します。

nlmefitnlmefitsa はどちらも尤度関数を最大化するパラメーター推定を見つけようとしますが、この計算は困難です。nlmefit は、尤度関数をさまざまな方法で近似し、近似関数を最大化することにより、この問題に対応します。この関数は、収束基準や反復制限などを使用する従来の最適化手法を使用します。

一方、nlmefitsa は、正確な尤度関数を最大化する値に結果的に収束する方法で、パラメーターのランダム値をシミュレートします。結果はランダムであり、従来の収束テストは適用されません。そのため、nlmefitsa には、シミュレーションが進行すると結果をプロットし、シミュレーションを複数回再開始するオプションがあります。これらの機能を使用して、必要な精度に結果が収束したかどうかを判断できます。

nlmefitsa 専用のパラメーター

以下のパラメーターは nlmefitsa 専用です。ほとんどのパラメーターは確率アルゴリズムを制御します。

  • Cov0 ― 共分散行列 PSI の初期値です。r 行 r 列の正定行列でなければなりません。空の場合、既定値は BETA0 の値によって決まります。

  • ComputeStdErrors ― 係数の推定の標準誤差を計算し、出力の STATS 構造体に格納する場合は true、この計算を省略する場合は false (既定値) です。

  • LogLikMethod ― 対数尤度を近似する方法を指定します。

  • NBurnIn ― パラメーター推定が再計算されない、最初のバーンイン反復の数です。既定値は 5 です。

  • NIterations ― アルゴリズムの 3 つの各段階の反復回数を制御します。

  • NMCMCIterations — マルコフ連鎖モンテカルロ (MCMC) 反復回数です。

モデルおよびデータの要件

nlmefitnlmefitsa の機能にはいくつかの相違点があります。そのため、どちらの関数でも使用できるデータおよびモデルもありますが、一部のデータおよびモデルでは関数を選択する必要がある場合があります。

  • "誤差モデル"nlmefitsa はさまざまな誤差モデルをサポートします。たとえば、応答の標準偏差は一定で、関数の値または 2 つの値の組み合わせに比例します。nlmefit は、応答の標準偏差が一定であるという仮定の下、モデルを近似します。誤差モデルの 1 つ、'exponential' は、応答の対数に一定の標準偏差があると指定します。入力として対数の応答を指定したり、モデル関数を再記述して非線形関数値の対数を生成したりすることによって、nlmefit を使用してこのようなモデルを近似できます。

  • "変量効果" ― どちらの関数もデータをパラメーターをもつ非線形関数に近似し、そのパラメーターが単純なスカラー値または共分散の線形関数である可能性があります。nlmefit を使用すると、線形関数の係数は固定効果と変量効果の両方をもつことができます。nlmefitsa は、線形関数の定数 (切片) 係数のみに対する変量効果をサポートし、勾配係数に対する変量効果はサポートしません。このため、「共分散モデルの指定」の例では、nlmefitsa は最初の 3 つのベータ値のみを変量効果として扱います。

  • "モデル形式"nlmefit は、固定係数とモデル パラメーターへの変量効果を関連付ける計画行列に対する制限がほとんどない、一般的なモデルの設定をサポートします。nlmefitsa には、さらに以下のような制限があります。

    • 固定効果の設計は、各グループ (各個体) で一定でなければならないため、観測に依存する設計はサポートされません。

    • 変量効果の設計はデータセット全体で一定でなければならないため、観測に依存する設計もグループに依存する設計もサポートされません。

    • "変量効果" で説明したとおり、変量効果の設計で勾配係数の変量効果を指定することはできません。つまり、設計が 0 と 1 で構成されなければなりません。

    • 変量効果の設計は複数の係数に対して同じランダムな影響を使用できません。また、1 つの係数に対して複数の変量効果を使用できません。

    • 固定効果の設計は、複数のパラメーターに対して同じ係数を使用できません。つまり、各列に 0 以外の値を最大 1 つしかもてません。

    共分散効果がランダムであるデータに nlmefitsa を使用するには、共分散を非線形モデル式に直接含めます。固定効果または変量効果の計画行列に共分散を含めないでください。

  • "収束"「モデル形式」で説明したとおり、nlmefitnlmefitsa は異なる手法で収束を測定します。nlmefit は従来の最適化測定を使用し、nlmefitsa はランダムなシミュレーションの収束の判定に役立つ診断を行います。

実際には、nlmefitsa はロバストで難しい問題でも失敗することはあまりありません。ただし、収束する場合は nlmefit は問題に対しすばやく収束できます。問題によっては、これらの手法を組み合わせると有効な場合もあります。たとえば、nlmefitsa をしばらく実行して適当なパラメーター推定を取得し、nlmefit を使用した次の反復の起点としてそれらを使用できます。

混合効果のモデルによる出力関数の使い方

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

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

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

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

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

     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 次条件付き推定に対するラプラシアン近似の最適化

iteration0 から始まる整数
innerフィールドで ALTLAP 手順内で内部反復の状況を説明する構造体
  • procedureprocedure'ALT' である場合

    • 'PNLS' (ペナルティが課された非線形最小二乗法)

    • 'LME' (線形混合効果の推定)

    • 'none'

    procedure'LAP' である場合、

    • 'PNLS' (ペナルティが課された非線形最小二乗法)

    • 'PLM' (プロファイルされた尤度最大化)

    • 'none'

  • state ― 以下のいずれかです。

    • 'init'

    • 'iter'

    • 'done'

    • 'none'

  • iteration ― 0 から始まる整数、またはNaN。バーンイン反復のある nlmefitsa では、出力関数は STATUS.iteration で負の値をもつ各反復の後で呼び出されます。

fval現在の対数尤度
Psi現在の変量効果の共分散行列
thetaPsi の現在の並列化
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 の [Stop] ボタンを押して出力関数を停止することができます。たとえば、以下のコードはダイアログを実装して計算をキャンセルします。

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 Toolbox 出力関数の例です。これにより、固定効果 (BETA) および変量効果 (diag(STATUS.Psi)) の分散をもつプロットを初期化および更新します。nlmefit では、プロットには対数尤度 (STATUS.fval) も含まれています。

nlmefitoutputfcnnlmefitsa の既定の設定の出力関数です。これを nlmefit で使用するには、options 構造体でその関数ハンドルを指定します。

opt = statset('OutputFcn', @nlmefitoutputfcn, …)
beta = nlmefit(…, 'Options', opt, …)

nlmefitsa がこの関数を使用しないようにするには、出力関数に空の値を指定します。

opt = statset('OutputFcn', [], …)
beta = nlmefitsa(…, 'Options', opt, …)

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

nlmefit と nlmefitsa を使用した混合効果のモデル

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

load indomethacin

indomethacin.mat のデータは、被験者 6 人の血流に含まれる薬物インドメタシンの濃度を 8 時間にわたって記録しています。

血流中のインドメタシンの散布図を被験者ごとにグループ化してプロットします。

gscatter(time,concentration,subject)
xlabel('Time (hours)')
ylabel('Concentration (mcg/ml)')
title('{\bf Indomethacin Elimination}')
hold on

Specifying Mixed-Effects Models」のページで、このタイプのデータに役立つモデルについて説明しています。

無名関数を使用してモデルを構成します。

model = @(phi,t)(phi(1)*exp(-exp(phi(2))*t) + ...
                 phi(3)*exp(-exp(phi(4))*t));

関数 nlinfit を使用して、被験者に固有の効果を無視し、モデルをデータのすべてに近似させます。

phi0 = [1 2 1 1];
[phi,res] = nlinfit(time,concentration,model,phi0);

二乗平均誤差を計算します。

numObs = length(time);
numParams = 4;
df = numObs-numParams;
mse = (res'*res)/df
mse =

    0.0304

データの散布図にこのモデルを重ね合わせます。

tplot = 0:0.01:8;
plot(tplot,model(phi,tplot),'k','LineWidth',2)
hold off

残差を被験者別に箱ひげ図で描画します。

colors = 'rygcbm';
h = boxplot(res,subject,'colors',colors,'symbol','o');
set(h(~isnan(h)),'LineWidth',2)
hold on
boxplot(res,subject,'colors','k','symbol','ko')
grid on
xlabel('Subject')
ylabel('Residual')
hold off

被験者別の残差の箱ひげ図は、ボックスがほとんど 0 の上か下にあることを示し、モデルが被験者に固有の効果を説明できなかったことを示します。

被験者固有の効果を考慮するには、被験者ごとにモデルを個別にデータに近似させます。

phi0 = [1 2 1 1];
PHI = zeros(4,6);
RES = zeros(11,6);
for I = 1:6
    tI = time(subject == I);
    cI = concentration(subject == I);
    [PHI(:,I),RES(:,I)] = nlinfit(tI,cI,model,phi0);
end
PHI
PHI =

    2.0293    2.8277    5.4683    2.1981    3.5661    3.0023
    0.5794    0.8013    1.7498    0.2423    1.0408    1.0882
    0.1915    0.4989    1.6757    0.2545    0.2915    0.9685
   -1.7878   -1.6354   -0.4122   -1.6026   -1.5069   -0.8731

二乗平均誤差を計算します。

numParams = 24;
df = numObs-numParams;
mse = (RES(:)'*RES(:))/df
mse =

    0.0057

データの散布図をプロットし、被験者別のモデルを重ね合わせます。

gscatter(time,concentration,subject)
xlabel('Time (hours)')
ylabel('Concentration (mcg/ml)')
title('{\bf Indomethacin Elimination}')
hold on
for I = 1:6
    plot(tplot,model(PHI(:,I),tplot),'Color',colors(I))
end
axis([0 8 0 3.5])
hold off

PHI は、6 人の被験者ごとに 4 つのモデル パラメーターの推定を与えます。推定はかなりばらつきますが、データの 24 パラメーター モデルの採用により、二乗平均誤差 0.0057 は、元の 4 パラメーター モデルの 0.0304 よりもかなり減少しています。

被験者別の残差の箱ひげ図を描画します。

h = boxplot(RES,'colors',colors,'symbol','o');
set(h(~isnan(h)),'LineWidth',2)
hold on
boxplot(RES,'colors','k','symbol','ko')
grid on
xlabel('Subject')
ylabel('Residual')
hold off

より大きなモデルによってほとんどの被験者固有の効果が説明できることが箱ひげ図で示されました。残差の広がり (箱ひげ図の垂直のスケール) は前の箱ひげ図よりずっと小さく、ボックスはほぼ 0 に集中するようになります。

24 パラメーター モデルは、研究対象の特定の被験者によって変動を説明できますが、この被験者がより大きい母集団の代表とは見なされません。被験者が描かれる標本分布は、おそらく標本自体より興味深いでしょう。混合効果のモデルの目的は、被験者に固有の変動を母集団の平均付近で変化する変量効果として、より大まかに説明することです。

関数 nlmefit を使用して、混合効果モデルをデータに近似します。nlmefit の代わりに nlmefitsa を使用することもできます。

次の無名関数 nlme_model は、異なるパラメーターを個々に使用できるようにすることで、nlinfit が使用する 4 パラメーターのモデルを nlmefit の呼び出し構文に適応させます。既定の設定では、nlmefit は変量効果をすべてのモデル パラメーターに割り当てます。また、既定の設定で、過剰なパラメーター化と関連収束問題を避けるために、nlmefit は対角線の共分散行列 (変量効果の間の無共分散) を仮定します。

nlme_model = @(PHI,t)(PHI(:,1).*exp(-exp(PHI(:,2)).*t) + ...
                      PHI(:,3).*exp(-exp(PHI(:,4)).*t));
phi0 = [1 2 1 1];
[phi,PSI,stats] = nlmefit(time,concentration,subject, ...
                          [],nlme_model,phi0)
phi =

    2.8277
    0.7729
    0.4606
   -1.3459


PSI =

    0.3264         0         0         0
         0    0.0250         0         0
         0         0    0.0124         0
         0         0         0    0.0000


stats = 

           dfe: 57
          logl: 54.5882
           mse: 0.0066
          rmse: 0.0787
    errorparam: 0.0815
           aic: -91.1765
           bic: -93.0506
          covb: [4x4 double]
        sebeta: [0.2558 0.1066 0.1092 0.2244]
          ires: [66x1 double]
          pres: [66x1 double]
         iwres: [66x1 double]
         pwres: [66x1 double]
         cwres: [66x1 double]

二乗平均誤差 0.0066 は、変量効果なしの 4 パラメーター モデルの 0.0304 よりかなり良く、変量効果なしの 24 パラメーター モデルの 0.0057 に匹敵します。

推定された共分散行列 PSI は 4 つ目の変量効果の分散が本質的に 0 であることを示し、モデルを簡略化するために削除できることを示唆します。これを行うために、名前と値のペア 'REParamsSelect' を使用して、nlmefit 内で変量効果でモデル化するパラメーターのインデックスを指定します。

[phi,PSI,stats] = nlmefit(time,concentration,subject, ...
                          [],nlme_model,phi0, ...
                          'REParamsSelect',[1 2 3])
phi =

    2.8277
    0.7728
    0.4605
   -1.3460


PSI =

    0.3270         0         0
         0    0.0250         0
         0         0    0.0124


stats = 

           dfe: 58
          logl: 54.5875
           mse: 0.0066
          rmse: 0.0780
    errorparam: 0.0815
           aic: -93.1750
           bic: -94.8410
          covb: [4x4 double]
        sebeta: [0.2560 0.1066 0.1092 0.2244]
          ires: [66x1 double]
          pres: [66x1 double]
         iwres: [66x1 double]
         pwres: [66x1 double]
         cwres: [66x1 double]

対数尤度 logl は、すべてのパラメーターが変量効果のものとほぼ同一です。赤池情報量基準 aic は -91.1765 から -93.1750 に減少します。そして、ベイズ情報量基準 bic は -93.0506 から -94.8410 に減少します。これらの方法は、4 つ目の変量効果を棄却する決定をサポートします。

完全な共分散行列で単純化したモデルの再近似により、変量効果の間で相関関係を識別することができます。これを行うには、CovPattern パラメーターを使用して、共分散行列内で非ゼロの要素のパターンを指定します。

[phi,PSI,stats] = nlmefit(time,concentration,subject, ...
                          [],nlme_model,phi0, ...
                          'REParamsSelect',[1 2 3], ...
                          'CovPattern',ones(3))
phi =

    2.8149
    0.8292
    0.5613
   -1.1409


PSI =

    0.4767    0.1151    0.0499
    0.1151    0.0321    0.0032
    0.0499    0.0032    0.0236


stats = 

           dfe: 55
          logl: 58.4726
           mse: 0.0061
          rmse: 0.0782
    errorparam: 0.0781
           aic: -94.9451
           bic: -97.2358
          covb: [4x4 double]
        sebeta: [0.3028 0.1103 0.1179 0.1662]
          ires: [66x1 double]
          pres: [66x1 double]
         iwres: [66x1 double]
         pwres: [66x1 double]
         cwres: [66x1 double]

推定された共分散行列 PSI は、最初の 2 つのパラメーターにおける変量効果が比較的強い相関関係をもっていること、そして、両方とも最後の無作為な効果で比較的弱い相関関係をもっていることを示します。corrcov を使用して、PSI を相関関係行列に変換すると、この共分散行列内の構造はより明白になります。

RHO = corrcov(PSI)
clf;
imagesc(RHO)
set(gca,'XTick',[1 2 3],'YTick',[1 2 3])
title('{\bf Random Effect Correlation}')
h = colorbar;
set(get(h,'YLabel'),'String','Correlation');
RHO =

    1.0000    0.9315    0.4706
    0.9315    1.0000    0.1175
    0.4706    0.1175    1.0000

共分散パターンの仕様をブロック対角に変更することにより、この構造をモデルに組み入れます。

P = [1 1 0;1 1 0;0 0 1] % Covariance pattern
[phi,PSI,stats,b] = nlmefit(time,concentration,subject, ...
                            [],nlme_model,phi0, ...
                            'REParamsSelect',[1 2 3], ...
                            'CovPattern',P)
P =

     1     1     0
     1     1     0
     0     0     1


phi =

    2.7830
    0.8981
    0.6581
   -1.0000


PSI =

    0.5180    0.1069         0
    0.1069    0.0221         0
         0         0    0.0454


stats = 

           dfe: 57
          logl: 58.0804
           mse: 0.0061
          rmse: 0.0768
    errorparam: 0.0782
           aic: -98.1608
           bic: -100.0350
          covb: [4x4 double]
        sebeta: [0.3171 0.1073 0.1384 0.1453]
          ires: [66x1 double]
          pres: [66x1 double]
         iwres: [66x1 double]
         pwres: [66x1 double]
         cwres: [66x1 double]


b =

   -0.8507   -0.1563    1.0427   -0.7559    0.5652    0.1550
   -0.1756   -0.0323    0.2152   -0.1560    0.1167    0.0320
   -0.2756    0.0519    0.2620    0.1064   -0.2835    0.1389

ブロック対角共分散構造は、対数尤度に大きな影響を与えずに、aic を -94.9462 から -98.1608 に減らし、bic を -97.2368 から -100.0350 に減らします。これらの方法は、最終モデル内で使用される共分散構造をサポートします。出力 b は、6 人の被験者ごとに 3 つの変量効果の予測を与えます。これらは phi における固定効果の推定と結合され、混合効果モデルを生成します。

6 人の被験者別に混合効果モデルをプロットします。比較するために、変量効果のないモデルも示されます。

PHI = repmat(phi,1,6) + ...                 % Fixed effects
      [b(1,:);b(2,:);b(3,:);zeros(1,6)];    % Random effects
RES = zeros(11,6); % Residuals
colors = 'rygcbm';
for I = 1:6
    fitted_model = @(t)(PHI(1,I)*exp(-exp(PHI(2,I))*t) + ...
                        PHI(3,I)*exp(-exp(PHI(4,I))*t));
    tI = time(subject == I);
    cI = concentration(subject == I);
    RES(:,I) = cI - fitted_model(tI);

    subplot(2,3,I)
    scatter(tI,cI,20,colors(I),'filled')
    hold on
    plot(tplot,fitted_model(tplot),'Color',colors(I))
    plot(tplot,model(phi,tplot),'k')
    axis([0 8 0 3.5])
    xlabel('Time (hours)')
    ylabel('Concentration (mcg/ml)')
    legend(num2str(I),'Subject','Fixed')
end

データ内の明らかな外れ値 (前の箱ひげ図を参照) が無視されると、残差の正規確率プロットは、誤差のモデル仮定と妥当な一致を示します。

clf; normplot(RES(:))

モデルの検証のための残差の検査

nlmefitnlmefitsa のどちらによっても返される stats 構造体を検査することによって、モデルの特性を判断できます。stats 構造体には、条件付きの重み付き残差 (cwres フィールド) のフィールドと、個別の重み付け残差 (iwres フィールド) のフィールドがあります。モデルでは、残差が正規分布されていると仮定されているため、残差を検証することによって、この仮定がどの程度有効であるかを確認できます。

この例では、正規分布を使用して合成されたデータを生成します。近似統計量がどのようになるかを次に示します。

  • データを生成したのと同じタイプのモデルに対するテストでは良好

  • 無効なデータ モデルに対するテストでは不良

  1. 100 個の個体で 2 次元モデルを初期化します。

    nGroups = 100; % 100 Individuals
    nlmefun = @(PHI,t)(PHI(:,1)*5 + PHI(:,2)^2.*t); % Regression fcn
    REParamSelect = [1  2]; % Both Parameters have random effect
    errorParam = .03; 
    beta0 = [ 1.5  5]; % Parameter means
    psi = [ 0.35  0; ...  % Covariance Matrix
           0   0.51 ];
    time =[0.25;0.5;0.75;1;1.25;2;3;4;5;6];
    nParameters = 2;
    rng(0,'twister') % for reproducibility
  2. 比例誤差モデルの近似に使用するデータを生成します。

    b_i = mvnrnd(zeros(1, numel(REParamSelect)), psi, nGroups);
    individualParameters = zeros(nGroups,nParameters);
    individualParameters(:, REParamSelect) = ...
         bsxfun(@plus,beta0(REParamSelect), b_i);
    
    groups = repmat(1:nGroups,numel(time),1);
    groups = vertcat(groups(:));
    
    y = zeros(numel(time)*nGroups,1);
    x = zeros(numel(time)*nGroups,1);
    for i = 1:nGroups
        idx = groups == i;
        f = nlmefun(individualParameters(i,:), time);
        % Make a proportional error model for y:
        y(idx) = f + errorParam*f.*randn(numel(f),1);
        x(idx) = time;
    end
    
    P = [ 1 0 ; 0 1 ];
  3. モデル ジェネレーターと同じ回帰関数と誤差モデルを使用してデータを近似します。

    [~,~,stats] = nlmefit(x,y,groups, ...
        [],nlmefun,[1 1],'REParamsSelect',REParamSelect,...
        'ErrorModel','Proportional','CovPattern',P);
  4. 次の関数定義をコピーしてプロット ルーチンを作成し、ファイル plotResiduals.m を MATLAB パス上に作成します。

    function plotResiduals(stats)
    pwres = stats.pwres;
    iwres = stats.iwres;
    cwres = stats.cwres;
    figure
    subplot(2,3,1);
    normplot(pwres); title('PWRES')
    subplot(2,3,4);
    createhistplot(pwres);
    
    subplot(2,3,2);
    normplot(cwres); title('CWRES')
    subplot(2,3,5);
    createhistplot(cwres);
    
    subplot(2,3,3);
    normplot(iwres); title('IWRES')
    subplot(2,3,6);
    createhistplot(iwres); title('IWRES')
    
    function createhistplot(pwres)
    [x, n] = hist(pwres);
    d = n(2)- n(1);
    x = x/sum(x*d);
    bar(n,x);
    ylim([0 max(x)*1.05]);
    hold on;
    x2 = -4:0.1:4;
    f2 = normpdf(x2,0,1);
    plot(x2,f2,'r');
    end
    
    end
  5. 関数 plotResiduals を使用して残差をプロットします。

    plotResiduals(stats);

    上段の確率プロットは直線に見えますが、これは残差が正規分布をしていることを意味します。下段のヒストグラム プロットは重ねて表示された正規密度プロットと一致しています。したがって、誤差モデルとデータが一致していると結論づけることができます。

  6. 比較のために、データ作成時と同じ比例モデルではなく、定誤差モデルを使用してデータを近似します。

    [~,~,stats] = nlmefit(x,y,groups, ...
        [],nlmefun,[0 0],'REParamsSelect',REParamSelect,...
        'ErrorModel','Constant','CovPattern',P);
    plotResiduals(stats);

    上段の確率プロットは直線ではなく、残差が正規分布をしていないことを示しています。下段のヒストグラム プロットは重ねて表示された正規密度プロットとほとんど一致しています。

  7. さらに別の比較をするために、データ作成時とは異なる構造モデルにデータを近似します。

    nlmefun2 = @(PHI,t)(PHI(:,1)*5 + PHI(:,2).*t.^4);
    [~,~,stats] = nlmefit(x,y,groups, ...
        [],nlmefun2,[0 0],'REParamsSelect',REParamSelect,...
        'ErrorModel','constant', 'CovPattern',P);
    plotResiduals(stats);

    上図の確率プロットは直線を示していないだけでなく、ヒストグラム プロットも、重ねて表示された正規密度と比較すると大きく歪んでいます。これらの残差は正規分布しておらず、モデルとは一致していません。

この情報は役に立ちましたか?