compare

クラス: GeneralizedLinearMixedModel

一般化線形混合効果モデルの比較

構文

results = compare(glme,altglme)

results = compare(glme,altglme,CheckNesting=checknesting)

説明

results = compare(glme,altglme) は、一般化線形混合効果モデル glme と altglme を比較する尤度比検定の結果を返します。有効な尤度比検定を実行するためには、どちらのモデルも近似で同じ応答ベクトルを使用しなければならず、glme を altglme の入れ子にしなければなりません。常に、小さい方のモデルを最初に入力し、次に大きい方のモデルを入力します。

compare は次の帰無仮説と対立仮説を検定します。

• H₀:観測した応答ベクトルは glme によって生成された。

• H₁:観測した応答ベクトルは altglme モデルによって生成された。

results = compare(glme,altglme,CheckNesting=checknesting) は、1 番目の入力モデル glme が 2 番目の入力モデル altglme に入れ子になっているかどうかを調べます。

入力引数

glme — 一般化線形混合効果モデル

GeneralizedLinearMixedModel オブジェクト

一般化線形混合効果モデル。GeneralizedLinearMixedModel オブジェクトとして指定します。このオブジェクトのプロパティとメソッドについては、GeneralizedLinearMixedModel を参照してください。

fitglme を使用して、一般化線形混合効果モデルを標本データに当てはめて GeneralizedLinearMixedModel オブジェクトを作成することができます。正規分布以外の応答分布をもつ 2 つのモデルに有効な尤度比検定を実行するには、両方のモデルを "ApproximateLaplace" または "Laplace" 近似メソッドを使用して当てはめなくてはなりません。"MPL" または "REMPL" を使用して当てはめる、正規分布以外の応答分布となるモデルは、尤度比検定を使用して比較できません。

altglme — 代替の一般化線形混合効果モデル

GeneralizedLinearMixedModel オブジェクト

代替一般化線形混合効果モデル。GeneralizedLinearMixedModel オブジェクトとして指定します。altglme は glme と同じ応答ベクトルに当てはめますが、モデル仕様は異なっています。altglme の一部のモデル パラメーターを 0 などの固定値に設定して、glme を altglme から取得できるように glme を altglme において入れ子にしなければなりません。

fitglme を使用して、一般化線形混合効果モデルを標本データに当てはめて GeneralizedLinearMixedModel オブジェクトを作成することができます。正規分布以外の応答分布をもつ 2 つのモデルに有効な尤度比検定を実行するには、両方のモデルを "ApproximateLaplace" または "Laplace" 近似メソッドを使用して当てはめなくてはなりません。"MPL" または "REMPL" を使用して当てはめる、正規分布以外の応答分布となるモデルは、尤度比検定を使用して比較できません。

checknesting — 2 つのモデルの間の入れ子を調べるためのインジケーター

true (既定値) | false

2 つのモデルの間の入れ子を調べるためのインジケーター。true または false のいずれかとして指定します。checknesting が true の場合、compare は小さい方のモデル glme が大きい方のモデル altglme に入れ子になっているかどうかを確認します。入れ子要件が満たされない場合、compare はエラーを返します。checknesting が false の場合、compare はこのチェックを実行しません。

例: CheckNesting=true

出力引数

results — 尤度比検定結果

table

尤度比検定結果。2 行の table として返されます。1 行目は glme に対応し、2 行目は altglme に対応しています。results の列は以下を含みます。

列名	説明
Model	モデルの名前
DF	自由度
AIC	モデルの赤池情報量基準
BIC	モデルのベイズ情報量基準
LogLik	モデルの最大化された対数尤度
LRStat	altglme と glme を比較するための尤度比検定統計
deltaDF	altglme の DF から glme の DF を引いた値
pValue	尤度比検定の p 値

例

混合効果モデルの比較

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

load mfr

このシミュレーションされたデータは、世界中で 50 の工場を操業している製造企業から取得しており、各工場が完成品の生産のためにバッチ処理を実行しています。同社は各バッチの欠陥数を減少させるために新たな製造プロセスを開発しました。新しいプロセスの効果をテストするため、同社は実験に参加させる 20 工場を無作為に選びました。10 工場では新プロセスを実施しますが、残りの 10 工場では旧プロセスの実行を続けます。各 20 工場で、同社は 5 つのバッチ (合計 100 バッチ) を実行し以下のデータを記録しました。

• 新しいプロセスがバッチに使用されたかどうかを示すフラグ (newprocess)

• 各バッチの処理時間。時間単位 (time)

• バッチの温度。摂氏 (temp)

• バッチに使用する化学薬品の供給業者を示すカテゴリカル変数 (supplier)

• バッチ内の欠陥数 (defects)

またデータに含まれる time_dev と temp_dev は、摂氏 20 度で 3 時間の標準プロセスから得られる時間と温度の絶対偏差をそれぞれ表します。

固定効果予測子として newprocess、time_dev、temp_dev および supplier を使用して固定効果専用モデルを当てはめます。応答分布にポアソンを、リンク関数に対数を、近似メソッドにラプラスを指定します。ダミー変数エンコードを 'effects' として指定すると、ダミー変数の係数の合計が 0 になります。

FEglme = fitglme(mfr,'defects ~ 1 + newprocess + time_dev + temp_dev + supplier','Distribution','Poisson','Link','log','FitMethod','Laplace','DummyVarCoding','effects');

同じ固定効果予測子、応答分布、リンク関数および近似メソッドを使用する 2 番目のモデルによる当てはめを実行します。今回は factory 別にグループ化された変量効果の切片を含めます。これは工場特有の変動に起因して品質に差がある可能性を考慮するためです。

欠陥数はポアソン分布を使用してモデル化できます

$${\text{defects}}_{ij}\sim \text{Poisson}({\mu }_{ij})$$

これは一般化線形混合効果モデルに対応します

$$\log ({\mu }_{ij})={\beta }_0+{\beta }_1{\text{newprocess}}_{ij}+{\beta }_2{\text{time}\text{\_}\text{dev}}_{ij}+{\beta }_3{\text{temp}\text{\_}\text{dev}}_{ij}+{\beta }_4{\text{supplier}\text{\_}\text{C}}_{ij}+{\beta }_5{\text{supplier}\text{\_}\text{B}}_{ij}+b_i,$$

ここで

• $${\text{defects}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ で実行されたバッチで観測された欠陥数です。

• $${\mu }_{ij}$$ は、バッチ $$j$$ ($$j=1,2,...,5$$) 処理中の工場 $$i$$ ($$i=1,2,...,20$$) に対応する欠陥の平均数です。

• $${\text{newprocess}}_{ij}$$、$${\text{time}\text{\_}\text{dev}}_{ij}$$ および $${\text{temp}\text{\_}\text{dev}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ に対応する各変数の測定値です。たとえば $${\text{newprocess}}_{ij}$$ は、工場 $$i$$ で実行されたバッチ $$j$$ 処理中に新プロセスが使用されたかどうかを示します。

• $${\text{supplier}\text{\_}\text{C}}_{ij}$$ および $${\text{supplier}\text{\_}\text{B}}_{ij}$$ はエフェクト (ゼロサム) コーディングを使用するダミー変数であり、バッチ $$j$$ 処理中に工場 $$i$$ で実行されたバッチに対して、それぞれ会社 C または B が加工化学薬品を供給したかどうかを示します。

• $$b_i\sim N(0,{\sigma }_b^2)$$ は、工場特有の品質変動に相当する、各工場 $$i$$ の変量効果の切片です。

glme = fitglme(mfr,'defects ~ 1 + newprocess + time_dev + temp_dev + supplier + (1|factory)','Distribution','Poisson','Link','log','FitMethod','Laplace','DummyVarCoding','effects');

理論的尤度比検定法を使用して 2 つのモデルを比較します。'CheckNesting' を true として指定すると、compare は入れ子要件が満たされていない場合に警告を返します。

results = compare(FEglme,glme,'CheckNesting',true)

results = 
    Theoretical Likelihood Ratio Test

    Model     DF    AIC       BIC       LogLik     LRStat    deltaDF    pValue    
    FEglme    6     431.02    446.65    -209.51                                   
    glme      7     416.35    434.58    -201.17    16.672    1          4.4435e-05

compare はエラーを返さなかったため、入れ子要件は満たされています。小さい $$p$$ 値は、観測された応答ベクトルがモデル FEglme により生成されたという帰無仮説が compare により棄却され、代わりに代替モデル glme が採択されることを示します。glme に対する AIC と BIC の値の方が小さいことも、glme の方が応答により当てはまるモデルを提供するという結論を裏付けています。

詳細

尤度比検定

"尤度比検定" は 2 つの入れ子モデルの仕様を比較するため、パラメーターによって制限されずに拡張されたモデルで制限の有意性を評価します。帰無仮説 H₀ の場合、尤度比検定統計は、自由度が deltaDF の近似カイ二乗参照分布になります。

2 つのモデルを比較するとき、compare は、観測される尤度比検定統計とこのカイ二乗参照分布を比較して、尤度比検定の p 値を計算します。p 値が小さい場合は H₁ を優先して H₀ は棄却され、代替モデル altglme が採択されます。一方、p 値が大きい場合は、H₀ が棄却できないことを示し、モデル altglme を採択するのに十分な証拠がないことを反映しています。

尤度比検定を使用して取得される p 値は、変量効果の項の有無について検定する場合は保守的になり、固定効果の項の有無について検定する場合は保守的でなくなる可能性もあります。代わりに、fixedEffects または coefTest メソッドを固定効果の検定に使用してください。

GLME モデルに有効な尤度比検定を実行するには、両方のモデルはラプラスまたは近似ラプラス近似メソッドを使用して当てはめなければなりません。疑似最尤 (MPL) または制限付き疑似最尤 (REMPL) メソッドを使用して当てはめたモデルは尤度比検定を使用して比較することはできません。MPL を使用して当てはめたモデルを比較する場合、尤度比検定に使用されるのは、疑似尤度の最後の反復からの疑似データが、最大化された対数尤度です。MPL を使用して当てはめた非正規分布のモデルを比較する場合、compare から発生する警告により、尤度比検定で使用されるのは、疑似尤度の最後の反復からの疑似データが、最大化された対数尤度であることを示します。尤度比検定で真の最大化された対数尤度を使用するには、モデルを比較する前に、近似ラプラスまたはラプラスを使用して glme と altglme の両方を当てはめなければなりません。

入れ子要件

有効な尤度比検定を実行するには、glme が altglme に入れ子になっていなければなりません。CheckNesting=true 引数は以下の要件を確認し、満たさないものが 1 つでもあるとエラーを返します。

• 両方のモデル (glme と altglme) を "ApproximateLaplace" または "Laplace" 近似メソッドを使用して当てはめなければなりません。"MPL" または "REMPL" を使用して当てはめた GLME モデルは、尤度比検定を使用して比較することはできません。

• どちらのモデルも同じ応答ベクトル、応答分布およびリンク関数を使用して当てはめなければなりません。

• 小さい方のモデル (glme) は大きい方のモデル (altglme) の内部に入れ子にしなければなりません。altglme の一部のモデル パラメーターを 0 などの固定値に設定すると、glme を altglme から取得できるようになります。

• 大きい方のモデル (altglme) の最大化された対数尤度は、小さい方のモデル (glme) の最大化された対数尤度以上でなければなりません。

• glme の当てはめに使用された重みベクトルと altglme は同一でなければなりません。

• 大きい方のモデル (altglme) の変量効果の計画行列には、小さい方のモデル (glme) の変量効果の計画行列が含まれなければなりません。

• 大きい方のモデル (altglme) の固定効果の計画行列には、小さい方のモデル (glme) の固定効果の計画行列が含まれなければなりません。

赤池およびベイズ情報量基準

"赤池情報量基準" (AIC) は AIC = –2logL_M + 2(param) です。

logL_M はモデルの当てはめに使用する方法に依存します。

• 'Laplace' または 'ApproximateLaplace' を使用する場合、logL_M は最大化された対数尤度です。

• 'MPL' を使用する場合、logL_M は疑似尤度の最後の反復による疑似データに対する最大化された対数尤度です。

• 'REMPL' を使用する場合、logL_M は疑似尤度の最後の反復による疑似データに対する最大化された制限付き対数尤度です。

param はモデル内で推定されたパラメーターの総数です。ほとんどの GLME モデルで、param は nc + p + 1 に等しくなります。ここで、nc は残差分散を除く変量効果共分散のパラメーターの総数、p は固定効果係数の個数です。ただし、二項分布またはポアソン分布で分散パラメーターが 1.0 に固定されている場合、param は (nc + p) に等しくなります。

"ベイズ情報量基準" (BIC) は BIC = –2*logL_M + ln(n_eff)(param) です。

logL_M はモデルの当てはめに使用する方法に依存します。

• 'Laplace' または 'ApproximateLaplace' を使用する場合、logL_M は最大化された対数尤度です。

• 'MPL' を使用する場合、logL_M は疑似尤度の最後の反復による疑似データに対する最大化された対数尤度です。

• 'REMPL' を使用する場合、logL_M は疑似尤度の最後の反復による疑似データに対する最大化された制限付き対数尤度です。

n_eff は有効な観測値の個数です。

• 'MPL'、'Laplace' または 'ApproximateLaplace' を使用する場合、n_eff = n になります。ここで、n は観測値の個数です。

• 'REMPL' を使用する場合、n_eff = n – p になります。

param はモデル内で推定されたパラメーターの総数です。ほとんどの GLME モデルで、param は nc + p + 1 に等しくなります。ここで、nc は残差分散を除く変量効果共分散のパラメーターの総数、p は固定効果係数の個数です。ただし、二項分布またはポアソン分布で分散パラメーターが 1.0 に固定されている場合、param は (nc + p) に等しくなります。

逸脱度の値が小さいほど、近似が優れていることを意味します。逸脱度の値が小さくなると、AIC および BIC も小さくなる傾向があります。AIC と BIC のどちらにも、推定されるパラメーターの数 p に基づくペナルティ項が含まれます。したがって、パラメーターの数が増えると、AIC および BIC の値も大きくなる傾向があります。異なるモデルを比較するときは、AIC または BIC の値が最も小さいモデルが最良近似のモデルと考えられます。

'MPL' および 'REMPL' を使用して当てはめたモデルの場合、AIC および BIC は疑似尤度の最後の反復から派生した疑似データの対数尤度 (または制限付き対数尤度) に基づきます。したがって、'MPL' と 'REMPL' を使用して当てはめたモデル間で AIC と BIC の値を直接比較するのは適切ではありません。

参考

GeneralizedLinearMixedModel | covarianceParameters | fixedEffects | randomEffects