Main Content

compare

クラス: LinearMixedModel

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

このページをすべて展開する

構文

results = compare(lme,altlme)
results = compare(___,Name,Value)
[results,siminfo] = compare(lme,altlme,'NSim',nsim)
[results,siminfo] = compare(___,Name,Value)

説明

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

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

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

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

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

固定効果を検定するには、lme と altlme を ML で近似してシミュレーションされた尤度比検定で compare を使用するか、fixedEffects、anova または coefTest メソッドを使用します。

例

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

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

例

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

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

例

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

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

入力引数

すべて展開する

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

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

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

例: 'NSim',1000

データ型: double | single

名前と値の引数

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

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

有意水準。'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モデルの最大化された対数尤度
    LRStataltlme と lme を比較するための尤度比検定統計
    deltaDFaltlme の DF から lme の DF を引いた値
    pValue尤度比検定の p 値

  • シミュレーションされた尤度比検定を使用している場合は、results に以下の列が含まれます。

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

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

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

例

すべて展開する

ライブ スクリプトを開く

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

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 つのモデルを比較します。また、lme2 が lme の入れ子になっているかどうかを調べます。 

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 という名前のデータセット配列に保存し、Tomato、Soil および 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 
    lmeSmall    10    511.06       532    -245.53                     
    lmeBig      23    522.57    570.74    -238.29    14.491    0.57343


    Lower      Upper  
                      
    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'});

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

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.094905    0.077462


    Upper  
           
    0.11477

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

この例では次を使用します:

ライブ スクリプトを開く

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

load('fertilizer.mat')

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

このデータを tbl という名前のテーブルに保存し、Tomato、Soil および 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 の並列計算でシミュレーションされた尤度比検定を使用し、lme2 と lme を比較します。 

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 に対する AIC と BIC の値が小さいこともこれを裏付けています。

詳細

すべて展開する

参照

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

拡張機能

参考

| | | | | |