ドキュメンテーション

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

fitlme

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

構文

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

説明

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

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

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

すべて折りたたむ

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

load imports-85

変数をテーブルに格納します。

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

テーブルの最初の 5 行を表示します。

tbl(1:5,:)
ans = 

    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    Lower 
    '(Intercept)'        '(Intercept)'        'std'        5.7338      2.3773


    Upper 
    13.829

Group: EngineType (7 Levels)
    Name1               Name2               Type         Estimate    Lower  
    'Horsepower'        'Horsepower'        'std'        0.050357    0.02307


    Upper  
    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    Lower   
    '(Intercept)'        '(Intercept)'        'std'           9.4952      4.7022
    'Horsepower'         '(Intercept)'        'corr'        -0.96843    -0.99568
    'Horsepower'         'Horsepower'         'std'         0.078874    0.039917


    Upper   
      19.174
    -0.78738
     0.15585

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

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

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

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 で変化するランダム切片で、線形混合効果モデルを近似します。このモデルは以下の式に対応します。

$${y_{im}} = {\beta _0} + {\beta _1}WtdIL{I_{im}} + {b_{0m}} +
{\varepsilon _{im}},\quad i = 1,2,...,468,\quad m = 1,2,...,52,$$

ここで、 $y_{im}$ はグループ化変数 Date のレベル $m$ の観測値 $i$ です。 $b_{0m}$ はグループ化変数 Date のレベル $m$ の変量効果、 $\epsilon_{im}$ は観測値 $i$ の観測誤差です。変量効果の事前分布は $b$ ~ N(0, $\sigma^{2}_{b}$ )、誤差項の分布は $\epsilon$ ~ 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    
    '(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    Lower  
    '(Intercept)'        '(Intercept)'        'std'        0.17146     0.13227


    Upper  
    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 は含まれていないため、変量効果項は有意です。LinearMixedModel.compare メソッドによる尤度比検定を使用すると、任意の変量効果項の有意性を正式に検定できます。

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

$$\begin{array}{l}
{{\hat y}_{28}} = {{\hat \beta }_0} + {{\hat \beta }_1}WtdIL{I_{28}} + {{\hat b}_{10/30/2005}}\\
\quad \quad = 0.1639 + 0.7236*(1.343) + 0.3318\\
\quad \quad = 1.46749,
\end{array}$$

ここで、 ${\hat 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

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

load(fullfile(matlabroot,'examples','stats','shift.mat'))

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

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

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

$$Shift{\rm{\_}}Evening = \left\{ {\begin{array}{*{20}{c}}
{0,\quad if\;Morning}\\
{1,\quad if\;Evening}\\
{ -1,\quad if\;Night}
\end{array}}\right.$$

および

$$Shift{\rm{\_}}Morning = \left\{ {\begin{array}{*{20}{c}}
{1,\quad if\;Morning}\\
{0,\quad if\;Evening}\\
{ - 1,\quad if\;Night }
\end{array}}\right..$$

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

$$\begin{array}{l}
{\rm{Morning Shift: }}QCDe{v_{im}} = {\beta _0} + {\beta _2}Shift{\rm{\_}}Mornin{g_i} + {b_{0m}} + {\varepsilon _{im}},\quad m = 1,2,...,5,\\
{\rm{Evening Shift: }}QCDe{v_{im}} = {\beta _0} + {\beta _1}Shift{\rm{\_}}Evenin{g_i} + {b_{0m}} + {\varepsilon _{im}},\\
{\rm{Night Shift:  }}\quad QCDe{v_{im}} = {\beta _0} - {\beta _1}Shift{\rm{\_}}Evenin{g_i} - {\beta _2}Shift{\rm{\_}}Mornin{g_i} + {b_{0m}} + {\varepsilon _{im}},
\end{array}$$

ここで、 $b$ ~ N(0, $\sigma^{2}_{b}$ ) および $\epsilon$ ~ 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   
    '(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    Lower  
    '(Intercept)'        '(Intercept)'        'std'        2.0457      0.98207


    Upper 
    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 番目の作業者の場合、品質目標の特性の推定絶対偏差は以下のとおりです。

$$ \begin{array}{l}
{{\hat y}_{Evening,Operator3}} = {{\hat \beta }_0} + {{\hat \beta }_1}Shift{\rm{\_}}Evening + {{\hat b}_{03}}\\
\quad \quad \quad \quad \quad \quad \quad \; = 3.6525 - 0.53293 - 2.1715\\
\quad \quad \quad \quad \quad \quad \quad \; = 0.94807.
\end{array}$$

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

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

    0.9481

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

$$ \begin{array}{l}
{{\hat y}_{Morning,Operator3}} = {{\hat \beta }_0} + {{\hat \beta }_2}Shift{\rm{\_}}Morning + {{\hat b}_{03}}\\
\quad \quad \quad \quad \quad \quad \quad \; = 3.6525 - 0.91973 - 2.1715\\
\quad \quad \quad \quad \quad \quad \quad \; = 0.56127.
\end{array} $$

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

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

    0.5613

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

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

load(fullfile(matlabroot,'examples','stats','fertilizer.mat'))

このデータセット配列には土壌の種類に基づいて土壌が 3 つのブロックに分けられている分割プロット試験のデータが含まれています。土壌の種類は砂質、シルトおよび粘土質です。各ブロックは 5 つのプロットに分割され、5 種類のトマトの苗木 (チェリー、エアルーム、グレープ、枝付き、プラム) がランダムにこれらのプロットに割り当てられます。その後、プロット内のトマトの苗木はサブプロットに分割され、それぞれのサブプロットが 4 つの肥料の中の 1 つにより処置されます。このデータは、シミュレーションされたものです。

このデータを ds という名前のデータセット配列に保存し、TomatoSoil および Fertilizer をカテゴリカル変数として定義します。

ds = fertilizer;
ds.Tomato = nominal(ds.Tomato);
ds.Soil = nominal(ds.Soil);
ds.Fertilizer = nominal(ds.Fertilizer);

線形混合効果モデルを近似します。Fertilizer および Tomato は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。

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

$$\begin{array}{l}
{y_{imjk}} = {\beta _0} + \sum\limits_{m = 2}^4 {{\beta _{1m}}I{{\left[ F \right]}_{im}}}  + \sum\limits_{j = 2}^5 {{\beta _{2j}}I{{\left[ T \right]}_{ij}}}  + \sum\limits_{j = 2}^5 {\sum\limits_{m = 2}^4 {{\beta _{3mj}}I{{\left[ F \right]}_{im}}I{{\left[ T \right]}_{ij}}} } \\
\quad \quad \quad  + {b_{0k}}{S_k} + {b_{0jk}}{(S*T)_{jk}} + {\varepsilon _{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^{2}_{S}$ )、 $b_{0jk}$~N(0, $\sigma^{2}_{S*T}$ ) および $\epsilon_{imjk}$ ~ N(0, $\sigma^{2}$ )

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    Lower   
    '(Intercept)'        '(Intercept)'        'std'        2.5028      0.027711


    Upper 
    226.05

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


    Upper 
    17.001

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

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

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

交互作用項 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    pValue    
    '(Intercept)'             77.733     7.3293      10.606    52    1.3108e-14
    'Tomato_Grape'           -9.1667     9.6045    -0.95441    52       0.34429
    'Tomato_Heirloom'        -12.583     9.6045     -1.3102    52        0.1959
    'Tomato_Plum'             28.833     9.6045      3.0021    52     0.0041138
    'Tomato_Vine'            -14.083     9.6045     -1.4663    52       0.14858
    'Fertilizer_2'            26.333     4.5004      5.8514    52    3.3024e-07
    'Fertilizer_3'                39     4.5004      8.6659    52    1.1459e-11
    'Fertilizer_4'            47.733     4.5004      10.607    52     1.308e-14


    Lower      Upper 
     63.026    92.441
    -28.439    10.106
    -31.856    6.6895
     9.5605    48.106
    -33.356    5.1895
     17.303    35.364
     29.969    48.031
     38.703    56.764

Random effects covariance parameters (95% CIs):
Group: Soil:Tomato (15 Levels)
    Name1                Name2                Type         Estimate    Lower 
    '(Intercept)'        '(Intercept)'        'std'        10.02       6.0812


    Upper 
    16.509

Group: Error
    Name             Estimate    Lower     Upper 
    'Res Std'        12.325      10.024    15.153

固定効果および変量効果の項がテストされるため、シミュレーションされた尤度比検定で LinearMixedModel.compare メソッドを使用して 2 つのモデルを比較できます。

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

load(fullfile(matlabroot,'examples','stats','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 を基準として使用して、必要なダミー変数 $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}\\
\quad \quad  + {\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)\\
\quad \quad  + b{}_{0m} + \,{b_{1m}}Wee{k_{im}} + {\varepsilon _{im}},
\end{array}$$

ここで、 $i$ = 1、2、...、120 および $m$ = 1、2、...、20 です。 $\beta_{j}$ は固定効果係数 ( $j$ = 0、1、...、8)、 $b_{1m}$ $b_{1m}$ は変量効果です。 $IW$ は初期体重を表し、I[.] はプログラムの種類を表すダミー変数です。たとえば、 $I[PB]_{i}$ はプログラム B を表すダミー変数です。変量効果と観測誤差の事前分布は次のようになります。 $b_{0m}$ ~ N(0, $\sigma^{2}_{0}$)、 $b_{1m}$ ~ N(0, $\sigma^{2}_{1}$) および $\epsilon_{im}$ ~ 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 
    '(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    Lower  
    '(Intercept)'        '(Intercept)'        'std'         0.18407     0.12281
    'Week'               '(Intercept)'        'corr'        0.66841     0.21076
    'Week'               'Week'               'std'         0.15033     0.11004


    Upper  
    0.27587
    0.88573
    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 メソッドを使用して、変量効果の有意性をテストすることもできます。

入力引数

すべて折りたたむ

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

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

モデル仕様の式。'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)' の形式の文字列として指定します。この文字列では大文字小文字を区別します。詳細は、を参照してください。

例: 'y ~ treatment + (1|block)'

名前/値のペアの引数

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

例: 'CovariancePattern','Diagonal','Optimizer','fminunc','OptimizerOptions',opt はモデルを指定します。ここで、変量効果項は対角共分散行列の構造体をもち、fitlme は、変数 opt で定義されたカスタム最適化パラメーターを使った fminunc 最適化アルゴリズムを使用します。

すべて折りたたむ

変量効果の共分散行列のパターン。'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','REML'

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

データ型: single | double

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

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

例: 'Exclude',[13,67]

データ型: single | double | logical

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

説明
'reference'既定の設定。0 に設定された最初のカテゴリの係数。
'effects'合計 0 の係数。
'full'各カテゴリに対して 1 つのダミー変数。

例: 'DummyVarCoding','effects'

最適化アルゴリズム。'Optimizer' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

'quasinewton'既定の設定。信頼領域ベースの準ニュートン オプティマイザ―を使用します。statset('LinearMixedModel') を使用して、アルゴリズムのオプションを変更します。オプションを指定しない場合、LinearMixedModel は、statset('LinearMixedModel') の既定のオプションを使用します。
'fminunc'このオプションを指定するには、Optimization Toolbox™ がなければなりません。optimoptions('fminunc') を使用して、アルゴリズムのオプションを変更します。オプションを指定しない場合、LinearMixedModel は、'Algorithm''quasi-newton' に設定された optimoptions('fminunc') の既定のオプションを使用します。

例: 'Optimizer','fminunc'

最適化アルゴリズムのオプション。'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') によって作成された構造体の以下のフィールドを使用します。

目的関数の勾配の相対許容誤差。正のスカラー値として指定します。

ステップ サイズの絶対許容誤差。正のスカラー値として指定します。

許容される最大反復回数。正のスカラー値として指定します。

表示のレベル。'off''iter''final' のいずれかとして指定します。

反復最適化を開始するメソッド。'StartMethod' と以下のいずれかで構成されるコンマ区切りのペアとして指定します。

説明
'default'内部で定義される既定値
'random'ランダムな初期値

例: 'StartMethod','random'

画面に最適化プロセスを表示するインジケーター。'Verbose' と、false または true で構成されるコンマ区切りのペアとして指定します。既定値は false です。

'Verbose' の設定は、'OptimizerOptions''Display' フィールドをオーバーライドします。

例: 'Verbose',true

収束における制約のないパラメーターに関して、目的関数のヘッシアンの正定性をチェックするインジケーター。'CheckHessian' と、false または true で構成されるコンマ区切りのペアとして指定します。既定値は false です。

解の最適性を確認する場合や、共分散パラメーターの数に関してモデルのパラメーターが多すぎるかどうかを判定する場合、'CheckHessian' として true を指定します。

例: 'CheckHessian',true

出力引数

すべて折りたたむ

LinearMixedModel オブジェクトとして返される線形混合効果モデル。

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

代替方法

式で表現することが困難なモデルの場合は、固定効果と変量効果を定義する行列を作成し、fitlmematrix(X,y,Z,G) を使用してモデルを近似させることができます。

詳細

すべて折りたたむ

一般に、モデル仕様の式は 'y ~ terms' の形式の文字列です。線形混合効果モデルでは、この式は 'y ~ fixed + (random1|grouping1) + ... + (randomR|groupingR)' の形式になります。ここで、fixed および random には固定効果および変量効果の項が含まれます。

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

  • 応答変数 y

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

  • グループ化変数 g1g2、...、gR

ここで、Xj および gr のグループ化変数は、カテゴリカル配列、論理配列、文字配列または文字列のセル配列です。

この場合、'y ~ fixed + (random1|g1) + ... + (randomR|gR)' の形式の式において、項 fixed は固定効果の計画行列 X の仕様に対応し、random1 はグループ化変数 g1 に対応する変量効果の計画行列 Z1 の仕様であり、同様に randomR はグループ化変数 gR に対応する変量効果の計画行列 ZR の仕様です。fixed 項および random 項はウィルキンソンの表記法で表現できます。

ウィルキンソンの表記法は、モデルに存在する因子を記述します。この表記法は、モデルに存在する因子に関係するものであり、それらの因子の乗数 (係数) に関係するものではありません。

ウィルキンソンの表記法標準表記の因子
1定数 (切片) 項
X^kk は正の整数XX2、...、Xk
X1 + X2X1X2
X1*X2X1X2X1.*X2 (X1 および X2 を要素ごとに乗算)
X1:X2X1.*X2 のみ
- X2X2 は含めない
X1*X2 + X3X1X2X3X1*X2
X1 + X2 + X3 + X1:X2X1X2X3X1*X2
X1*X2*X3 - X1:X2:X3X1X2X3X1*X2X1*X3X2*X3
X1*(X2 + X3)X1X2X3X1*X2X1*X3

Statistics and Machine Learning Toolbox™ 表記は、-1 を使用して項を明示的に削除しない限り、常に定数項を含みます。線形混合効果モデルの仕様例を次にいくつか挙げます。

次に例を示します。

説明
'y ~ X1 + X2'切片 X1 および X2 の固定効果。これは、'y ~ 1 + X1 + X2' と等価です。
'y ~ -1 + X1 + X2'X1X2 の切片と固定効果はありません。-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)'g1g2 に対する独立したメイン効果のあるランダムな切片モデル + 独立した交互作用効果。

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

線形混合効果モデルの仮定の 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)].

参照

[1] Pinherio, J. C., and D. M. Bates. “Unconstrained Parametrizations for Variance-Covariance Matrices”. Statistics and Computing, Vol. 6, 1996, pp. 289–296.

R2013b で導入

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