このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。
fitlme
線形混合効果モデルの当てはめ
説明
例
線形混合効果モデルの当てはめ
標本データを読み込みます。
load imports-85
変数をテーブルに格納します。
tbl = table(X(:,12),X(:,14),X(:,24),'VariableNames',{'Horsepower','CityMPG','EngineType'});
テーブルの最初の 5 行を表示します。
tbl(1:5,:)
ans=5×3 table
Horsepower CityMPG EngineType
__________ _______ __________
111 21 13
111 21 13
154 19 37
102 24 35
115 18 35
市内でのガロンあたりの走行マイル数について、線形混合効果モデルを近似します。ここでは、馬力の固定効果とエンジンの種類でグループ化した切片と馬力の無相関な変量効果を使用します。
lme = fitlme(tbl,'CityMPG~Horsepower+(1|EngineType)+(Horsepower-1|EngineType)');
このモデルでは、CityMPG
は応答変数、馬力は予測子変数、エンジンの種類はグループ化変数です。既定の設定では切片が含まれているため、モデルの固定効果部分は 1 + Horsepower
に対応します。
切片と馬力の変量効果項は無相関であるため、これらの項は個別に指定されます。2 番目の変量効果項は馬力しか表していないため、–1
を使用して 2 番目の変量効果項から切片を消去しなければなりません。
モデルを表示します。
lme
lme = Linear mixed-effects model fit by ML Model information: Number of observations 203 Fixed effects coefficients 2 Random effects coefficients 14 Covariance parameters 3 Formula: CityMPG ~ 1 + Horsepower + (1 | EngineType) + (Horsepower | EngineType) Model fit statistics: AIC BIC LogLikelihood Deviance 1099.5 1116 -544.73 1089.5 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 37.276 2.8556 13.054 201 1.3147e-28 {'Horsepower' } -0.12631 0.02284 -5.53 201 9.8848e-08 Lower Upper 31.645 42.906 -0.17134 -0.081269 Random effects covariance parameters (95% CIs): Group: EngineType (7 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 5.7338 Lower Upper 2.3773 13.829 Group: EngineType (7 Levels) Name1 Name2 Type Estimate {'Horsepower'} {'Horsepower'} {'std'} 0.050357 Lower Upper 0.02307 0.10992 Group: Error Name Estimate Lower Upper {'Res Std'} 3.226 2.9078 3.5789
切片と馬力の変量効果共分散パラメーターが別個に表示されています。
市内でのガロンあたりの走行マイル数について、線形混合効果モデルを近似します。今回は、エンジンの種類でグループ化した切片と馬力の同じ固定効果項および相関が考えられる変量効果を使用します。
lme2 = fitlme(tbl,'CityMPG~Horsepower+(Horsepower|EngineType)');
既定の設定では変量効果項に切片が含まれているため、1
を追加する必要はありません。変量効果項は (1 + Horsepower|EngineType)
と等しくなります。
モデルを表示します。
lme2
lme2 = Linear mixed-effects model fit by ML Model information: Number of observations 203 Fixed effects coefficients 2 Random effects coefficients 14 Covariance parameters 4 Formula: CityMPG ~ 1 + Horsepower + (1 + Horsepower | EngineType) Model fit statistics: AIC BIC LogLikelihood Deviance 1089 1108.9 -538.52 1077 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 33.824 4.0181 8.4178 201 7.1678e-15 {'Horsepower' } -0.1087 0.032912 -3.3029 201 0.0011328 Lower Upper 25.901 41.747 -0.1736 -0.043806 Random effects covariance parameters (95% CIs): Group: EngineType (7 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 9.4952 {'Horsepower' } {'(Intercept)'} {'corr'} -0.96843 {'Horsepower' } {'Horsepower' } {'std' } 0.078874 Lower Upper 4.7022 19.174 -0.99568 -0.78738 0.039917 0.15585 Group: Error Name Estimate Lower Upper {'Res Std'} 3.1845 2.8774 3.5243
切片と馬力の変量効果共分散パラメーターが同時に表示され、さらに切片と馬力の間の相関 ('corr'
) も含まれています。
切片がランダムな LME モデルの当てはめ
標本データを読み込みます。
load flu
データセット配列 flu
には、変数 Date
と、インフルエンザ推定罹患率 (Google® 検索から推定される 9 地域の値と疾病対策センター (CDC) による全国の推定値) が格納されている 10 個の変数が含まれています。
線形混合効果モデルを近似するには、データが適切な形式のデータセット配列になっていなければなりません。インフルエンザ羅患率を応答として使用することにより線形混合効果モデルを当てはめるには、地域に対応する 9 個の列を 1 つの配列にまとめます。新しいデータセット配列 flu2
には、新しい応答変数 FluRate
、各推定の元になっている地域を示すノミナル変数 Region
、全国の推定値 WtdILI
およびグループ化変数 Date
が含まれなければなりません。
flu2 = stack(flu,2:10,'NewDataVarName','FluRate', ... 'IndVarName','Region'); flu2.Date = nominal(flu2.Date);
flu2
の最初の 6 行を表示します。
flu2(1:6,:)
ans = Date WtdILI Region FluRate 10/9/2005 1.182 NE 0.97 10/9/2005 1.182 MidAtl 1.025 10/9/2005 1.182 ENCentral 1.232 10/9/2005 1.182 WNCentral 1.286 10/9/2005 1.182 SAtl 1.082 10/9/2005 1.182 ESCentral 1.457
全国の推定値 WtdILI
に対する固定効果の項と、Date
で変化するランダム切片で、線形混合効果モデルを近似します。このモデルは以下に対応します。
ここで、 はグループ化変数 Date
のレベル に対する観測値 、 はグループ化変数 Date
のレベル に対する変量効果、 は観測値 の観測誤差です。変量効果の事前分布は次のようになります。
誤差項の分布は次のようになります。
lme = fitlme(flu2,'FluRate ~ 1 + WtdILI + (1|Date)')
lme = Linear mixed-effects model fit by ML Model information: Number of observations 468 Fixed effects coefficients 2 Random effects coefficients 52 Covariance parameters 2 Formula: FluRate ~ 1 + WtdILI + (1 | Date) Model fit statistics: AIC BIC LogLikelihood Deviance 286.24 302.83 -139.12 278.24 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)'} 0.16385 0.057525 2.8484 466 0.0045885 {'WtdILI' } 0.7236 0.032219 22.459 466 3.0502e-76 Lower Upper 0.050813 0.27689 0.66028 0.78691 Random effects covariance parameters (95% CIs): Group: Date (52 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 0.17146 Lower Upper 0.13227 0.22226 Group: Error Name Estimate Lower Upper {'Res Std'} 0.30201 0.28217 0.32324
"Random effects covariance parameters" というタイトルのセクションに、推定共分散パラメーターが表示されます。 の推定値は 0.17146、95% 信頼区間は [0.13227, 0.22226] です。この間隔には 0 は含まれていないため、変量効果項は有意です。compare
メソッドによる尤度比検定を使用すると、任意の変量効果項の有意性を正式に検定できます。
観測値について推定される応答は、その観測値に対応するグループ化変数レベルの固定効果および変量効果の値の合計です。たとえば、観測値 28 のインフルエンザ推定羅患率は次のようになります。
ここで、 は切片の変量効果について推定した最良線形不偏予測量 (BLUP) です。この値を以下のように計算することもできます。
beta = fixedEffects(lme); [~,~,STATS] = randomEffects(lme); % Compute the random-effects statistics (STATS) STATS.Level = nominal(STATS.Level); y_hat = beta(1) + beta(2)*flu2.WtdILI(28) + STATS.Estimate(STATS.Level=='10/30/2005')
y_hat = 1.4674
fitted
メソッドを使用して、近似値を表示できます。
F = fitted(lme); F(28)
ans = 1.4674
乱塊法用の LME モデル
標本データを読み込みます。
load('shift.mat')
このデータは 5 人の作業者が 3 つのシフトの間に製造した製品から計測された品質目標の特性の絶対偏差を示します。3 つのシフトとは朝、夕方、夜です。これは作業者をブロックとする乱塊法です。この実験は、シフトの時間によるパフォーマンスへの影響の調査を意図しています。パフォーマンスの測定は、目標値からの品質特性の絶対偏差です。このデータは、シミュレーションされたものです。
シフトの時間によってパフォーマンスに有意差があるかどうかを評価するために、作業者別のランダムな切片をもつ線形混合効果モデルに当てはめます。制限付き最尤法と 'effects'
対比を使用します。
'effects'
対比は係数の合計が 0 であることを意味し、fitlme
はシフトの効果を表す、"固定効果の計画行列" と呼ばれる行列を作成します。この行列には、 および という 2 つの列があります。ここで、
このモデルは以下に対応します。
ここで、 は観測値を、 は作業者を表します。 = 1、2、...、15 および = 1、2、...、5 です。変量効果と観測値の誤差の分布は次のとおりです。
および
lme = fitlme(shift,'QCDev ~ Shift + (1|Operator)',... 'FitMethod','REML','DummyVarCoding','effects')
lme = Linear mixed-effects model fit by REML Model information: Number of observations 15 Fixed effects coefficients 3 Random effects coefficients 5 Covariance parameters 2 Formula: QCDev ~ 1 + Shift + (1 | Operator) Model fit statistics: AIC BIC LogLikelihood Deviance 58.913 61.337 -24.456 48.913 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF pValue {'(Intercept)' } 3.6525 0.94109 3.8812 12 0.0021832 {'Shift_Evening'} -0.53293 0.31206 -1.7078 12 0.11339 {'Shift_Morning'} -0.91973 0.31206 -2.9473 12 0.012206 Lower Upper 1.6021 5.703 -1.2129 0.14699 -1.5997 -0.23981 Random effects covariance parameters (95% CIs): Group: Operator (5 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 2.0457 Lower Upper 0.98207 4.2612 Group: Error Name Estimate Lower Upper {'Res Std'} 0.85462 0.52357 1.395
変量効果の最良線形不偏予測量 (BLUP) の推定値を計算します。
B = randomEffects(lme)
B = 5×1
0.5775
1.1757
-2.1715
2.3655
-1.9472
夕方シフトで作業している 3 番目の作業者の場合、品質目標の特性の推定絶対偏差は以下のとおりです。
この値を以下のように表示することもできます。
F = fitted(lme); F(shift.Shift=='Evening' & shift.Operator=='3')
ans = 0.9481
同様に、朝シフトで作業している 3 番目の作業者の場合、品質目標の特性の推定絶対偏差を以下のように計算できます。
この値を以下のように表示することもできます。
F(shift.Shift=='Morning' & shift.Operator=='3')
ans = 0.5613
この作業者は、朝シフトのときに誤差の大きさが減少する傾向にあります。
分割プロット試験用の 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
は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。
このモデルは以下に対応します。
ここで = 1、2、...、60 で、インデックス は肥料の種類、 はトマトの種類に対応し、 = 1、2、3 はブロック (土壌) に対応します。 は 番目の土壌の種類を表し、 は 番目の土壌の種類で入れ子にされている 番目のトマトの種類を表します。 は肥料のレベル を表すダミー変数です。同様に、 は、トマトの種類のレベル を表すダミー変数です。
変量効果と観測誤差の事前分布は、~N(0, )、~N(0, )、および ~ N(0, ) です。
lme = fitlme(ds,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)')
lme = Linear mixed-effects model fit by ML Model information: Number of observations 60 Fixed effects coefficients 20 Random effects coefficients 18 Covariance parameters 3 Formula: Yield ~ 1 + Tomato*Fertilizer + (1 | Soil) + (1 | Soil:Tomato) Model fit statistics: AIC BIC LogLikelihood Deviance 522.57 570.74 -238.29 476.57 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF {'(Intercept)' } 77 8.5836 8.9706 40 {'Tomato_Grape' } -16 11.966 -1.3371 40 {'Tomato_Heirloom' } -6.6667 11.966 -0.55714 40 {'Tomato_Plum' } 32.333 11.966 2.7022 40 {'Tomato_Vine' } -13 11.966 -1.0864 40 {'Fertilizer_2' } 34.667 8.572 4.0442 40 {'Fertilizer_3' } 33.667 8.572 3.9275 40 {'Fertilizer_4' } 47.667 8.572 5.5607 40 {'Tomato_Grape:Fertilizer_2' } -2.6667 12.123 -0.21997 40 {'Tomato_Heirloom:Fertilizer_2'} -8 12.123 -0.65992 40 {'Tomato_Plum:Fertilizer_2' } -15 12.123 -1.2374 40 {'Tomato_Vine:Fertilizer_2' } -16 12.123 -1.3198 40 {'Tomato_Grape:Fertilizer_3' } 16.667 12.123 1.3748 40 {'Tomato_Heirloom:Fertilizer_3'} 3.3333 12.123 0.27497 40 {'Tomato_Plum:Fertilizer_3' } 3.6667 12.123 0.30246 40 {'Tomato_Vine:Fertilizer_3' } 3 12.123 0.24747 40 {'Tomato_Grape:Fertilizer_4' } 13.333 12.123 1.0999 40 {'Tomato_Heirloom:Fertilizer_4'} -19 12.123 -1.5673 40 {'Tomato_Plum:Fertilizer_4' } -2.6667 12.123 -0.21997 40 {'Tomato_Vine:Fertilizer_4' } 8.6667 12.123 0.71492 40 pValue Lower Upper 4.0206e-11 59.652 94.348 0.18873 -40.184 8.1837 0.58053 -30.85 17.517 0.010059 8.1496 56.517 0.28379 -37.184 11.184 0.00023272 17.342 51.991 0.00033057 16.342 50.991 1.9567e-06 30.342 64.991 0.82701 -27.167 21.834 0.51309 -32.501 16.501 0.22317 -39.501 9.5007 0.19439 -40.501 8.5007 0.17683 -7.8341 41.167 0.78476 -21.167 27.834 0.76387 -20.834 28.167 0.80581 -21.501 27.501 0.27796 -11.167 37.834 0.12492 -43.501 5.5007 0.82701 -27.167 21.834 0.47881 -15.834 33.167 Random effects covariance parameters (95% CIs): Group: Soil (3 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 2.5028 Lower Upper 0.027711 226.05 Group: Soil:Tomato (15 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 10.225 Lower Upper 6.1497 17.001 Group: Error Name Estimate Lower Upper {'Res Std'} 10.499 8.5389 12.908
固定効果係数の表示の最後の 12 行に対応する 値 (0.82701 ~ 0.47881) は、トマトの種類と肥料の種類の間の交互作用係数が有意ではないことを示しています。トマトと肥料の間の全体的な交互作用を検定するには、'effects'
対比を使用してモデルを再度当てはめた後でanova
メソッドを使用します。
変量効果項の標準偏差 () の信頼区間は非常に大きくなっています。切片は土壌でグループ化されています。この項には有意性がありません。
交互作用項 Tomato:Fertilizer
および変量効果の項 (1 | Soil)
を削除した後、モデルを再近似します。
lme = fitlme(ds,'Yield ~ Fertilizer + Tomato + (1|Soil:Tomato)')
lme = Linear mixed-effects model fit by ML Model information: Number of observations 60 Fixed effects coefficients 8 Random effects coefficients 15 Covariance parameters 2 Formula: Yield ~ 1 + Tomato + Fertilizer + (1 | Soil:Tomato) Model fit statistics: AIC BIC LogLikelihood Deviance 511.06 532 -245.53 491.06 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF {'(Intercept)' } 77.733 7.3293 10.606 52 {'Tomato_Grape' } -9.1667 9.6045 -0.95441 52 {'Tomato_Heirloom'} -12.583 9.6045 -1.3102 52 {'Tomato_Plum' } 28.833 9.6045 3.0021 52 {'Tomato_Vine' } -14.083 9.6045 -1.4663 52 {'Fertilizer_2' } 26.333 4.5004 5.8514 52 {'Fertilizer_3' } 39 4.5004 8.6659 52 {'Fertilizer_4' } 47.733 4.5004 10.607 52 pValue Lower Upper 1.3108e-14 63.026 92.441 0.34429 -28.439 10.106 0.1959 -31.856 6.6895 0.0041138 9.5605 48.106 0.14858 -33.356 5.1895 3.3024e-07 17.303 35.364 1.1459e-11 29.969 48.031 1.308e-14 38.703 56.764 Random effects covariance parameters (95% CIs): Group: Soil:Tomato (15 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std'} 10.02 Lower Upper 6.0812 16.509 Group: Error Name Estimate Lower Upper {'Res Std'} 12.325 10.024 15.153
固定効果項と変量効果項の両方が検定されるので、シミュレートされた尤度比検定でcompare
メソッドを使用して 2 つのモデルを比較できます。
共変量を使用した長期間の調査
標本データを読み込みます。
load('weight.mat')
weight
には長期間の調査によるデータが含まれています。そこには 20 人の被験者が 4 つの運動プログラム (A、B、C、D) にランダムに割り当てられ、体重の減少が 6 回の 2 週間の期間にわたって記録されています。このデータは、シミュレーションされたものです。
データをテーブルに保存します。Subject
および Program
をカテゴリカル変数として定義します。
tbl = table(InitialWeight,Program,Subject,Week,y); tbl.Subject = nominal(tbl.Subject); tbl.Program = nominal(tbl.Program);
線形混合効果モデルを近似します。初期体重、プログラムの種類、週、週とプログラムの種類の間の交互作用は固定効果です。切片と週は被験者ごとに異なります。
fitlme
はプログラム A を基準として使用して、必要なダミー変数 [.] を作成します。モデルには既に切片があるので、fitlme
は、プログラム B、C、D に対してダミー変数のみを作成します。これは、ダミー変数をコーディングする 'reference'
メソッドとも呼ばれます。
このモデルは以下に対応します。
ここで、 = 1、2、...、120、 = 1、2、...、20 です。 は固定効果係数 ( = 0、1、...、8)、 と は変量効果です。 は初期体重を表し、 はプログラムの種類を表すダミー変数です。たとえば、 はプログラムの種類 B を表すダミー変数です。変量効果と観測誤差の事前分布は次のとおりです。
lme = fitlme(tbl,'y ~ InitialWeight + Program*Week + (Week|Subject)')
lme = Linear mixed-effects model fit by ML Model information: Number of observations 120 Fixed effects coefficients 9 Random effects coefficients 40 Covariance parameters 4 Formula: y ~ 1 + InitialWeight + Program*Week + (1 + Week | Subject) Model fit statistics: AIC BIC LogLikelihood Deviance -22.981 13.257 24.49 -48.981 Fixed effects coefficients (95% CIs): Name Estimate SE tStat DF {'(Intercept)' } 0.66105 0.25892 2.5531 111 {'InitialWeight' } 0.0031879 0.0013814 2.3078 111 {'Program_B' } 0.36079 0.13139 2.746 111 {'Program_C' } -0.033263 0.13117 -0.25358 111 {'Program_D' } 0.11317 0.13132 0.86175 111 {'Week' } 0.1732 0.067454 2.5677 111 {'Program_B:Week'} 0.038771 0.095394 0.40644 111 {'Program_C:Week'} 0.030543 0.095394 0.32018 111 {'Program_D:Week'} 0.033114 0.095394 0.34713 111 pValue Lower Upper 0.012034 0.14798 1.1741 0.022863 0.00045067 0.0059252 0.0070394 0.10044 0.62113 0.80029 -0.29319 0.22666 0.39068 -0.14706 0.3734 0.011567 0.039536 0.30686 0.68521 -0.15026 0.2278 0.74944 -0.15849 0.21957 0.72915 -0.15592 0.22214 Random effects covariance parameters (95% CIs): Group: Subject (20 Levels) Name1 Name2 Type Estimate {'(Intercept)'} {'(Intercept)'} {'std' } 0.18407 {'Week' } {'(Intercept)'} {'corr'} 0.66841 {'Week' } {'Week' } {'std' } 0.15033 Lower Upper 0.12281 0.27587 0.21076 0.88573 0.11004 0.20537 Group: Error Name Estimate Lower Upper {'Res Std'} 0.10261 0.087882 0.11981
値 0.022863 および 0.011567 は、被験者の初期体重と時間が体重の減少量に対して有意な影響を与えることを示しています。プログラム B の被験者の体重の減少は、プログラム A の被験者の体重の減少に対して有意差があります。変量効果の共分散パラメーターの下限または上限に 0 は含まれないため、この差は有意です。compare
メソッドを使用して、変量効果の有意性をテストすることもできます。
入力引数
formula
— モデル仕様の式
'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)'
という形式の文字ベクトルまたは string スカラー
モデル仕様の式。'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)'
という形式の文字ベクトルまたは string スカラーを指定します。この式では大文字小文字が区別されます。詳細は、式を参照してください。
例: 'y ~ treatment + (1|block)'
名前と値の引数
オプションの引数のペアを Name1=Value1,...,NameN=ValueN
として指定します。ここで Name
は引数名、Value
は対応する値です。名前と値の引数は他の引数の後ろにする必要がありますが、ペアの順序は関係ありません。
R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name
を引用符で囲みます。
例: 'CovariancePattern','Diagonal','Optimizer','fminunc','OptimizerOptions',opt
はモデルを指定します。ここで、変量効果項は対角共分散行列の構造体をもち、fitlme
は、変数 opt
で定義されたカスタム最適化パラメーターを使った fminunc
最適化アルゴリズムを使用します。
CovariancePattern
— 共分散行列のパターン
'FullCholesky'
(既定値) | 文字ベクトル | string スカラー | 対称正方 logical 行列 | string 配列 | 文字ベクトルまたは logical 行列の cell 配列
変量効果の共分散行列のパターン。'CovariancePattern'
と文字ベクトル、string スカラー、対称正方 logical 行列、string 配列、文字ベクトルの cell 配列、または logical 行列の cell 配列から構成されるコンマ区切りのペアとして指定します。
変量効果の項が R 個ある場合、'CovariancePattern'
の値は長さ R の string 配列または cell 配列でなければなりません。配列の各要素 r では、r 番目の変量効果の項に関連付けられている変量効果ベクトルの共分散行列のパターンを指定します。各要素のオプションは以下のとおりです。
'FullCholesky' | 既定の設定。コレスキー パラメーター表現を使用したフルの共分散行列。fitlme は、共分散行列のすべての要素を推定します。 |
'Full' | 対数コレスキー パラメーター表現を使用したフルの共分散行列。fitlme は、共分散行列のすべての要素を推定します。 |
'Diagonal' | 対角共分散行列。つまり、共分散行列の非対角要素は 0 に制約されます。 |
'Isotropic' | 分散が等しい対角共分散行列。つまり、共分散行列の非対角要素は 0 に制約され、対角要素は等価に制約されます。たとえば、等方性共分散構造をもつ変量効果の項が 3 つある場合、この共分散行列は次のようになります。 σ2b は、変量効果項の共通分散です。 |
'CompSymm' | 複合対称構造。つまり、対角線上の共通分散とすべての変量効果間の等しい相関です。たとえば、複合対称構造の共分散行列をもつ変量効果の項が 3 つある場合、この共分散行列は次のようになります。 σ2b1 は変量効果項の共通分散、σb1,b2 は任意の 2 つの変量効果項の間の共通共分散です。 |
PAT | 対称正方 logical 行列。'CovariancePattern' が行列 PAT によって定義されており、PAT(a,b) = false の場合、対応する共分散行列の要素 (a,b) は 0 に制約されます。 |
例: 'CovariancePattern','Diagonal'
例: 'CovariancePattern',{'Full','Diagonal'}
データ型: char
| string
| logical
| cell
FitMethod
— パラメーターを推定するメソッド
'ML'
(既定値) | 'REML'
線形混合効果モデルのパラメーターを推定するメソッド。'FitMethod'
と以下のいずれかで構成されるコンマ区切りのペアとして指定します。
'ML' | 既定の設定。最尤推定法 |
'REML' | 制限付き最尤推定法 |
例: 'FitMethod','REML'
Weights
— 観測値の重み
スカラー値のベクトル
観測値の重み。'Weights'
と、長さ n のベクトルで構成されるコンマ区切りのペアとして指定します。ここで、n は観測値の数です。
データ型: single
| double
Exclude
— 除外する行のインデックス
NaNs
を含まないすべての行を使用 (既定値) | 整数または論理値のベクトル
データ内の線形混合効果モデルから除外する行のインデックス。'Exclude'
と、整数または論理値のベクトルで構成されるコンマ区切りのペアとして指定します。
たとえば、以下のようにして、近似から 13 番目と 67 番目の行を除外できます。
例: 'Exclude',[13,67]
データ型: single
| double
| logical
DummyVarCoding
— ダミー変数に対して使用するコーディング
'reference'
(既定値) | 'effects'
| 'full'
カテゴリカル変数から作成されたダミー変数に対して使用するコーディング。'DummyVarCoding'
と次の表に記載された変数のいずれかで構成されるコンマ区切りのペアとして指定します。
値 | 説明 |
---|---|
'reference' (既定の設定) | fitlme は、基準グループを使用してダミー変数を作成します。この方式では、最初のカテゴリを基準グループとして扱い、カテゴリの数よりも 1 つ少ないダミー変数を作成します。カテゴリカル変数のカテゴリの順序は、関数 categories を使用してチェックできます。順序を変更するには、関数 reordercats を使用します。 |
'effects' | fitlme は、エフェクト コーディングを使用してダミー変数を作成します。この方式では、–1 を使用して最後のカテゴリを表します。この方式では、カテゴリの数よりも 1 つ少ないダミー変数を作成します。 |
'full' | fitlme は、完全なダミー変数を作成します。この方式では、各カテゴリに対して 1 つのダミー変数を作成します。 |
ダミー変数の作成に関する詳細については、ダミー変数の自動作成を参照してください。
例: 'DummyVarCoding','effects'
Optimizer
— 最適化アルゴリズム
'quasinewton'
(既定値) | 'fminunc'
最適化アルゴリズム。'Optimizer'
と以下のいずれかで構成されるコンマ区切りのペアとして指定します。
'quasinewton' | 既定の設定。信頼領域ベースの準ニュートン オプティマイザ―を使用します。statset('LinearMixedModel') を使用して、アルゴリズムのオプションを変更します。オプションを指定しない場合、LinearMixedModel は、statset('LinearMixedModel') の既定のオプションを使用します。 |
'fminunc' | このオプションを指定するには、Optimization Toolbox™ がなければなりません。optimoptions('fminunc') を使用して、アルゴリズムのオプションを変更します。オプションを指定しない場合、LinearMixedModel は、'Algorithm' が 'quasi-newton' に設定された optimoptions('fminunc') の既定のオプションを使用します。 |
例: 'Optimizer','fminunc'
OptimizerOptions
— 最適化アルゴリズムのオプション
statset
によって返される構造体 | optimoptions
によって返されるオブジェクト
最適化アルゴリズムのオプション。'OptimizerOptions'
と、statset('LinearMixedModel')
で返される構造体または optimoptions('fminunc')
で返されるオブジェクトで構成されるコンマ区切りのペアとして指定します。
'Optimizer'
が'fminunc'
の場合、optimoptions('fminunc')
を使用して、最適化アルゴリズムのオプションを変更します。'fminunc'
が使用するオプションについては、optimoptions
を参照してください。'Optimizer'
が'fminunc'
で、'OptimizerOptions'
を指定しない場合、LinearMixedModel
の既定の設定は、'Algorithm'
が'quasi-newton'
に設定されたoptimoptions('fminunc')
によって作成された既定のオプションになります。'Optimizer'
が'quasinewton'
の場合は、statset('LinearMixedModel')
を使用して、最適化パラメーターを変更します。最適化パラメーターを変更しない場合、LinearMixedModel
は、statset('LinearMixedModel')
によって作成された既定のオプションを使用します。
'quasinewton'
オプティマイザ―は、statset('LinearMixedModel')
によって作成された構造体の以下のフィールドを使用します。
TolFun
— 目的関数の勾配の相対許容誤差
1e-6
(既定値) | 正のスカラー値
目的関数の勾配の相対許容誤差。正のスカラー値として指定します。
TolX
— ステップ サイズの絶対許容誤差
1e-12
(既定値) | 正のスカラー値
ステップ サイズの絶対許容誤差。正のスカラー値として指定します。
MaxIter
— 許容される最大反復回数
10000
(既定値) | 正のスカラー値
許容される最大反復回数。正のスカラー値として指定します。
Display
— 表示のレベル
'off'
(既定値) | 'iter'
| 'final'
表示のレベル。'off'
、'iter'
、'final'
のいずれかとして指定します。
StartMethod
— 反復最適化を開始するメソッド
'default'
(既定値) | 'random'
反復最適化を開始するメソッド。'StartMethod'
と以下のいずれかで構成されるコンマ区切りのペアとして指定します。
値 | 説明 |
---|---|
'default' | 内部で定義される既定値 |
'random' | ランダムな初期値 |
例: 'StartMethod','random'
Verbose
— 画面に最適化プロセスを表示するインジケーター
false
(既定値) | true
画面に最適化プロセスを表示するインジケーター。'Verbose'
と、false
または true
で構成されるコンマ区切りのペアとして指定します。既定値は false
です。
'Verbose'
の設定は、'OptimizerOptions'
の 'Display'
フィールドをオーバーライドします。
例: 'Verbose',true
CheckHessian
— ヘッシアンの正定性をチェックするインジケーター
false
(既定値) | true
収束における制約のないパラメーターに関して、目的関数のヘッシアンの正定性をチェックするインジケーター。'CheckHessian'
と、false
または true
で構成されるコンマ区切りのペアとして指定します。既定値は false
です。
解の最適性を確認する場合や、共分散パラメーターの個数に関してモデルのパラメーターが多すぎるかどうかを判定する場合、'CheckHessian'
として true
を指定します。
例: 'CheckHessian',true
出力引数
lme
— 線形混合効果モデル
LinearMixedModel
オブジェクト
線形混合効果モデル。LinearMixedModel
オブジェクトとして返されます。
詳細
式
一般に、モデル仕様の式は 'y ~ terms'
という形式の文字ベクトルまたは string スカラーです。線形混合効果モデルでは、この式は 'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)'
の形式になります。ここで、fixed
および random
には固定効果および変量効果の項が含まれます。
テーブル tbl
に以下のものが格納されていると仮定します。
応答変数
y
連続変数またはグループ化変数である予測子変数
Xj
グループ化変数
g1
、g2
、...、gR
ここで、Xj
および gr
のグループ化変数は、categorical 配列、logical 配列、文字配列、string 配列、または文字ベクトルの cell 配列が可能です。
この場合、'y ~ fixed + (random1|g1) + ... + (randomR|gR)'
の形式の式において、項 fixed
は固定効果の計画行列 X
の仕様に対応し、random
1 はグループ化変数 g
1 に対応する変量効果の計画行列 Z
1 の仕様であり、同様に random
R はグループ化変数 g
R に対応する変量効果の計画行列 Z
R の仕様です。fixed
項および random
項はウィルキンソンの表記法で表現できます。
ウィルキンソンの表記法は、モデルに存在する因子を記述します。この表記法は、モデルに存在する因子に関係するものであり、それらの因子の乗数 (係数) に関係するものではありません。
ウィルキンソンの表記法 | 標準表記の因子 |
---|---|
1 | 定数 (切片) 項 |
X^k 、k は正の整数 | X , X2 , ..., Xk |
X1 + X2 | X1 , X2 |
X1*X2 | X1 , X2 , X1.*X2 (elementwise multiplication of X1 and X2) |
X1:X2 | X1.*X2 のみ |
- X2 | X2 は含めない |
X1*X2 + X3 | X1 , X2 , X3 , X1*X2 |
X1 + X2 + X3 + X1:X2 | X1 , X2 , X3 , X1*X2 |
X1*X2*X3 - X1:X2:X3 | X1 , X2 , X3 , X1*X2 , X1*X3 , X2*X3 |
X1*(X2 + X3) | X1 , X2 , X3 , X1*X2 , X1*X3 |
Statistics and Machine Learning Toolbox™ 表記は、-1
を使用して項を明示的に削除しない限り、常に定数項を含みます。線形混合効果モデルの仕様例を次にいくつか挙げます。
次に例を示します。
式 | 説明 |
---|---|
'y ~ X1 + X2' | 切片 X1 および X2 の固定効果。これは、'y ~ 1 + X1 + X2' と等価です。 |
'y ~ -1 + X1 + X2' | X1 と X2 の切片と固定効果はありません。-1 を含めることによって暗黙的な切片の項は抑制されます。 |
'y ~ 1 + (1 | g1)' | グループ化変数 g1 のレベルごとの切片の固定効果と切片の変量効果の和。 |
'y ~ X1 + (1 | g1)' | 固定勾配のランダム切片モデル。 |
'y ~ X1 + (X1 | g1)' | 相関があり得るランダムな切片と勾配。これは、'y ~ 1 + X1 + (1 + X1|g1)' と等価です。 |
'y ~ X1 + (1 | g1) + (-1 + X1 | g1)' | 切片と勾配の独立した変量効果項。 |
'y ~ 1 + (1 | g1) + (1 | g2) + (1 | g1:g2)' | g1 と g2 に対する独立したメイン効果のあるランダムな切片モデル + 独立した交互作用効果。 |
コレスキー パラメーター化
線形混合効果モデルの仮定の 1 つに、変量効果の事前分布が次のようになっているというものがあります。
ここで、D は分散成分ベクトル θ によってパラメーター表現された q 行 q 列の対称な半正定値行列、q は変量効果項の変数の数、σ2 は観測誤差の分散です。変量効果 D の共分散行列は対称なので、q(q+1)/2 の自由パラメーターをもちます。L が次の条件を満たす D(θ) の下三角コレスキー因子であるとします。
すると、制約がない q*(q+1)/2 行 1 列のパラメーター ベクトル θ が L の下三角部分にある要素から形成されます。
たとえば、次のようになっているとします。
この場合、次のようになります。
対数コレスキー パラメーター表現
コレスキー パラメーター表現の L の対角要素が正に制約されている場合、L の解は一意になります。対数コレスキー パラメーター表現は、一意のパラメーター表現を保証するために L の対角要素の対数を使用するという点を除いて、コレスキー パラメーター表現と同じです。
たとえば、3 行 3 列のコレスキー パラメーター表現の例では、Lii ≥ 0 を適用すると次のようになります。
代替方法
式で表現することが困難なモデルの場合は、固定効果と変量効果を定義する行列を作成し、fitlmematrix(X,y,Z,G)
を使用してモデルを当てはめることができます。
参照
[1] Pinherio, J. C., and D. M. Bates. “Unconstrained Parametrizations for Variance-Covariance Matrices”. Statistics and Computing, Vol. 6, 1996, pp. 289–296.
バージョン履歴
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)