ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

fitlmematrix

線形混合効果モデルの近似

構文

  • lme = fitlmematrix(X,y,Z,[])
  • lme = fitlmematrix(X,y,Z,G)
  • lme = fitlmematrix(___,Name,Value)

説明

lme = fitlmematrix(X,y,Z,[]) は、固定効果の計画行列 XZ 内の 1 つ以上の変量効果の計画行列を使用して、応答 y の線形混合効果モデルを作成します。

[] は、1 つのグループがあることを意味しています。つまり、グループ化変数 Gones(n,1) になります。ここで、n は観測値の数です。共分散パターンを指定せずに fitlmematrix(X,Y,Z,[]) を使用すると、認識されないモデルになる可能性が高くなります。グループ化情報を変量効果計画 Z にビルドし、名前と値のペアの引数 'CovariancePattern' を使用して変量効果の共分散パターンを指定した場合にのみ、この構文は推奨されます。

lme = fitlmematrix(X,y,Z,G) は、固定効果の計画行列 X、変量効果の計画行列 Z または Z 内の複数の変量効果の計画行列、あるいは G 内の 1 つまたは複数のグループ化変数を使用して、応答 y の線形混合効果モデルを作成します。

また、lme = fitlmematrix(___,Name,Value) は、1 つまたは複数の Name,Value の引数のペアによって指定された追加オプションを使用して線形混合効果モデルを作成します。

たとえば、応答変数、予測子変数、グループ化変数の名前を指定できます。共分散パターン、近似法または最適化アルゴリズムを指定することもできます。

すべて折りたたむ

グループ化変数は指定されない

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

load carsmall

線形混合効果モデルを近似します。ここで、ガロンあたりの走行マイル数 (MPG) が応答で、重みが予測子変数であり、切片はモデル年度に基づいて変化します。最初に、計画行列を定義します。次に、指定した計画行列を使用してモデルを近似します。

y = MPG;
X = [ones(size(Weight)), Weight];
Z = ones(size(y));
lme = fitlmematrix(X,y,Z,Model_Year)
lme = 


Linear mixed-effects model fit by ML

Model information:
    Number of observations              94
    Fixed effects coefficients           2
    Random effects coefficients          3
    Covariance parameters                2

Formula:
    y ~ x1 + x2 + (z11 | g1)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    486.09    496.26    -239.04          478.09  

Fixed effects coefficients (95% CIs):
    Name        Estimate      SE           tStat      DF    pValue        Lower         Upper     
    'x1'            43.575       2.3038     18.915    92    1.8371e-33            39        48.151
    'x2'        -0.0067097    0.0004242    -15.817    92    5.5373e-28    -0.0075522    -0.0058672

Random effects covariance parameters (95% CIs):
Group: g1 (3 Levels)
    Name1        Name2        Type         Estimate    Lower     Upper 
    'z11'        'z11'        'std'        3.301       1.4448    7.5421

Group: Error
    Name             Estimate    Lower     Upper 
    'Res Std'        2.8997      2.5075    3.3532

ここで、グループを行列 Z にビルドすることにより、同じモデルを近似します。

Z = double([Model_Year==70, Model_Year==76, Model_Year==82]);
lme = fitlmematrix(X,y,Z,[],'Covariancepattern','Isotropic')
lme = 


Linear mixed-effects model fit by ML

Model information:
    Number of observations              94
    Fixed effects coefficients           2
    Random effects coefficients          3
    Covariance parameters                2

Formula:
    y ~ x1 + x2 + (z11 + z12 + z13 | g1)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    486.09    496.26    -239.04          478.09  

Fixed effects coefficients (95% CIs):
    Name        Estimate      SE           tStat      DF    pValue        Lower         Upper     
    'x1'            43.575       2.3038     18.915    92    1.8371e-33            39        48.151
    'x2'        -0.0067097    0.0004242    -15.817    92    5.5373e-28    -0.0075522    -0.0058672

Random effects covariance parameters (95% CIs):
Group: g1 (1 Levels)
    Name1        Name2        Type         Estimate    Lower     Upper 
    'z11'        'z11'        'std'        3.301       1.4448    7.5421

Group: Error
    Name             Estimate    Lower     Upper 
    'Res Std'        2.8997      2.5075    3.3532

共分散を使用した長期間の調査

標本データが含まれたフォルダーに移動します。

cd(matlabroot)
cd('help/toolbox/stats/examples')

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

load weight

weight には長期間の調査によるデータが含まれています。そこには 20 人の被験者が 4 つの運動プログラム (A、B、C、D) にランダムに割り当てられ、体重の減少が 6 回の 2 週間の期間にわたって記録されています。このデータは、シミュレーションされたものです。

Subject および Program をカテゴリカル変数として定義します。線形混合効果モデルの計画行列を作成します。ここで、初期体重、プログラムの種類、週、週とプログラムの種類の間の交互作用は固定効果です。切片と週の係数は被験者ごとに異なります。

このモデルは以下の式に対応します。

yim=β0+β1IWi+β2Weeki+β3I[PB]i+β4I[PC]i+β5I[PD]i+β6(Weeki*I[PB]i)+β7(Weeki*I[PC]i)+β8(Weeki*I[PD]i)+b+0mb1mWeekim+εim,

ここで、 i = 1, 2, ..., 120 および m = 1, 2, ..., 20. βj は固定効果係数で、j = 0, 1, ..., 8、b0m および b1m は変量効果です。IW は初期体重を表し、I[.] はプログラムのタイプを示すダミー変数です。たとえば、I[PB]i はプログラム タイプ B を表すダミー変数です。変量効果と観測値の誤差の事前分布は次のとおりです。b0m~N(0,σ20)、b1m~N(0,σ21)、εim ~ N(0,σ2) です。

Subject = nominal(Subject);
Program = nominal(Program);
D = dummyvar(Program); % Create dummy variables for Program
X = [ones(120,1), InitialWeight, D(:,2:4), Week,...
		 D(:,2).*Week, D(:,3).*Week, D(:,4).*Week];
Z = [ones(120,1), Week];
G = Subject;

モデルには切片があるので、プログラム B、C、D のダミー変数のみが必要です。これは、ダミー変数をコーディングする 'reference' メソッドとも呼ばれます。

定義した計画行列とグループ化変数で fitlmematrix を使用してモデルを近似します。

lme = fitlmematrix(X,y,Z,G,'FixedEffectPredictors',...
{'Intercept','InitWeight','PrgB','PrgC','PrgD','Week','Week_PrgB','Week_PrgC','Week_PrgD'},...
'RandomEffectPredictors',{{'Intercept','Week'}},'RandomEffectGroups',{'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:
    Linear Mixed Formula with 10 predictors.

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     pValue       Lower         Upper    
    'Intercept'           0.66105      0.25892      2.5531    111     0.012034       0.14798       1.1741
    'InitWeight'        0.0031879    0.0013814      2.3078    111     0.022863    0.00045067    0.0059252
    'PrgB'                0.36079      0.13139       2.746    111    0.0070394       0.10044      0.62113
    'PrgC'              -0.033263      0.13117    -0.25358    111      0.80029      -0.29319      0.22666
    'PrgD'                0.11317      0.13132     0.86175    111      0.39068      -0.14706       0.3734
    'Week'                 0.1732     0.067454      2.5677    111     0.011567      0.039536      0.30686
    'Week_PrgB'          0.038771     0.095394     0.40644    111      0.68521      -0.15026       0.2278
    'Week_PrgC'          0.030543     0.095394     0.32018    111      0.74944      -0.15849      0.21957
    'Week_PrgD'          0.033114     0.095394     0.34713    111      0.72915      -0.15592      0.22214

Random effects covariance parameters (95% CIs):
Group: Subject (20 Levels)
    Name1              Name2              Type          Estimate    Lower      Upper  
    'Intercept'        'Intercept'        'std'         0.18407     0.12281    0.27587
    'Week'             'Intercept'        'corr'        0.66841     0.21076    0.88573
    'Week'             'Week'             'std'         0.15033     0.11004    0.20537

Group: Error
    Name             Estimate    Lower       Upper  
    'Res Std'        0.10261     0.087882    0.11981

p 値 0.0228 および 0.0115 は、減少した体重に関して被験者の初期の体重と時間因子は有意な影響があることを示しています。プログラム B の被験者の体重の減少は、プログラム A の被験者の体重の減少に対して有意差があります。変量効果の共分散パラメーターの下限または上限に 0 は含まれないため、この差は有意であるように見えます。compare メソッドを使用して、変量効果の有意性をテストすることもできます。

ランダムな切片モデル

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

load flu

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

インフルエンザ罹患率を応答として線形混合効果モデルを近似するには、地域に対応する 9 個の列を、単一の応答変数 FluRate、ノミナル変数 Region、各推定の元になっている地域を示す全国の推定値 WtdILI、グループ化変数 Date をもつ 1 つの高さ配列にまとめます。

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

ランダムな切片の線形混合効果モデルの計画行列を定義します。ここで、切片は Date で変化します。対応するモデルは以下のとおりです。

yim=β0+β1WtdILIim+b0m+εim,i=1,2,...,468,m=1,2,...,52,

ここで、yim は、グループ化変数 Date のレベル m の観測値 i です。b0m はグループ化変数 Date のレベル m の変量効果で、εim は観測値 i の観測誤差です。変量効果は事前分布 b0m ~ N(0,σ2FR) をもち、誤差項は分布 εim ~ N(0,σ2) をもちます。

y = flu2.FluRate;
X = [ones(468,1) flu2.WtdILI];
Z = [ones(468,1)];
G = flu2.Date;

線形混合効果モデルを近似します。

lme = fitlmematrix(X,y,Z,G,'FixedEffectPredictors',{'Intercept','NationalRate'},...
'RandomEffectPredictors',{{'Intercept'}},'RandomEffectGroups',{'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:
    y ~ Intercept + NationalRate + (Intercept | 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        Lower       Upper  
    'Intercept'           0.16385     0.057525    2.8484    466     0.0045885    0.050813    0.27689
    'NationalRate'         0.7236     0.032219    22.459    466    3.0502e-76     0.66028    0.78691

Random effects covariance parameters (95% CIs):
Group: Date (52 Levels)
    Name1              Name2              Type         Estimate    Lower      Upper  
    'Intercept'        'Intercept'        'std'        0.17146     0.13227    0.22226

Group: Error
    Name             Estimate    Lower      Upper  
    'Res Std'        0.30201     0.28217    0.32324

変量効果の項 σ2b0m の標準偏差の信頼限界には、変量効果の項が有意であることを示す 0 (0.13227, 0.22226) が含まれません。compare メソッドを使用して、変量効果の有意性をテストすることもできます。

観測値の推定値は、その観測値に対応するグループ化変数レベルの固定効果値および変量効果値の合計です。たとえば、観測値 28 のインフルエンザ推定羅患率は以下のとおりです。

y^28=β^0+β^1WtdILI28+b^10/30/2005=0.1639+0.7236*(1.343)+0.3318=1.46749,

ここで、b^ は切片の変量効果の最良線形不偏予測子 (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(lme) メソッドを使用して、近似値を簡単に表示できます。

F = fitted(lme);
F(28)
ans =

    1.4674

乱塊法

標本データが含まれたフォルダーに移動します。

cd(matlabroot)
cd('help/toolbox/stats/examples')

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

load shift

このデータは 5 人の作業者が 3 つのシフトの間に製造した製品から計測された品質目標の特性の絶対偏差を示します。3 つのシフトとは朝、夕方、夜です。これは作業者をブロックとする乱塊法です。この実験は、シフトの時間によるパフォーマンスへの影響の調査を意図しています。パフォーマンスの測定は、目標値からの品質特性の偏差です。このデータは、シミュレーションされたものです。

作業者でグループ化されたランダムな切片と固定効果としてのシフトをもつ線形混合効果モデルの計画行列を定義します。'effects' 対比を使用します。'effects' 対比は係数の合計が 0 であることを示します。固定効果の計画行列 X1 および X2 内に 2 つの対比コード変数を作成する必要があります。ここで以下のようになります。

Shift_Evening={0,ifMorning1,ifEvening1,ifNightandShift_Morning={1,ifMorning0,ifEvening1,ifNight.

このモデルは以下の式に対応します。

Morning Shift: QCDevim=β0+β2Shift_Morningi+b0m+εim,m=1,2,...,5,Evening Shift: QCDevim=β0+β1Shift_Eveningi+b0m+εim,Night Shift:     QCDevim=β0β1Shift_Eveningiβ2Shift_Morningi+b0m+εim,

ここで、i は観測値を表し、m は作業者を表します。i = 1、2、...、15 で、m = 1、2、...、5 です。変量効果と観測値の誤差の分布は b0m ~ N(0,σ2b0m) および εim ~ N(0,σ2) です。

S = shift.Shift;
X1 = (S=='Morning') - (S=='Night');
X2 = (S=='Evening') - (S=='Night');
X = [ones(15,1), X1, X2];
y = shift.QCDev;
Z = ones(15,1);
G = shift.Operator;

指定した計画行列と制限付き最尤法を使用して、線形混合効果モデルを近似します。

lme = fitlmematrix(X,y,Z,G,'FitMethod','REML','FixedEffectPredictors',....
{'Intercept','S_Morning','S_Evening'},'RandomEffectPredictors',{{'Intercept'}},...
'RandomEffectGroups',{'Operator'},'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:
    y ~ Intercept + S_Morning + S_Evening + (Intercept | 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       Lower      Upper   
    'Intercept'          3.6525    0.94109     3.8812    12    0.0021832     1.6021       5.703
    'S_Morning'        -0.91973    0.31206    -2.9473    12     0.012206    -1.5997    -0.23981
    'S_Evening'        -0.53293    0.31206    -1.7078    12      0.11339    -1.2129     0.14699

Random effects covariance parameters (95% CIs):
Group: Operator (5 Levels)
    Name1              Name2              Type         Estimate    Lower      Upper 
    'Intercept'        'Intercept'        'std'        2.0457      0.98207    4.2612

Group: Error
    Name             Estimate    Lower      Upper
    'Res Std'        0.85462     0.52357    1.395

変量効果の最良線形不偏予測子 (BLUP) の推定値を計算します。

B = randomEffects(lme)
B =

    0.5775
    1.1757
   -2.1715
    2.3655
   -1.9472

夕方シフトで作業している 3 番目の作業者の場合、品質目標の特性の推定偏差は以下のとおりです。

y^Evening,Operator3=β^0+β^1Shift_Evening+b^03=3.65250.532932.1715=0.94807.

この値を以下のように表示することもできます。

F = fitted(lme);
F(shift.Shift=='Evening' & shift.Operator=='3')
ans =

    0.9481

相関および無相関の変量効果の項

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

load carbig

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

MPGim=β0+β1Acci+β2HP+b0m+b1mAccim+εim,m=1,2,3,

また、変量効果項の分布は b0m ~ N(0,σ20)、b1m ~ N(0,σ21) です。m はモデル年を表します。

最初に、線形混合効果モデルを近似するための計画行列を準備します。

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

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

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


Linear mixed-effects model fit by ML

Model information:
    Number of observations             392
    Fixed effects coefficients           3
    Random effects coefficients         26
    Covariance parameters                3

Formula:
    y ~ Intercept + Acceleration + Horsepower + (Intercept | Model_Year) + (Acceleration | Model_Year)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    2194.5    2218.3    -1091.3          2182.5  

Fixed effects coefficients (95% CIs):
    Name                  Estimate    SE           tStat      DF     pValue        Lower       Upper   
    'Intercept'             49.839       2.0518     24.291    389    5.6168e-80      45.806      53.873
    'Acceleration'        -0.58565      0.10846    -5.3995    389    1.1652e-07     -0.7989     -0.3724
    'Horsepower'          -0.16534    0.0071227    -23.213    389    1.9755e-75    -0.17934    -0.15133

Random effects covariance parameters (95% CIs):
Group: Model_Year (13 Levels)
    Name1              Name2              Type         Estimate      Lower    Upper
    'Intercept'        'Intercept'        'std'        8.0669e-07    NaN      NaN  

Group: Model_Year (13 Levels)
    Name1                 Name2                 Type         Estimate    Lower      Upper  
    'Acceleration'        'Acceleration'        'std'        0.18783     0.12523    0.28172

Group: Error
    Name             Estimate    Lower     Upper 
    'Res Std'        3.7258      3.4698    4.0007

切片と加速度の変量効果共分散パラメーターが個別に表示されています。切片の変量効果の標準偏差は有意ではないようです。

切片と加速度について相関する可能性がある変量効果でモデルを再近似します。この場合、変量効果の項は以下の事前分布をもちます。

bm=(b0mb1m)~N(0,(σ02σ0,1σ0,1σ12)),

ここで、m はモデル年度を表します。

最初に、変量効果の計画行列とグループ化変数を準備します。

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

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


Linear mixed-effects model fit by ML

Model information:
    Number of observations             392
    Fixed effects coefficients           3
    Random effects coefficients         26
    Covariance parameters                4

Formula:
    y ~ Intercept + Acceleration + Horsepower + (Intercept + Acceleration | Model_Year)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    2193.5    2221.3    -1089.7          2179.5  

Fixed effects coefficients (95% CIs):
    Name                  Estimate    SE           tStat      DF     pValue        Lower       Upper   
    'Intercept'             50.133       2.2652     22.132    389    7.7727e-71      45.679      54.586
    'Acceleration'        -0.58327      0.13394    -4.3545    389    1.7075e-05    -0.84661    -0.31992
    'Horsepower'          -0.16954    0.0072609     -23.35    389     5.188e-76    -0.18382    -0.15527

Random effects covariance parameters (95% CIs):
Group: Model_Year (13 Levels)
    Name1                 Name2                 Type          Estimate    Lower       Upper   
    'Intercept'           'Intercept'           'std'           3.3475      1.2862      8.7119
    'Acceleration'        'Intercept'           'corr'        -0.87971    -0.98501    -0.29675
    'Acceleration'        'Acceleration'        'std'          0.33789      0.1825     0.62558

Group: Error
    Name             Estimate    Lower     Upper 
    'Res Std'        3.6874      3.4298    3.9644

切片と加速の変量効果共分散パラメーターがあわせて表示され、さらに切片と加速度の間の相関も含まれています。標準偏差の信頼区間と、切片に対する変量効果と加速度の間の相関には 0 は含まれていません。そのため、これらは有意と考えられます。compare メソッドを使用して、これらの 2 つのモデルを比較できます。

共分散パターンの指定

標本データが含まれたフォルダーに移動します。

cd(matlabroot)
cd('help/toolbox/stats/examples')

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

load weight

weight には長期間の調査によるデータが含まれています。そこには 20 人の被験者が 4 つの運動プログラムにランダムに割り当てられ、体重の減少が 6 回の 2 週間の期間にわたって記録されています。このデータは、シミュレーションされたものです。

Subject および Program をカテゴリカル変数として定義します。

Subject = nominal(Subject);
Program = nominal(Program);

初期体重、プログラムの種類および週を固定効果としてもつ線形混合効果モデルの計画行列を作成します。

D = dummyvar(Program);
X = [ones(120,1), InitialWeight, D(:,2:4), Week];
Z = [ones(120,1) Week];
G = Subject;

このモデルは以下の式に対応します。

yim=β0+β1IWi+β2Weeki+β3I[PB]i+β4I[PC]i+β5I[PD]i+b+0mb1mWeek2im+b2mWeek4im+b3mWeek6im+b4mWeek8im+b5mWeek10im+b6mWeek12im+εim,

ここで、i = 1, 2, ..., 120 および m = 1, 2, ..., 20 です。

βj は固定効果係数で、j = 0, 1, ..., 8、b1m および b1m は変量効果です。IW は初期体重を表し、I[.] はプログラムの種類を表すダミー変数です。たとえば、I[PB]i はプログラム タイプ B を表すダミー変数です。変量効果と観測値の誤差の事前分布は次のとおりです。b0m~N(0,σ20)、b1m~N(0,σ21)、εim ~ N(0,σ2) です。

定義した計画行列とグループ化変数で fitlmematrix を使用してモデルを近似します。1 人の被験者に関して繰り返し収集された観測値が対角線上の共通分散をもつと仮定します。

lme = fitlmematrix(X,y,Z,G,'FixedEffectPredictors',...
{'Intercept','InitWeight','PrgB','PrgC','PrgD','Week'},...
'RandomEffectPredictors',{{'Intercept','Week'}},...
'RandomEffectGroups',{'Subject'},'CovariancePattern','Isotropic')
lme = 


Linear mixed-effects model fit by ML

Model information:
    Number of observations             120
    Fixed effects coefficients           6
    Random effects coefficients         40
    Covariance parameters                2

Formula:
    y ~ Intercept + InitWeight + PrgB + PrgC + PrgD + Week + (Intercept + Week | Subject)

Model fit statistics:
    AIC        BIC       LogLikelihood    Deviance
    -24.783    -2.483    20.391           -40.783 

Fixed effects coefficients (95% CIs):
    Name                Estimate     SE           tStat       DF     pValue        Lower        Upper    
    'Intercept'            0.4208      0.28169      1.4938    114       0.13799     -0.13723      0.97883
    'InitWeight'        0.0045552    0.0015338      2.9699    114     0.0036324    0.0015168    0.0075935
    'PrgB'                0.36993      0.12119      3.0525    114     0.0028242      0.12986         0.61
    'PrgC'              -0.034009       0.1209    -0.28129    114       0.77899     -0.27351       0.2055
    'PrgD'                  0.121      0.12111     0.99911    114       0.31986     -0.11891      0.36091
    'Week'                0.19881     0.037134      5.3538    114    4.5191e-07      0.12525      0.27237

Random effects covariance parameters (95% CIs):
Group: Subject (20 Levels)
    Name1              Name2              Type         Estimate    Lower      Upper  
    'Intercept'        'Intercept'        'std'        0.16561     0.12896    0.21269

Group: Error
    Name             Estimate    Lower       Upper  
    'Res Std'        0.10272     0.088014    0.11987 

入力引数

すべて折りたたむ

X — 固定効果の計画行列n 行 p 列の行列

n 行 p 列の行列として指定される固定効果の計画行列。ここで、n は観測値の数、p は固定効果の予測子変数の数です。X の各行は 1 つの観測値に対応し、X の各列は 1 つの変数に対応します。

データ型: single | double

y — 応答値n 行 1 列のベクトル

n 行 1 列のベクトルとして指定された応答値。ここで、n は観測値の数です。

データ型: single | double

Z — 変量効果の計画n 行 q 列の行列 | R n 行 q(r) 列の行列のセル配列、r = 1, 2, ..., R

変量効果の計画。以下のいずれかとして指定します。

  • モデルに 1 個の変量効果項がある場合、Z は n 行 q 列の行列でなければなりません。n は観測値の数、q は変量効果項の変数の数です。

  • R 個の変量効果項がある場合、Z は長さ R のセル配列でなければなりません。Z の各セルには、各変量効果項に対応する、n 行 q (r) 列の計画行列 Z{r}、r = 1, 2, ..., R が含まれています。ここで、q (r) は、r 番目の変量効果の計画行列 Z{r} の変量効果項の数です。

データ型: single | double | cell

G — グループ化変数n 行 1 列のベクトル | R n 行 1 列のベクトルのセル配列

グループ化変数。次のいずれかとして指定します。

  • 1 個の変量効果項がある場合、G は、M 個のレベルまたはグループをもつ単一のグループ化変数に対応する n 行 1 列のベクトルでなければなりません。

    G には、カテゴリカル ベクトル、数値ベクトル、文字配列または文字列のセル配列のいずれかを指定できます。

  • 複数の変量効果項がある場合、G は長さ R のセル配列でなければなりません。G の各セルには、M(r) 個のレベルをもつグループ化変数 G{r}、r = 1, 2, ..., R が含まれています。

    G{r} には、カテゴリカル ベクトル、数値ベクトル、文字配列または文字列のセル配列のいずれかを指定できます。

データ型: single | double | char | cell

名前/値のペアの引数

オプションの Name,Value の引数ペアをコンマ区切りで指定します。ここで、Name は引数名で、Value は対応する値です。Name は単一引用符 (' ') で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を任意の順番で指定できます。

例: 'CovariancePattern','Diagonal','DummyVarCoding','full','Optimizer','fminunc' は、0 の非対角要素をもつ変量効果の共分散パターンを指定し、カテゴリカル変数の各レベルに対してダミー変数を作成し、最適化アルゴリズム fminunc を使用します。

'FixedEffectPredictors' — 固定効果の計画行列の列名{'x1','x2',...,'xP'} (既定値) | 長さ p のセル配列

固定効果の計画行列 X の列名。'FixedEffectPredictors' と、長さ p のセル配列で構成されるコンマ区切りのペアとして指定します。

たとえば、固定効果として、1 つの定数項と 2 つの予測子、つまり TimeSpent および Gender がある場合 (ここで、FemaleGender の基準レベルです)、以下の方法で固定効果の名前を指定できます。Gender_Male は、カテゴリ Male に対して作成しなければならないダミー変数です。これらの変数に対して、別々の名前を選択できます。

例: 'FixedEffectPredictors',{'Intercept','TimeSpent','Gender_Male'},

データ型: cell

'RandomEffectPredictors' — 変量効果の計画行列の列名またはセル配列長さ q のセル配列 | 長さ q(r) の要素をもつ長さ R のセル配列、r = 1, 2, ..., R

変量効果の計画行列の列名またはセル配列 Z'RandomEffectPredictors' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

  • Z が n 行 q 列の計画行列であるときは、長さ q のセル配列。この場合、既定の設定は {'z1','z2',...,'zQ'} です。

  • Z が、長さ q(r)、r = 1, 2, ..., R の各要素 Z{r} をもつ長さ R のセル配列であるときは、長さ R のセル配列。この場合、既定の設定は {'z11','z12',...,'z1Q(1)'},...,{'zr1','zr2',...,'zrQ(r)'} です。

たとえば、切片と Acceleration という名前の変数に対する変量効果に相関があるとします。この場合、以下のように変量効果の予測子名を指定できます。

例: 'RandomEffectPredictors',{'Intercept','Acceleration'}

2 つの変量効果項、つまり変数 g1 によってグループ化された切片と変数 Acceleration に対する変量効果項と、変数 g2 によってグループ化された切片に対する変量効果項がある場合、変量効果の予測子名を以下のように指定できます。

例: 'RandomEffectPredictors',{{'Intercept','Acceleration'},{'Intercept'}}

データ型: cell

'ResponseVarName' — 応答変数の名前'y' (既定値) | 文字列

応答変数の名前。'ResponseVarName' と文字列で構成されるコンマ区切りのペアとして指定します。

たとえば、応答変数の名前が score の場合、以下のように指定できます。

例: 'ResponseVarName','score'

データ型: char

'RandomEffectGroups' — 変量効果のグループ化変数の名前'g' または {'g1','g2',...,'gR'} (既定値) | 文字列 | 文字列のセル配列

変量効果のグループ化変数の名前。'RandomEffectGroups' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

  • 文字列 — 変量効果項が 1 つしかない場合、つまり G がベクトルの場合、'RandomEffectGroups' の値はグループ化変数 G の名前を含む文字列になります。既定値は 'g' です。

  • 文字列のセル配列 — 複数の変量効果項がある場合、つまり G が長さ R のセル配列の場合、'RandomEffectGroups' の値は長さ R のセル配列になります。ここで、各セルにはグループ化変数 G{r} の名前が含まれています。既定値は {'g1','g2',...,'gR'} です。

たとえば、グループ化変数 sex および subject でグループ化された 2 つの変量効果項 z1 および z2 がある場合、グループ化変数の名前を以下のように指定できます。

例: 'RandomEffectGroups',{'sex','subject'}

データ型: char | cell

'CovariancePattern' — 共分散行列のパターン'FullCholesky' (既定値) | 文字列 | 対称正方論理行列 | 文字列または論理行列のセル配列

変量効果の共分散行列のパターン。'CovariancePattern' および文字列と、対称正方論理行列、あるいは文字列または論理行列のセル配列で構成されるコンマ区切りのペアとして指定します。

変量効果の項 R がある場合、'CovariancePattern' の値は長さ R のセル配列にしなければなりません。このセル配列の各要素 r は、r 番目の変量効果の項に関連付けられた変量効果ベクトルの共分散行列のパターンを指定します。各要素のオプションは以下のとおりです。

'FullCholesky'既定の設定。コレスキー パラメーター化を使用した完全な共分散行列。fitlme は、共分散行列のすべての要素を推定します。
'Full'対数コレスキー パラメーター化を使用した完全な共分散行列。fitlme は、共分散行列のすべての要素を推定します。
'Diagonal'対角共分散行列。つまり、共分散行列の非対角要素は 0 に制約されます。

(σb12000σb22000σb32)

'Isotropic'分散が等しい対角共分散行列。つまり、共分散行列の非対角要素は 0 に制約され、対角要素は等価に制約されます。たとえば、等方性共分散構造をもつ変量効果の項が 3 つある場合、この共分散行列は以下のようになります。

(σb2000σb2000σb2)

ここで、σ2b は、変量効果の項の共通分散です。

'CompSymm'複合対称構造。つまり、対角線上の共通分散とすべての変量効果間の等しい相関です。たとえば、複合対称構造をもつ共分散行列を使用した変量効果の項が 3 つある場合、この共分散行列は以下のようになります。

(σb12σb1,b2σb1,b2σb1,b2σb12σb1,b2σb1,b2σb1,b2σb12)

ここで、σ2b1 は変量効果の項の共通分散で、σb1,b2 は 2 つの任意の変量効果の項間の共通共分散です。
PAT対称正方論理行列。'CovariancePattern' が行列 PAT によって定義されていて、PAT(a,b) = false の場合は、対応する共分散行列の要素 (a,b) が 0 に制約されます。

例: 'CovariancePattern','Diagonal'

例: 'CovariancePattern',{'Full','Diagonal'}

'FitMethod' — パラメーターを推定するメソッド'ML' (既定値) | 'REML'

線形混合効果モデルのパラメーターを推定するメソッド。'FitMethod' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

'ML'既定の設定。最尤推定法
'REML'制限付き最尤推定法

例: 'FitMethod','REML'

'Weights' — 観測値の重みスカラー値のベクトル

観測値の重み。'Weights' と、長さ n のベクトルで構成されるコンマ区切りのペアとして指定します。ここで、n は観測値の数です。

データ型: single | double

'Exclude' — 除外する行のインデックスNaN を含まないすべての行を使用 (既定値) | 整数または論理値のベクトル

データ内の線形混合効果モデルから除外する行のインデックス。'Exclude' と、整数または論理値のベクトルで構成されるコンマ区切りのペアとして指定します。

たとえば、以下のようにして、近似から 13 番目と 67 番目の行を除外できます。

例: 'Exclude',[13,67]

データ型: single | double | logical

'DummyVarCoding' — ダミー変数に対して使用するコーディング'reference' (既定値) | 'effects' | 'full'

カテゴリカル変数から作成されたダミー変数に対して使用するコーディング。'DummyVarCoding' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

'reference'既定の設定。0 に設定された最初のカテゴリの係数。
'effects'合計 0 の係数。
'full'各カテゴリに対して 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 オブジェクトとして返される線形混合効果モデル。

このオブジェクトのプロパティとメソッドについては、LinearMixedModel を参照してください。

代替機能

fitlme(tbl,formula) を使用して、線形混合効果モデルを近似することもできます。ここで、tbl は応答 y を含むテーブルまたはデータセット配列であり、予測子変数 X、グループ化変数、formula'y ~ fixed + (random1|g1) + ... + (randomR|gR)' の形式です。

詳細

すべて折りたたむ

コレスキー パラメーター化

線形混合効果モデルの仮定の 1 つは、変量効果が以下の事前分布をもつということです。

b~N(0,σ2D(θ)),

ここで、D は q 行 q 列の、分散成分ベクトル θ でパラメーター化された対称半正定行列、q は変量効果項の変数の数、σ2 は観測誤差分散です。変量効果 D の共分散行列は対称なので、q(q+1)/2 の自由パラメーターをもちます。L が D(θ) の下三角コレスキー因子であり、次の条件をみたすものとします。

D(θ)=L(θ)L(θ)T,

この場合、q*(q+1)/2 行 1 列の制約のないパラメーター ベクトル θ は、L の下三角部分の要素から成ります。

たとえば、

L=[L1100L21L220L31L32L33],

の場合、

θ=[L11L21L31L22L32L33].

になります。

対数コレスキー パラメーター化

コレスキー パラメーター化の L の対角要素が正に制約されている場合、L の解は一意になります。対数コレスキー パラメーター化は、一意のパラメーター化を保証するために L の対角要素の対数を使用するという点を除いて、コレスキー パラメーター化と同じです。

たとえば、コレスキー パラメーター化の 3 行 3 列の例では、次のように強制します。Lii ≥ 0、

θ=[log(L11)L21L31log(L22)L32log(L33)].

この情報は役に立ちましたか?