compare
線形混合効果モデルの比較
構文
説明
は、線形混合効果モデル 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(___,Name,Value
)lme
と altlme
を比較するシミュレーションされた尤度比検定の結果を返しますが、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 つのモデルを比較します。また、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 という小さい 値は、モデル 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 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]
大きい 値は、大きいモデル 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.095904 0.078373 Upper 0.11585
大きい 値は、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
大きい 値は、大きいモデル lme
が小さいモデル lme2
より有意には優れていないことを示します。lme2
に対する AIC
と BIC
の値が小さいこともこれを裏付けています。
入力引数
lme
— 線形混合効果モデル
LinearMixedModel
オブジェクト
線形混合効果モデル。fitlme
または fitlmematrix
を使用して構築した LinearMixedModel
オブジェクトとして指定します。
altlme
— 代替線形混合効果モデル
LinearMixedModel
オブジェクト
同じ応答ベクトルに異なるモデル仕様によって当てはめる代替線形混合効果モデル。LinearMixedModel
オブジェクトとして指定されます。lme
は altlme
で入れ子にされている必要があります。つまり、一部のパラメーターを 0 などの固定値に設定することにより lme
を altlme
から取得する必要があります。線形混合効果オブジェクトは、fitlme
または fitlmematrix
を使用して作成できます。
nsim
— シミュレーションの反復数
正の整数
シミュレーションされた尤度比検定でのシミュレーションの反復数。正の整数値として指定します。シミュレーションされた尤度比検定を実行するには nsim
を指定しなければなりません。
例: 'NSim',1000
データ型: double
| single
名前と値の引数
オプションの引数のペアを Name1=Value1,...,NameN=ValueN
として指定します。ここで Name
は引数名、Value
は対応する値です。名前と値の引数は他の引数の後ろにする必要がありますが、ペアの順序は関係ありません。
R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name
を引用符で囲みます。
例: results = compare(lme,altlme,'CheckNesting',true)
Alpha
— 有意水準
0.05 (既定値) | 0 ~ 1 の範囲のスカラー値
有意水準。'Alpha'
と 0 ~ 1 の範囲にあるスカラー値から構成されるコンマ区切りのペアとして指定します。値が α の場合、信頼水準は 100 × (1 – α)% です。
たとえば、99% の信頼区間の場合は、次のように信頼水準を指定できます。
例: 'Alpha',0.01
データ型: single
| double
Options
— シミュレーションされた尤度比検定を実行するためのオプション
構造体
シミュレーションされた尤度比検定を並列実行するためのオプション。'Options'
と、statset('LinearMixedModel')
によって作成される構造体から構成されるコンマ区切りのペアとして指定します。
これらのオプションには Parallel Computing Toolbox™ が必要です。
compare
は次のフィールドを使用します。
'UseParallel' |
並列計算には Parallel Computing Toolbox が必要です。 |
'UseSubstreams' |
|
'Streams' |
|
コマンド ラインでの並列統計計算の詳細を表示するには、次のように入力します。
help parallelstats
データ型: struct
CheckNesting
— 2 つのモデルの間の入れ子を調べるためのインジケーター
false
(既定値) | true
2 つのモデルの間の入れ子を確認するためのインジケーター。'CheckNesting'
と以下のいずれかで構成されるコンマ区切りのペアとして指定します。
false | 既定の設定。確認しません。 |
true | compare は、小さい方のモデル lme が大きい方のモデル altlme に入れ子になっているかどうかを確認します。 |
理論的尤度比検定を有効にするためには、lme
を代替モデル altlme
で入れ子にしなければなりません。入れ子要件が満たされていない場合、compare
はエラー メッセージを返します。
どちらの検定に対しても有効ですが、シミュレーションされた尤度比検定の方が入れ子の要件は弱くなります。
例: 'CheckNesting',true
データ型: single
| double
出力引数
results
— 尤度比検定またはシミュレーションされた尤度比検定の結果
データセット配列
尤度比検定またはシミュレーションされた尤度比検定の結果。2 行のデータセット配列として返します。1 行目は lme
に対応し、2 行目は altlme
に対応しています。results
の列は、検定が尤度比検定かシミュレーションされた尤度比検定かに依存します。
尤度比検定を使用している場合は、
results
に以下の列が含まれます。Model
モデルの名前 DF
自由度、つまりモデル内の自由パラメーターの数 AIC
モデルの赤池情報量基準 BIC
モデルのベイズ情報量基準 LogLik
モデルの最大化された対数尤度 LRStat
altlme
とlme
を比較するための尤度比検定統計deltaDF
altlme
のDF
からlme
のDF
を引いた値pValue
尤度比検定の p 値 シミュレーションされた尤度比検定を使用している場合は、
results
に以下の列が含まれます。Model
モデルの名前 DF
自由度、つまりモデル内の自由パラメーターの数 LogLik
モデルの最大化された対数尤度 LRStat
altlme
とlme
を比較するための尤度比検定統計pValue
尤度比検定の p 値 Lower
pValue
の信頼区間の下限Upper
pValue
の信頼区間の上限
siminfo
— シミュレーション出力
構造体
シミュレーション出力。以下のフィールドを含む構造体として返します。
nsim | nsim の値セット。 |
alpha | 'Alpha' の値セット。 |
pValueSim | シミュレーション ベースの p 値。 |
pValueSimCI | pValueSim の信頼区間。ベクトルの最初の要素は下限であり、ベクトルの 2 番目の要素は上限です。 |
deltaDF | altlme の自由パラメーターの数から lme の自由パラメーターの数を引いた値。altlme の DF から lme の DF を引いた値。 |
THO | モデル lme によって観測済みの応答ベクトル y が生成されたという帰無仮説の下でシミュレーションされた尤度比検定統計のベクトル。 |
詳細
尤度比検定
帰無仮説 H0 の場合、観測される尤度比検定統計は、自由度が deltaDF
の近似カイ二乗参照分布になります。2 つのモデルを比較するとき、compare
は、観測される尤度比検定統計とこのカイ二乗参照分布を比較して、尤度比検定の p 値を計算します。
尤度比検定を使用して取得される p 値は、変量効果の項の有無について検定する場合は保守的になり、固定効果の項の有無について検定する場合は保守的でなくなる可能性もあります。したがって、固定効果の検定のときは、fixedEffects
、anova
または coefTest
メソッドか、シミュレーションされた尤度比検定を使用します。
シミュレーションされた尤度比検定
シミュレーションされた尤度比検定を実行するため、compare
は最初に帰無仮説での尤度比検定統計の参照分布を生成します。次に、この参照分布を観測された尤度比検定統計と比較することで、代替モデルの統計的な有意性を評価します。
compare
は、次のようにして、帰無仮説での尤度比検定統計のシミュレーションされた参照分布を生成します。
当てはめられたモデル
lme
から乱数データysim
を生成します。lme
および代替モデルaltlme
で指定されているモデルをシミュレーションされたデータysim
に当てはめます。手順 2 の結果を使用して尤度比検定統計を計算し、値を格納します。
手順 1 から 3 を
nsim
回繰り返します。
次に、compare
は、観測された尤度比検定統計とシミュレーションされた参照分布を比較して、シミュレーションされた尤度比検定の p 値を計算します。p 値の推定は、シミュレーションされた尤度比検定統計が観測された値プラス 1 以上である回数と、反復回数プラス 1 の比です。
観測された尤度比検定統計が T で、シミュレーションされた参照分布がベクトル TH0 に格納されているとします。このとき次のようになります。
シミュレーションされた参照分布の不確実性を考慮するため、compare
は true の p 値の 100*(1 – α)% 信頼区間を計算します。
シミュレーションされた尤度比検定を使用して、任意の線形混合効果モデルを比較できます。つまり、シミュレーションされた尤度比検定を使用する場合は、lme
を altlme
内に入れ子にする必要はなく、ML (最尤) 法または REML (制限付き最尤) 法を使用して lme
および altlme
を当てはめることができます。REML (制限付き最尤) 法を使用してモデルを当てはめる場合は、両方のモデルが同じ固定効果計画行列をもっていなければなりません。
入れ子要件
'CheckNesting','True'
名前と値のペアの引数は、以下の要件を確認します。
シミュレーションされた尤度比検定の場合:
同じ方法を使用して両方のモデル (
lme
とaltlme
) を当てはめなければなりません。compare
では、ML を使用するモデル近似と REML を使用するモデル近似を比較することはできません。両方のモデルを同じ応答ベクトルに当てはめなければなりません。
REML を使用して
lme
とaltlme
を当てはめる場合は、両方のモデルが同じ固定効果計画行列をもっていなければなりません。大きい方のモデル (
altlme
) の最大化された対数尤度または制限付き対数尤度は、小さい方のモデル (lme
) の最大化された対数尤度または制限付き対数尤度より大きいか等しくなければなりません。
理論的検定の場合、'CheckNesting','True'
はシミュレーションされた尤度比検定の要件と以下の要件をすべて確認します。
lme
およびaltlme
の当てはめに使用する重みベクトルは、同一でなければなりません。ML を使用して
lme
とaltlme
を当てはめる場合、大きいモデル (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 の逸脱度は次のようになります。
逸脱度の値が小さいほど、近似が優れていることを意味します。M1 と M2 は 2 つの異なるモデルであり、M1 は M2 に入れ子になっているものとします。この場合、モデルの当てはめはこれらのモデルの逸脱度 Dev1 と Dev2 を比較することによって評価できます。逸脱度の差異は以下のとおりです。
通常、この差異の漸近分布はカイ二乗分布であり、その自由度 v は、1 つのモデルで推定され、他のモデルでは固定されている (通常は 0) パラメーターの数に等しくなります。つまり、M1 と M2 で推定されるパラメーターの数の差異と等しくなります。この検定の p 値は 1 – chi2cdf(Dev,V)
を使用して得ることができます。ここで、Dev = Dev2 - Dev1 です。
ただし、混合効果モデルでは、一部の分散成分がパラメーター空間の境界上になると、この差異の漸近分布はさらに複雑になります。たとえば、次の仮説について考えます。
H0: とした場合、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.
拡張機能
自動並列サポート
Parallel Computing Toolbox™ を使用して自動的に並列計算を実行することで、コードを高速化します。
並列実行するには、この関数を呼び出すときに名前と値の引数 Options
を指定し、statset
を使用してオプション構造体の UseParallel
フィールドを true
に設定します。
Options=statset(UseParallel=true)
並列計算の詳細については、自動並列サポートを使用した MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。
バージョン履歴
R2013b で導入
MATLAB コマンド
次の MATLAB コマンドに対応するリンクがクリックされました。
コマンドを MATLAB コマンド ウィンドウに入力して実行してください。Web ブラウザーは MATLAB コマンドをサポートしていません。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)