fitlme

線形混合効果モデルの当てはめ

構文

lme = fitlme(tbl,formula)

lme = fitlme(tbl,formula,Name,Value)

説明

lme = fitlme(tbl,formula) は線形混合効果モデルを返します。formula によって指定されるこのモデルは、table またはデータセット配列 tbl の変数に当てはめられています。

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

たとえば、変量効果の項の共分散パターン、パラメーターの推定に使用するメソッド、最適化アルゴリズムのオプションなどを指定できます。

例

線形混合効果モデルの当てはめ

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

load imports-85

変数を table に格納します。

tbl = table(X(:,12),X(:,14),X(:,24),'VariableNames',{'Horsepower','CityMPG','EngineType'});

table の最初の 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        Lower       Upper    
    {'(Intercept)'}          37.276     2.8556    13.054    201    1.3147e-28      31.645       42.906
    {'Horsepower' }        -0.12631    0.02284     -5.53    201    9.8848e-08    -0.17134    -0.081269

Random effects covariance parameters (95% CIs):
Group: EngineType (7 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        5.7338      2.3773    13.829

Group: EngineType (7 Levels)
    Name1                 Name2                 Type           Estimate    Lower      Upper  
    {'Horsepower'}        {'Horsepower'}        {'std'}        0.050357    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        Lower      Upper    
    {'(Intercept)'}         33.824       4.0181     8.4178    201    7.1678e-15     25.901       41.747
    {'Horsepower' }        -0.1087     0.032912    -3.3029    201     0.0011328    -0.1736    -0.043806

Random effects covariance parameters (95% CIs):
Group: EngineType (7 Levels)
    Name1                  Name2                  Type            Estimate    Lower       Upper   
    {'(Intercept)'}        {'(Intercept)'}        {'std' }          9.4952      4.7022      19.174
    {'Horsepower' }        {'(Intercept)'}        {'corr'}        -0.96843    -0.99568    -0.78738
    {'Horsepower' }        {'Horsepower' }        {'std' }        0.078874    0.039917     0.15585

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        3.1845      2.8774    3.5243

切片と馬力の変量効果共分散パラメーターが同時に表示され、さらに切片と馬力の間の相関 ('corr') も含まれています。

切片がランダムな LME モデルの当てはめ

標本データを読み込み、table の形式に変換します。

load flu
flu = dataset2table(flu);

table flu には、変数 Date と、インフルエンザ推定罹患率 (Google® 検索から推定される 9 地域の値と疾病対策センター (CDC) による全国の推定値) が格納されている 10 個の変数が含まれています。

線形混合効果モデルを当てはめるには、データが適切な形式の table になっていなければなりません。インフルエンザ羅患率を応答として使用することにより線形混合効果モデルを当てはめるには、地域に対応する 9 個の列を 1 つの table 変数にまとめます。新しい table flu2 には、新しい応答変数 FluRate、各推定の元になっている地域を示すノミナル変数 Region、全国の推定値 WtdILI、およびグループ化変数 Date が含まれなければなりません。

flu2 = stack(flu,2:10,NewDataVariableName="FluRate", ...
    IndexVariableName="Region");
flu2.Date = categorical(flu2.Date);

flu2 の最初の 6 行を表示します。

flu2(1:6,:)

ans=6×4 table
      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 で変化するランダム切片で、線形混合効果モデルを当てはめます。このモデルは以下に対応します。

$$y_{im}={\beta }_0+{\beta }_1{WtdILI}_{im}+b_{0m}+{\epsilon }_{im},i=1,2,...,468,m=1,2,...,52,$$

ここで、$$y_{im}$$ はグループ化変数 Date の水準 $$m$$ に対する観測値 $$i$$、$$b_{0m}$$ はグループ化変数 Date の水準 $$m$$ に対する変量効果、$${\epsilon }_{im}$$ は観測値 $$i$$ の観測誤差です。変量効果の事前分布は次のようになります。

$$b_{0m}\sim N(0,{\sigma }_b^2),$$

誤差項の分布は次のようになります。

$${\epsilon }_{im}\sim N(0,{\sigma }^2).$$

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        Lower       Upper  
    {'(Intercept)'}        0.16385     0.057525    2.8484    466     0.0045885    0.050813    0.27689
    {'WtdILI'     }         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

"Random effects covariance parameters" というタイトルのセクションに、推定共分散パラメーターが表示されます。$${\sigma }_b$$ の推定値は 0.17146、95% 信頼区間は [0.13227, 0.22226] です。この間隔には 0 は含まれていないため、変量効果項は有意です。compareメソッドによる尤度比検定を使用すると、任意の変量効果項の有意性を正式に検定できます。

観測値について推定される応答は、その観測値に対応するグループ化変数水準の固定効果および変量効果の値の合計です。たとえば、観測値 28 のインフルエンザ推定羅患率は次のようになります。

$$\begin{array}{l}{\underset{}{\overset{ˆ}{y}}}_{28}={\underset{}{\overset{ˆ}{\beta }}}_0+{\underset{}{\overset{ˆ}{\beta }}}_1{WtdILI}_{28}+{\underset{}{\overset{ˆ}{b}}}_{10/30/2005}\\ =0.1639+0.7236*(1.343)+0.3318\\ =1.46749,\end{array}$$

ここで、$$\underset{}{\overset{ˆ}{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 メソッドを使用して、近似値を表示できます。

F = fitted(lme);
F(28)

ans = 
1.4674

乱塊法用の LME モデル

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

load('shift.mat')

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

シフトの時間によってパフォーマンスに有意差があるかどうかを評価するために、作業者別のランダムな切片をもつ線形混合効果モデルを当てはめます。制限付き最尤法と 'effects' 対比を使用します。

'effects' 対比は係数の合計が 0 であることを意味し、fitlme はシフトの効果を表す、"固定効果の計画行列" と呼ばれる行列を作成します。この行列には、$$Shift\_Evening$$ および $$Shift\_Morning$$ という 2 つの列があります。ここで、

$$\text{Shift}\text{\_}\text{Evening}=\{\begin{array}{ll}0,& \text{if}\text{Morning}\\ 1,& \text{if}\text{Evening}\\ -1,& \text{if}\text{Night}\end{array}$$

$$\text{Shift}\text{\_}\text{Morning}=\{\begin{array}{ll}1,& \text{if}\text{Morning}\\ 0,& \text{if}\text{Evening}\\ -1,& \text{if}\text{Night}\end{array}$$

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

$$\begin{array}{lrl}\text{Morning}\text{Shift:}& {\text{QCDev}}_{im}& ={\beta }_0+{\beta }_2{\text{Shift}\text{\_}\text{Morning}}_i+b_{0m}+{ϵ}_{im},\\ \text{Evening}\text{Shift:}& {\text{QCDev}}_{im}& ={\beta }_0+{\beta }_1{\text{Shift}\text{\_}\text{Evening}}_i+b_{0m}+{ϵ}_{im},\\ \text{Night}\text{Shift:}& {\text{QCDev}}_{im}& ={\beta }_0-{\beta }_1{\text{Shift}\text{\_}\text{Evening}}_i-{\beta }_2{\text{Shift}\text{\_}\text{Morning}}_i+b_{0m}+{ϵ}_{im},\end{array}$$

ここで、$$i$$ は観測値を、$$m$$ は作業者を表します。$$i$$ = 1、2、...、15 および $$m$$ = 1、2、...、5 です。変量効果と観測値の誤差の分布は次のとおりです。

$$b_{0m}\sim N(0,{\sigma }_b^2)$$

および

$${\epsilon }_{im}\sim N(0,{\sigma }^2).$$

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       Lower      Upper   
    {'(Intercept)'  }          3.6525    0.94109     3.8812    12    0.0021832     1.6021       5.703
    {'Shift_Evening'}        -0.53293    0.31206    -1.7078    12      0.11339    -1.2129     0.14699
    {'Shift_Morning'}        -0.91973    0.31206    -2.9473    12     0.012206    -1.5997    -0.23981

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 = 5×1

    0.5775
    1.1757
   -2.1715
    2.3655
   -1.9472

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

$$\begin{array}{l}{\underset{}{\overset{ˆ}{y}}}_{\text{Evening},\text{Operator}3}={\underset{}{\overset{ˆ}{\beta }}}_0+{\underset{}{\overset{ˆ}{\beta }}}_1\text{Shift}\text{\_}\text{Evening}+{\underset{}{\overset{ˆ}{b}}}_{03}\\ =3.6525-0.53293-2.1715\\ =0.94807.\end{array}$$

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

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

ans = 
0.9481

同様に、朝シフトで作業している 3 番目の作業者の場合、品質目標の特性の推定絶対偏差を以下のように計算できます。

$$\begin{array}{l}{\underset{}{\overset{ˆ}{y}}}_{\text{Morning},\text{Operator}3}={\underset{}{\overset{ˆ}{\beta }}}_0+{\underset{}{\overset{ˆ}{\beta }}}_2\text{Shift}\text{\_}\text{Morning}+{\underset{}{\overset{ˆ}{b}}}_{03}\\ =3.6525-0.91973-2.1715\\ =0.56127.\end{array}$$

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

F(shift.Shift=='Morning' & shift.Operator=='3')

ans = 
0.5613

この作業者は、朝シフトのときに誤差の大きさが減少する傾向にあります。

分割プロット試験用の LME モデル

標本データを読み込んで表示します。

load fertilizer.mat
tbl

tbl=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 は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。

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

$$\begin{array}{l}{y}_{imjk}={\beta }_0+\sum _{m=2}^4{\beta }_{1m}I{\left[F\right]}_{im}+\sum _{j=2}^5{\beta }_{2j}I{\left[T\right]}_{ij}+\sum _{j=2}^5\sum _{m=2}^4{\beta }_{3mj}I{\left[F\right]}_{im}I{\left[T\right]}_{ij}\\ +b_{0k}{S}_k+b_{0jk}(S*T{)}_{jk}+{\epsilon }_{imjk},\end{array}$$

ここで $$i$$ = 1、2、...、60 で、インデックス $$m$$ は肥料の種類、$$j$$ はトマトの種類に対応し、$$k$$ = 1、2、3 はブロック (土壌) に対応します。$$S_k$$ は $$k$$ 番目の土壌の種類を表し、$$(S*T{)}_{jk}$$ は $$k$$ 番目の土壌の種類で入れ子にされている $$j$$ 番目のトマトの種類を表します。$$I[F{]}_{im}$$ は肥料の水準 $$m$$ を表すダミー変数です。同様に、$$I[T{]}_{ij}$$ は、トマトの種類の水準 $$j$$ を表すダミー変数です。

変量効果と観測誤差の事前分布は、$$b_{0k}$$~N(0, $${\sigma }_S^2$$ )、$$b_{0jk}$$~N(0, $${\sigma }_{S*T}^2$$ )、および $${ϵ}_{imjk}$$ ~ N(0, $${\sigma }^2$$) です。

lme = fitlme(tbl,"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    pValue        Lower      Upper 
    {'(Intercept)'                 }             77     8.5836      8.9706    40    4.0206e-11     59.652    94.348
    {'Tomato_Grape'                }            -16     11.966     -1.3371    40       0.18873    -40.184    8.1837
    {'Tomato_Heirloom'             }        -6.6667     11.966    -0.55714    40       0.58053     -30.85    17.517
    {'Tomato_Plum'                 }         32.333     11.966      2.7022    40      0.010059     8.1496    56.517
    {'Tomato_Vine'                 }            -13     11.966     -1.0864    40       0.28379    -37.184    11.184
    {'Fertilizer_2'                }         34.667      8.572      4.0442    40    0.00023272     17.342    51.991
    {'Fertilizer_3'                }         33.667      8.572      3.9275    40    0.00033057     16.342    50.991
    {'Fertilizer_4'                }         47.667      8.572      5.5607    40    1.9567e-06     30.342    64.991
    {'Tomato_Grape:Fertilizer_2'   }        -2.6667     12.123    -0.21997    40       0.82701    -27.167    21.834
    {'Tomato_Heirloom:Fertilizer_2'}             -8     12.123    -0.65992    40       0.51309    -32.501    16.501
    {'Tomato_Plum:Fertilizer_2'    }            -15     12.123     -1.2374    40       0.22317    -39.501    9.5007
    {'Tomato_Vine:Fertilizer_2'    }            -16     12.123     -1.3198    40       0.19439    -40.501    8.5007
    {'Tomato_Grape:Fertilizer_3'   }         16.667     12.123      1.3748    40       0.17683    -7.8341    41.167
    {'Tomato_Heirloom:Fertilizer_3'}         3.3333     12.123     0.27497    40       0.78476    -21.167    27.834
    {'Tomato_Plum:Fertilizer_3'    }         3.6667     12.123     0.30246    40       0.76387    -20.834    28.167
    {'Tomato_Vine:Fertilizer_3'    }              3     12.123     0.24747    40       0.80581    -21.501    27.501
    {'Tomato_Grape:Fertilizer_4'   }         13.333     12.123      1.0999    40       0.27796    -11.167    37.834
    {'Tomato_Heirloom:Fertilizer_4'}            -19     12.123     -1.5673    40       0.12492    -43.501    5.5007
    {'Tomato_Plum:Fertilizer_4'    }        -2.6667     12.123    -0.21997    40       0.82701    -27.167    21.834
    {'Tomato_Vine:Fertilizer_4'    }         8.6667     12.123     0.71492    40       0.47881    -15.834    33.167

Random effects covariance parameters (95% CIs):
Group: Soil (3 Levels)
    Name1                  Name2                  Type           Estimate    Lower       Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        2.5028      0.027711    226.04

Group: Soil:Tomato (15 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        10.225      6.1497    17.001

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        10.499      8.5389    12.908

固定効果係数の表示の最後の 12 行に対応する $$p$$ 値 (0.82701 ～ 0.47881) は、トマトの種類と肥料の種類の間の交互作用係数が有意ではないことを示しています。トマトと肥料の間の全体的な交互作用を検定するには、effects 対比を使用してモデルを再度当てはめた後でanovaメソッドを使用します。

変量効果項の標準偏差 ($${\sigma }_S^2$$) の信頼区間は非常に大きくなっています。切片は土壌でグループ化されています。この項には有意性がありません。

交互作用項 Tomato:Fertilizer および変量効果の項 (1 | Soil) を削除した後、モデルを再度当てはめます。

lme = fitlme(tbl,"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    pValue        Lower      Upper 
    {'(Intercept)'    }         77.733     7.3293      10.606    52    1.3108e-14     63.026    92.441
    {'Tomato_Grape'   }        -9.1667     9.6045    -0.95441    52       0.34429    -28.439    10.106
    {'Tomato_Heirloom'}        -12.583     9.6045     -1.3102    52        0.1959    -31.856    6.6895
    {'Tomato_Plum'    }         28.833     9.6045      3.0021    52     0.0041138     9.5605    48.106
    {'Tomato_Vine'    }        -14.083     9.6045     -1.4663    52       0.14858    -33.356    5.1895
    {'Fertilizer_2'   }         26.333     4.5004      5.8514    52    3.3024e-07     17.303    35.364
    {'Fertilizer_3'   }             39     4.5004      8.6659    52    1.1459e-11     29.969    48.031
    {'Fertilizer_4'   }         47.733     4.5004      10.607    52     1.308e-14     38.703    56.764

Random effects covariance parameters (95% CIs):
Group: Soil:Tomato (15 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        10.02       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 週間の期間にわたって記録されています。このデータは、シミュレーションされたものです。

データを table に保存します。Subject および Program をカテゴリカル変数として定義します。

tbl = table(InitialWeight,Program,Subject,Week,y);
tbl.Subject = nominal(tbl.Subject);
tbl.Program = nominal(tbl.Program);

線形混合効果モデルを当てはめます。初期体重、プログラムの種類、週、週とプログラムの種類の間の交互作用は固定効果です。切片と週は被験者ごとに異なります。

fitlme はプログラム A を基準として使用して、必要なダミー変数 $$I$$[.] を作成します。モデルには既に切片があるので、fitlme は、プログラム B、C、D に対してダミー変数のみを作成します。これは、ダミー変数をコーディングする 'reference' メソッドとも呼ばれます。

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

$$\begin{array}{l}{y}_{im}={\beta }_0+{\beta }_{1}I{W}_i+{\beta }_{2}Wee{k}_i+{\beta }_{3}I{\left[PB\right]}_i+{\beta }_{4}I{\left[PC\right]}_i+{\beta }_{5}I{\left[PD\right]}_i\\ +{\beta }_6\left(Wee{k}_i*I{\left[PB\right]}_i\right)+{\beta }_7\left(Wee{k}_i*I{\left[PC\right]}_i\right)+{\beta }_8\left(Wee{k}_i*I{\left[PD\right]}_i\right)\\ +b{}_{0m}+b_{1m}Wee{k}_{im}+{\epsilon }_{im},\end{array}$$

ここで、$$i$$ = 1、2、...、120、$$m$$ = 1、2、...、20 です。$${\beta }_j$$ は固定効果係数 ($$j$$ = 0、1、...、8)、$$b_{0m}$$ と $$b_{1m}$$ は変量効果です。$$IW$$ は初期体重を表し、$$I[\cdot ]$$ はプログラムの種類を表すダミー変数です。たとえば、$$I[PB{]}_i$$ はプログラムの種類 B を表すダミー変数です。変量効果と観測誤差の事前分布は次のとおりです。

$$b_{0m}\sim N(0,{\sigma }_0^2)$$

$$b_{1m}\sim N(0,{\sigma }_1^2)$$

$${\epsilon }_{im}\sim N(0,{\sigma }^2).$$

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     pValue       Lower         Upper    
    {'(Intercept)'   }          0.66105      0.25892      2.5531    111     0.012034       0.14798       1.1741
    {'InitialWeight' }        0.0031879    0.0013814      2.3078    111     0.022863    0.00045067    0.0059252
    {'Program_B'     }          0.36079      0.13139       2.746    111    0.0070394       0.10044      0.62113
    {'Program_C'     }        -0.033263      0.13117    -0.25358    111      0.80029      -0.29319      0.22666
    {'Program_D'     }          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
    {'Program_B:Week'}         0.038771     0.095394     0.40644    111      0.68521      -0.15026       0.2278
    {'Program_C:Week'}         0.030543     0.095394     0.32018    111      0.74944      -0.15849      0.21957
    {'Program_D:Week'}         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.21077    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.022863 および 0.011567 は、被験者の初期体重と時間が体重の減少量に対して有意な影響を与えることを示しています。プログラム B の被験者の体重の減少は、プログラム A の被験者の体重の減少に対して有意差があります。変量効果の共分散パラメーターの下限または上限に 0 は含まれないため、この差は有意です。compare メソッドを使用して、変量効果の有意性をテストすることもできます。

入力引数

tbl — 入力データ

table | dataset 配列

入力データ。応答変数、予測子変数およびグループ化変数を含みます。table または dataset 配列として指定します。予測子変数は、連続変数またはグループ化変数にすることができます (グループ化変数を参照してください)。応答変数は、数値でなければなりません。formula を使用して、変数のモデルを指定しなければなりません。

データ型: table

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 に制約されます。 $$\left(\begin{array}{ccc}{\sigma }_{b1}^2& 0& 0\\ 0& {\sigma }_{b2}^2& 0\\ 0& 0& {\sigma }_{b3}^2\end{array}\right)$$
'Isotropic'	分散が等しい対角共分散行列。つまり、共分散行列の非対角要素は 0 に制約され、対角要素は等価に制約されます。たとえば、等方性共分散構造をもつ変量効果の項が 3 つある場合、この共分散行列は次のようになります。 $$\left(\begin{array}{ccc}{\sigma }_b^2& 0& 0\\ 0& {\sigma }_b^2& 0\\ 0& 0& {\sigma }_b^2\end{array}\right)$$ σ2b は、変量効果項の共通分散です。
'CompSymm'	複合対称構造。つまり、対角線上の共通分散とすべての変量効果間の等しい相関です。たとえば、複合対称構造の共分散行列をもつ変量効果の項が 3 つある場合、この共分散行列は次のようになります。 $$\left(\begin{array}{ccc}{\sigma }_{b1}^2& {\sigma }_{b1,b2}& {\sigma }_{b1,b2}\\ {\sigma }_{b1,b2}& {\sigma }_{b1}^2& {\sigma }_{b1,b2}\\ {\sigma }_{b1,b2}& {\sigma }_{b1,b2}& {\sigma }_{b1}^2\end{array}\right)$$ σ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 には固定効果および変量効果の項が含まれます。

table tbl に以下のものが格納されていると仮定します。

• 応答変数 y

• 連続変数またはグループ化変数である予測子変数 Xj

• グループ化変数 g1、g2、...、gR

ここで、Xj および gr のグループ化変数は、categorical 配列、logical 配列、文字配列、string 配列、または文字ベクトルの cell 配列が可能です。

この場合、'y ~ fixed + (random1|g1) + ... + (randomR|gR)' の形式の式において、項 fixed は固定効果の計画行列 X の仕様に対応し、random₁ はグループ化変数 g₁ に対応する変量効果の計画行列 Z₁ の仕様であり、同様に 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 つに、変量効果の事前分布が次のようになっているというものがあります。

$$b~N\left(0,{\sigma }^{2}D\left(\theta \right)\right),$$

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

$$D\left(\theta \right)=L\left(\theta \right)L{\left(\theta \right)}^T,$$

すると、制約がない q*(q+1)/2 行 1 列のパラメーター ベクトル θ が L の下三角部分にある要素から形成されます。

たとえば、次のようになっているとします。

$$L=\left[\begin{array}{ccc}{L}_{11}& 0& 0\\ L_{21}& L_{22}& 0\\ L_{31}& L_{32}& L_{33}\end{array}\right],$$

この場合、次のようになります。

$$\theta =\left[\begin{array}{c}{L}_{11}\\ L_{21}\\ L_{31}\\ L_{22}\\ L_{32}\\ L_{33}\end{array}\right].$$

対数コレスキー パラメーター表現

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

たとえば、3 行 3 列のコレスキー パラメーター表現の例では、L_ii ≥ 0 を適用すると次のようになります。

$$\theta =\left[\begin{array}{c}\log \left(L_{11}\right)\\ L_{21}\\ L_{31}\\ \log \left(L_{22}\right)\\ L_{32}\\ \log \left(L_{33}\right)\end{array}\right].$$

代替方法

式で表現することが困難なモデルの場合は、固定効果と変量効果を定義する行列を作成し、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 で導入

参考

LinearMixedModel | fitlmematrix

トピック

• 式と計画行列の関係