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 値の信頼区間の信頼水準を変更したりできます。
例
標本データを読み込み、table の形式に変換します。
load flu
flu = dataset2table(flu)flu=52×11 table
Date NE MidAtl ENCentral WNCentral SAtl ESCentral WSCentral Mtn Pac WtdILI
______________ _____ ______ _________ _________ _____ _________ _________ _____ _____ ______
{'10/9/2005' } 0.97 1.025 1.232 1.286 1.082 1.457 1.1 0.981 0.971 1.182
{'10/16/2005'} 1.136 1.06 1.228 1.286 1.146 1.644 1.123 0.976 0.917 1.22
{'10/23/2005'} 1.135 1.172 1.278 1.536 1.274 1.556 1.236 1.102 0.895 1.31
{'10/30/2005'} 1.52 1.489 1.576 1.794 1.59 2.252 1.612 1.321 1.082 1.343
{'11/6/2005' } 1.365 1.394 1.53 1.825 1.62 2.059 1.471 1.453 1.118 1.586
{'11/13/2005'} 1.39 1.477 1.506 1.9 1.683 1.813 1.464 1.388 1.204 1.47
{'11/20/2005'} 1.212 1.231 1.295 1.495 1.347 1.794 1.303 1.371 1.137 1.611
{'11/27/2005'} 1.477 1.546 1.557 1.855 1.678 2.159 1.739 1.628 1.443 1.827
{'12/4/2005' } 1.285 1.43 1.482 1.635 1.577 1.903 1.53 1.701 1.516 1.776
{'12/11/2005'} 1.354 1.45 1.46 1.794 1.583 1.894 1.831 2.364 2.094 1.941
{'12/18/2005'} 1.502 1.622 1.638 1.988 1.947 2.22 2.577 3.89 2.66 2.34
{'12/25/2005'} 1.86 1.915 1.955 2.38 2.343 3.027 3.219 4.862 2.595 3.086
{'1/1/2006' } 2.114 2.174 2.065 2.557 2.275 2.498 2.644 3.352 2.181 3.26
{'1/8/2006' } 1.815 1.932 1.822 2.046 1.969 1.805 2.189 2.132 1.717 2.613
{'1/15/2006' } 1.541 1.695 1.581 2.008 1.718 1.662 2.156 1.694 1.351 2.247
{'1/22/2006' } 1.632 1.758 1.711 2.217 1.866 2.194 2.268 1.826 1.384 2.352
⋮
table flu には、変数 Date と、インフルエンザ推定罹患率 (Google® 検索から推定される 9 地域の値と CDC による全国の推定値) が格納されている 10 個の変数が含まれています。
線形混合効果モデルを当てはめるには、データが適切な形式の table になっていなければなりません。インフルエンザ罹患率を応答として、地域を予測子変数として線形混合効果モデルを当てはめるため、地域に対応する 9 個の列を 1 つの配列にまとめます。新しい table flu2 には、応答変数 FluRate、各推定の元になっている地域を示すノミナル変数 Region、およびグループ化変数 Date が含まれなければなりません。
flu2 = stack(flu,2:10,NewDataVariableName="FluRate",... IndexVariableName="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;
tbltbl=60×4 table
Soil Tomato Fertilizer Yield
_________ ____________ __________ _____
{'Sandy'} {'Plum' } 1 104
{'Sandy'} {'Plum' } 2 136
{'Sandy'} {'Plum' } 3 158
{'Sandy'} {'Plum' } 4 174
{'Sandy'} {'Cherry' } 1 57
{'Sandy'} {'Cherry' } 2 86
{'Sandy'} {'Cherry' } 3 89
{'Sandy'} {'Cherry' } 4 98
{'Sandy'} {'Heirloom'} 1 65
{'Sandy'} {'Heirloom'} 2 62
{'Sandy'} {'Heirloom'} 3 113
{'Sandy'} {'Heirloom'} 4 84
{'Sandy'} {'Grape' } 1 54
{'Sandy'} {'Grape' } 2 86
{'Sandy'} {'Grape' } 3 89
{'Sandy'} {'Grape' } 4 115
⋮
table tbl には土壌の種類に基づいて土壌が 3 つのブロックに分けられている分割プロット試験のデータが含まれています。土壌の種類は砂質、シルト、および粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。その後、プロット内のトマトの苗木はサブプロットに分割され、それぞれのサブプロットが 4 つの肥料の中の 1 つにより処置されます。このデータは、シミュレーションされたものです。
Tomato、Soil、および Fertilizer をカテゴリカル変数に変換します。
tbl.Tomato = nominal(tbl.Tomato); tbl.Soil = nominal(tbl.Soil); tbl.Fertilizer = nominal(tbl.Fertilizer);
線形混合効果モデルを当てはめます。Fertilizer および Tomato は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。
lmeBig = fitlme(tbl,"Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)");交互作用項 Tomato:Fertilizer および変量効果の項 (1 | Soil) を削除した後、モデルを再度当てはめます。
lmeSmall = fitlme(tbl,"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: [1000×1 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
tbltbl=60×4 table
Soil Tomato Fertilizer Yield
_________ ____________ __________ _____
{'Sandy'} {'Plum' } 1 104
{'Sandy'} {'Plum' } 2 136
{'Sandy'} {'Plum' } 3 158
{'Sandy'} {'Plum' } 4 174
{'Sandy'} {'Cherry' } 1 57
{'Sandy'} {'Cherry' } 2 86
{'Sandy'} {'Cherry' } 3 89
{'Sandy'} {'Cherry' } 4 98
{'Sandy'} {'Heirloom'} 1 65
{'Sandy'} {'Heirloom'} 2 62
{'Sandy'} {'Heirloom'} 3 113
{'Sandy'} {'Heirloom'} 4 84
{'Sandy'} {'Grape' } 1 54
{'Sandy'} {'Grape' } 2 86
{'Sandy'} {'Grape' } 3 89
{'Sandy'} {'Grape' } 4 115
⋮
table tbl には土壌の種類に基づいて土壌が 3 つのブロックに分けられている分割プロット試験のデータが含まれています。土壌の種類は砂質、シルト、および粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。その後、プロット内のトマトの苗木はサブプロットに分割され、それぞれのサブプロットが 4 つの肥料の中の 1 つにより処置されます。このデータは、シミュレーションされたものです。
Tomato、Soil、および 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 'Processes' profile ... 13-Nov-2024 10:18:44: Job Queued. Waiting for parallel pool job with ID 4 to start ... Connected to parallel pool with 4 workers.
反復数 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.52747 0.496 0.55878
大きい 値は、大きいモデル 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.
拡張機能
バージョン履歴
R2013b で導入
