Main Content

compare

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

説明

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

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

H0:観測した応答ベクトルは lme によって生成された。

H1:観測した応答ベクトルは altlme モデルによって生成された。

モデルの比較に先立ち、ML (最尤) 法を使用して lmealtlme を当てはめることをお勧めします。REML (制限付き最尤) 法を使用する場合は、両方のモデルに含まれる固定効果の計画行列を同じにしなければなりません。

固定効果を検定するには、lmealtlme を ML で当てはめてるか、シミュレーションされた尤度比検定compare を使用するか、fixedEffectsanova または coefTest メソッドを使用します。

results = compare(___,Name,Value) も、線形混合効果モデル lmealtlme を比較する尤度比検定の結果を返しますが、1 つ以上の Name,Value ペア引数で指定する追加オプションがあります。

たとえば、最初の入力モデルが 2 番目の入力モデルに入れ子になっているかどうかを調べることができます。

[results,siminfo] = compare(lme,altlme,'NSim',nsim) は、線形混合効果モデル lmealtlme を比較するシミュレーションされた尤度比検定の結果を返します。

lmealtlme は ML または REML を使用して当てはめることができます。また、lmealtlme の入れ子にする必要はありません。REML (制限付き最尤) 法を使用してモデルを当てはめる場合は、両方のモデルが同じ固定効果計画行列をもっていなければなりません。

[results,siminfo] = compare(___,Name,Value) も、線形混合効果モデル lmealtlme を比較するシミュレーションされた尤度比検定の結果を返しますが、1 つ以上の Name,Value ペア引数で指定する追加オプションがあります。

たとえば、シミュレーションされた尤度比検定の実行に関するオプションを変更したり、p 値の信頼区間の信頼水準を変更したりできます。

すべて折りたたむ

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

load flu

データセット配列 flu には、変数 Date と、インフルエンザ推定罹患率 (Google® 検索から推定される 9 地域の値と CDC による全国の推定値) が格納されている 10 個の変数が含まれています。

線形混合効果モデルを当てはめるには、データが適切な形式のデータセット配列になっていなければなりません。インフルエンザ罹患率を応答として、地域を予測子変数として線形混合効果モデルを当てはめるため、地域に対応する 9 個の列を 1 つの配列にまとめます。新しいデータセット配列 flu2 には、応答変数 FluRate、各推定の元になっている地域を示すノミナル変数 Region およびグループ化変数 Date が含まれなければなりません。

flu2 = stack(flu,2:10,'NewDataVarName','FluRate',...
    'IndVarName','Region');
flu2.Date = nominal(flu2.Date);

地域ごとに切片と勾配を変化させ、Date でグループ化して、線形混合効果モデルを当てはめます。

altlme = fitlme(flu2,'FluRate ~ 1 + Region + (1 + Region|Date)');

地域に対する固定効果と、Date で変化するランダム切片で、線形混合効果モデルを当てはめます。

lme = fitlme(flu2,'FluRate ~ 1 + Region + (1|Date)');

2 つのモデルを比較します。また、lme2lme の入れ子になっているかどうかを調べます。

compare(lme,altlme,'CheckNesting',true)
ans = 
    THEORETICAL LIKELIHOOD RATIO TEST

    Model     DF    AIC        BIC        LogLik     LRStat    deltaDF    pValue
    lme       11     318.71     364.35    -148.36                               
    altlme    55    -305.51    -77.346     207.76    712.22    44         0     

0 という小さい p 値は、モデル altlme が単純なモデル lme より有意に優れていることを示します。

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

load('fertilizer.mat');

このデータセット配列には土壌の種類に基づいて土壌が 3 つのブロックに分けられている分割プロット試験のデータが含まれています。土壌の種類は砂質、シルトおよび粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。その後、プロット内のトマトの苗木はサブプロットに分割され、それぞれのサブプロットが 4 つの肥料の中の 1 つにより処置されます。このデータは、シミュレーションされたものです。

実用目的でこのデータを ds という名前のデータセット配列に保存し、TomatoSoil および Fertilizer をカテゴリカル変数として定義します。

ds = fertilizer;
ds.Tomato = nominal(ds.Tomato);
ds.Soil = nominal(ds.Soil);
ds.Fertilizer = nominal(ds.Fertilizer);

線形混合効果モデルを当てはめます。Fertilizer および Tomato は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。

lmeBig = fitlme(ds,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)');

交互作用項 Tomato:Fertilizer および変量効果の項 (1 | Soil) を削除した後、モデルを再度当てはめます。

lmeSmall = fitlme(ds,'Yield ~ Fertilizer + Tomato + (1|Soil:Tomato)');

シミュレーションされた尤度比検定を 1000 反復で使用して 2 つのモデルを比較します。この検定を使用して、固定効果の項と変量効果の項の両方を検定しなければなりません。既定の近似法 ML を使用して両方のモデルを当てはめることに注意してください。このため、固定効果計画行列に対して制限はありません。REML (制限付き最尤) 法を使用する場合は、両方のモデルが同じ固定効果計画行列をもっていなければなりません。

[table,siminfo] = compare(lmeSmall,lmeBig,'nsim',1000)
table = 
    SIMULATED LIKELIHOOD RATIO TEST: NSIM = 1000, ALPHA = 0.05

    Model       DF    AIC       BIC       LogLik     LRStat    pValue     Lower      Upper  
    lmeSmall    10    511.06       532    -245.53                                           
    lmeBig      23    522.57    570.74    -238.29    14.491    0.57343    0.54211    0.60431

siminfo = struct with fields:
           nsim: 1000
          alpha: 0.0500
      pvalueSim: 0.5734
    pvalueSimCI: [0.5421 0.6043]
        deltaDF: 13
            TH0: [1000x1 double]

大きい p 値は、大きいモデル lme が小さいモデル lme2 より有意には優れていないことを示します。lme2 の方が赤池およびベイズ情報量基準の値が小さいこともこれを裏付けています。

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

load carbig

ガロンあたりの走行マイル数 (MPG) の線形混合効果モデルを当てはめます。この近似では、加速度、馬力、気筒数に対する固定効果と、モデル年度によってグループ化される切片と加速度に対する、相関された可能性がある変量効果を使用します。

最初に、計画行列を準備します。

X = [ones(406,1) Acceleration Horsepower];
Z = [ones(406,1) Acceleration];
Model_Year = nominal(Model_Year);
G = Model_Year;

次に、定義した計画行列とグループ化変数で fitlmematrix を使用してモデルを当てはめます。

lme = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Horsepower'},'RandomEffectPredictors',...
{{'Intercept','Acceleration'}},'RandomEffectGroups',{'Model_Year'});

切片と加速度について無相関の変量効果でモデルを再度当てはめます。最初に、変量効果計画と変量効果グループ化変数を準備します。

Z = {ones(406,1),Acceleration};
G = {Model_Year,Model_Year};

lme2 = fitlmematrix(X,MPG,Z,G,'FixedEffectPredictors',....
{'Intercept','Acceleration','Horsepower'},'RandomEffectPredictors',...
{{'Intercept'},{'Acceleration'}},'RandomEffectGroups',...
{'Model_Year','Model_Year'});

シミュレーションされた尤度比検定を使用して lmelme2 を比較します。

compare(lme2,lme,'CheckNesting',true,'NSim',1000)
ans = 


    SIMULATED LIKELIHOOD RATIO TEST: NSIM = 1000, ALPHA = 0.05

    Model    DF    AIC       BIC       LogLik     LRStat    pValue      Lower   
    lme2     6     2194.5    2218.3    -1091.3                                  
    lme      7     2193.5    2221.3    -1089.7    3.0323    0.095904    0.078373


    Upper  
           
    0.11585

大きい $p$ 値は、lme2lme より有意に優れた当てはめではないことを示します。

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

load('fertilizer.mat')

このデータセット配列には土壌の種類に基づいて土壌が 3 つのブロックに分けられている分割プロット試験のデータが含まれています。土壌の種類は砂質、シルトおよび粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。その後、プロット内のトマトの苗木はサブプロットに分割され、それぞれのサブプロットが 4 つの肥料の中の 1 つにより処置されます。このデータは、シミュレーションされたものです。

このデータを tbl という名前のテーブルに保存し、TomatoSoil および Fertilizer をカテゴリカル変数として定義します。

tbl = dataset2table(fertilizer);
tbl.Tomato = categorical(tbl.Tomato);
tbl.Soil = categorical(tbl.Soil);
tbl.Fertilizer = categorical(tbl.Fertilizer);

線形混合効果モデルを当てはめます。Fertilizer および Tomato は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。

lme = fitlme(tbl,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)');

交互作用項 Tomato:Fertilizer および変量効果の項 (1|Soil) を削除した後、モデルを再度当てはめます。

lme2 = fitlme(tbl,'Yield ~ Fertilizer + Tomato + (1|Soil:Tomato)');

LinearMixedModel のオプション構造体を作成します。

opt = statset('LinearMixedModel')
opt = struct with fields:
          Display: 'off'
      MaxFunEvals: []
          MaxIter: 10000
           TolBnd: []
           TolFun: 1.0000e-06
       TolTypeFun: []
             TolX: 1.0000e-12
         TolTypeX: []
          GradObj: []
         Jacobian: []
        DerivStep: []
      FunValCheck: []
           Robust: []
     RobustWgtFun: []
           WgtFun: []
             Tune: []
      UseParallel: []
    UseSubstreams: []
          Streams: {}
        OutputFcn: []

並列検定用に options を変更します。

opt.UseParallel = true;

並列環境を開始します。

mypool = parpool();
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

反復数 1000 の並列計算でシミュレーションされた尤度比検定を使用し、lme2lme を比較します。

compare(lme2,lme,'nsim',1000,'Options',opt)
ans = 
    Simulated Likelihood Ratio Test: Nsim = 1000, Alpha = 0.05

    Model    DF    AIC       BIC       LogLik     LRStat    pValue     Lower    Upper  
    lme2     10    511.06       532    -245.53                                         
    lme      23    522.57    570.74    -238.29    14.491    0.53447    0.503    0.56573

大きい p 値は、大きいモデル lme が小さいモデル lme2 より有意には優れていないことを示します。lme2 に対する AICBIC の値が小さいこともこれを裏付けています。

入力引数

すべて折りたたむ

線形混合効果モデル。fitlme または fitlmematrix を使用して構築した LinearMixedModel オブジェクトとして指定します。

同じ応答ベクトルに異なるモデル仕様によって当てはめる代替線形混合効果モデル。LinearMixedModel オブジェクトとして指定されます。lmealtlme で入れ子にされている必要があります。つまり、一部のパラメーターを 0 などの固定値に設定することにより lmealtlme から取得する必要があります。線形混合効果オブジェクトは、fitlme または fitlmematrix を使用して作成できます。

シミュレーションされた尤度比検定でのシミュレーションの反復数。正の整数値として指定します。シミュレーションされた尤度比検定を実行するには nsim を指定しなければなりません。

例: 'NSim',1000

データ型: double | single

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後ろにする必要がありますが、ペアの順序は関係ありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: results = compare(lme,altlme,'CheckNesting',true)

有意水準。'Alpha' と 0 ~ 1 の範囲にあるスカラー値から構成されるコンマ区切りのペアとして指定します。値が α の場合、信頼水準は 100 × (1 – α)% です。

たとえば、99% の信頼区間の場合は、次のように信頼水準を指定できます。

例: 'Alpha',0.01

データ型: single | double

シミュレーションされた尤度比検定を並列実行するためのオプション。'Options' と、statset('LinearMixedModel') によって作成される構造体から構成されるコンマ区切りのペアとして指定します。

これらのオプションには Parallel Computing Toolbox™ が必要です。

compare は次のフィールドを使用します。

'UseParallel'
  • 逐次計算の場合は False。既定の設定。

  • 並列計算の場合は True

並列計算には Parallel Computing Toolbox が必要です。

'UseSubstreams'
  • 各反復に対して乱数発生器の個別のサブストリームを使用しない場合は False。既定の設定。

  • 各反復に対して乱数発生器の個別のサブストリームを使用する場合は True。このオプションは、サブストリームをサポートする乱数ストリーム タイプでのみ使用できます。

'Streams'
  • 'UseSubstreams'True の場合、'Streams' は単一の乱数ストリームであるか、単一のストリームを含むスカラー cell 配列でなければなりません。

  • 'UseSubstreams'False の場合は、次のようになります。

    • 'UseParallel'False の場合は、'Streams' は単一の乱数ストリームであるか、単一のストリームを含むスカラー cell 配列でなければなりません。

    • 'UseParallel'True の場合、'Streams' は使用されるプロセッサ数と等しくなければなりません。並列プールが開いている場合、'Streams' は並列プールのサイズと同じ長さです。'UseParallel'True の場合、並列プールは自動的に開いている可能性があります。ただし、'Streams' を使用されるプロセッサの数と等しくしなければならないので、parpool コマンドを使用してプールを明示的に開き、次に 'UseParallel','True' オプションを指定して compare を呼び出すのが最善の方法です。

コマンド ラインでの並列統計計算の詳細を表示するには、次のように入力します。

help parallelstats

データ型: struct

2 つのモデルの間の入れ子を確認するためのインジケーター。'CheckNesting' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

false既定の設定。確認しません。
truecompare は、小さい方のモデル lme が大きい方のモデル altlme に入れ子になっているかどうかを確認します。

理論的尤度比検定を有効にするためには、lme を代替モデル altlme で入れ子にしなければなりません。入れ子要件が満たされていない場合、compare はエラー メッセージを返します。

どちらの検定に対しても有効ですが、シミュレーションされた尤度比検定の方が入れ子の要件は弱くなります。

例: 'CheckNesting',true

データ型: single | double

出力引数

すべて折りたたむ

尤度比検定またはシミュレーションされた尤度比検定の結果。2 行のデータセット配列として返します。1 行目は lme に対応し、2 行目は altlme に対応しています。results の列は、検定が尤度比検定かシミュレーションされた尤度比検定かに依存します。

  • 尤度比検定を使用している場合は、results に以下の列が含まれます。

    Modelモデルの名前
    DF自由度、つまりモデル内の自由パラメーターの数
    AICモデルの赤池情報量基準
    BICモデルのベイズ情報量基準
    LogLikモデルの最大化された対数尤度
    LRStataltlmelme を比較するための尤度比検定統計
    deltaDFaltlmeDF から lmeDF を引いた値
    pValue尤度比検定の p 値
  • シミュレーションされた尤度比検定を使用している場合は、results に以下の列が含まれます。

    Modelモデルの名前
    DF自由度、つまりモデル内の自由パラメーターの数
    LogLikモデルの最大化された対数尤度
    LRStataltlmelme を比較するための尤度比検定統計
    pValue尤度比検定の p 値
    LowerpValue の信頼区間の下限
    UpperpValue の信頼区間の上限

シミュレーション出力。以下のフィールドを含む構造体として返します。

nsimnsim の値セット。
alpha'Alpha' の値セット。
pValueSimシミュレーション ベースの p 値。
pValueSimCIpValueSim の信頼区間。ベクトルの最初の要素は下限であり、ベクトルの 2 番目の要素は上限です。
deltaDFaltlme の自由パラメーターの数から lme の自由パラメーターの数を引いた値。altlmeDF から lmeDF を引いた値。
THOモデル lme によって観測済みの応答ベクトル y が生成されたという帰無仮説の下でシミュレーションされた尤度比検定統計のベクトル。

詳細

すべて折りたたむ

尤度比検定

帰無仮説 H0 の場合、観測される尤度比検定統計は、自由度が deltaDF の近似カイ二乗参照分布になります。2 つのモデルを比較するとき、compare は、観測される尤度比検定統計とこのカイ二乗参照分布を比較して、尤度比検定の p 値を計算します。

尤度比検定を使用して取得される p 値は、変量効果の項の有無について検定する場合は保守的になり、固定効果の項の有無について検定する場合は保守的でなくなる可能性もあります。したがって、固定効果の検定のときは、fixedEffectsanova または coefTest メソッドか、シミュレーションされた尤度比検定を使用します。

シミュレーションされた尤度比検定

シミュレーションされた尤度比検定を実行するため、compare は最初に帰無仮説での尤度比検定統計の参照分布を生成します。次に、この参照分布を観測された尤度比検定統計と比較することで、代替モデルの統計的な有意性を評価します。

compare は、次のようにして、帰無仮説での尤度比検定統計のシミュレーションされた参照分布を生成します。

  • 当てはめられたモデル lme から乱数データ ysim を生成します。

  • lme および代替モデル altlme で指定されているモデルをシミュレーションされたデータ ysim に当てはめます。

  • 手順 2 の結果を使用して尤度比検定統計を計算し、値を格納します。

  • 手順 1 から 3 を nsim 回繰り返します。

次に、compare は、観測された尤度比検定統計とシミュレーションされた参照分布を比較して、シミュレーションされた尤度比検定の p 値を計算します。p 値の推定は、シミュレーションされた尤度比検定統計が観測された値プラス 1 以上である回数と、反復回数プラス 1 の比です。

観測された尤度比検定統計が T で、シミュレーションされた参照分布がベクトル TH0 に格納されているとします。このとき次のようになります。

pvalue=[j=1nsimI(TH0(j)T)]+1nsim+1.

シミュレーションされた参照分布の不確実性を考慮するため、compare は true の p 値の 100*(1 – α)% 信頼区間を計算します。

シミュレーションされた尤度比検定を使用して、任意の線形混合効果モデルを比較できます。つまり、シミュレーションされた尤度比検定を使用する場合は、lmealtlme 内に入れ子にする必要はなく、ML (最尤) 法または REML (制限付き最尤) 法を使用して lme および altlme を当てはめることができます。REML (制限付き最尤) 法を使用してモデルを当てはめる場合は、両方のモデルが同じ固定効果計画行列をもっていなければなりません。

入れ子要件

'CheckNesting','True' 名前と値のペアの引数は、以下の要件を確認します。

シミュレーションされた尤度比検定の場合:

  • 同じ方法を使用して両方のモデル (lmealtlme) を当てはめなければなりません。compare では、ML を使用するモデル近似と REML を使用するモデル近似を比較することはできません。

  • 両方のモデルを同じ応答ベクトルに当てはめなければなりません。

  • REML を使用して lmealtlme を当てはめる場合は、両方のモデルが同じ固定効果計画行列をもっていなければなりません。

  • 大きい方のモデル (altlme) の最大化された対数尤度または制限付き対数尤度は、小さい方のモデル (lme) の最大化された対数尤度または制限付き対数尤度より大きいか等しくなければなりません。

理論的検定の場合、'CheckNesting','True' はシミュレーションされた尤度比検定の要件と以下の要件をすべて確認します。

  • lme および altlme の当てはめに使用する重みベクトルは、同一でなければなりません。

  • ML を使用して lmealtlme を当てはめる場合、大きいモデル (altlme) の固定効果の計画行列には、小さいモデル (lme) の固定効果の計画行列が含まれていなければなりません。

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

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

AIC (赤池情報量基準) は、AIC = –2*logLM + 2*(nc + p + 1) です。logLM はモデルの最大化された対数尤度 (または最大化された制限付き対数尤度) であり、nc + p + 1 はモデルで推定されるパラメーターの数です。p は固定効果係数の数であり、nc は残差分散を除く変量効果の共分散のパラメーターの総数です。

ベイズ情報量基準 (BIC) は、BIC = –2*logLM + ln(neff)*(nc + p + 1) です。logLM はモデルの最大化された対数尤度 (つまり最大化された制限付き対数尤度)、neff は有効な観測値数、(nc + p + 1) はモデルで推定されているパラメーターの数です。

  • 近似法が最尤 (ML) 法の場合は neff = n で、n は観測の数です。

  • 近似法が制限付き最尤 (REML) 法の場合は neff = n-p です。

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

逸脱度

LinearMixedModel は、モデル M の逸脱度を、そのモデルの対数尤度のマイナス 2 倍として計算します。LM がモデル M の尤度関数の最大値を示しているものとします。この場合、モデル M の逸脱度は次のようになります。

2*logLM.

逸脱度の値が小さいほど、近似が優れていることを意味します。M1 と M2 は 2 つの異なるモデルであり、M1 は M2 に入れ子になっているものとします。この場合、モデルの当てはめはこれらのモデルの逸脱度 Dev1 と Dev2 を比較することによって評価できます。逸脱度の差異は以下のとおりです。

Dev=Dev1Dev2=2(logLM2logLM1).

通常、この差異の漸近分布はカイ二乗分布であり、その自由度 v は、1 つのモデルで推定され、他のモデルでは固定されている (通常は 0) パラメーターの数に等しくなります。つまり、M1 と M2 で推定されるパラメーターの数の差異と等しくなります。この検定の p 値は 1 – chi2cdf(Dev,V) を使用して得ることができます。ここで、Dev = Dev2 - Dev1 です。

ただし、混合効果モデルでは、一部の分散成分がパラメーター空間の境界上になると、この差異の漸近分布はさらに複雑になります。たとえば、次の仮説について考えます。

H0: D=(D11000), とした場合、D は q 行 q 列の対称な半正定値行列である。

H1: D は、(q+1) 行 (q+1) 列の対称な半正定値行列です。

つまり、H1 は、D の最後の行と列がゼロとは異なることを記述します。ここで、大きい方のモデル M2 には q + 1 個のパラメーターがあり、小さい方のモデル M1 には q 個のパラメーターがあります。また、Dev では、χ2q 分布と χ2(q + 1) 分布が 50:50 で混合されています (Stram and Lee、1994)。

参照

[1] Hox, J. Multilevel Analysis, Techniques and Applications. Lawrence Erlbaum Associates, Inc., 2002.

[2] Stram D. O. and J. W. Lee. “Variance components testing in the longitudinal mixed-effects model”. Biometrics, Vol. 50, 4, 1994, pp. 1171–1177.

拡張機能

バージョン履歴

R2013b で導入