fitlmematrix

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

構文

lme = fitlmematrix(X,y,Z,[])

lme = fitlmematrix(X,y,Z,G)

lme = fitlmematrix(___,Name,Value)

説明

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

[] は、1 つのグループがあることを意味しています。つまり、グループ化変数 G は ones(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

共変量を使用した長期間の調査

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

load('weight.mat');

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

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

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

$$\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).$$

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:
    y ~ Intercept + InitWeight + PrgB + PrgC + PrgD + Week + Week_PrgB + Week_PrgC + Week_PrgD + (Intercept + 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
    {'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.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

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

ランダムな切片モデル

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

load flu
flu = dataset2table(flu)

flu=52×11 table
         Date          NE      MidAtl    ENCentral    WNCentral    SAtl     ESCentral    WSCentral     Mtn      Pac     WtdILI
    ______________    _____    ______    _________    _________    _____    _________    _________    _____    _____    ______

    {'10/9/2005' }     0.97    1.025       1.232        1.286      1.082      1.457          1.1      0.981    0.971    1.182 
    {'10/16/2005'}    1.136     1.06       1.228        1.286      1.146      1.644        1.123      0.976    0.917     1.22 
    {'10/23/2005'}    1.135    1.172       1.278        1.536      1.274      1.556        1.236      1.102    0.895     1.31 
    {'10/30/2005'}     1.52    1.489       1.576        1.794       1.59      2.252        1.612      1.321    1.082    1.343 
    {'11/6/2005' }    1.365    1.394        1.53        1.825       1.62      2.059        1.471      1.453    1.118    1.586 
    {'11/13/2005'}     1.39    1.477       1.506          1.9      1.683      1.813        1.464      1.388    1.204     1.47 
    {'11/20/2005'}    1.212    1.231       1.295        1.495      1.347      1.794        1.303      1.371    1.137    1.611 
    {'11/27/2005'}    1.477    1.546       1.557        1.855      1.678      2.159        1.739      1.628    1.443    1.827 
    {'12/4/2005' }    1.285     1.43       1.482        1.635      1.577      1.903         1.53      1.701    1.516    1.776 
    {'12/11/2005'}    1.354     1.45        1.46        1.794      1.583      1.894        1.831      2.364    2.094    1.941 
    {'12/18/2005'}    1.502    1.622       1.638        1.988      1.947       2.22        2.577       3.89     2.66     2.34 
    {'12/25/2005'}     1.86    1.915       1.955         2.38      2.343      3.027        3.219      4.862    2.595    3.086 
    {'1/1/2006'  }    2.114    2.174       2.065        2.557      2.275      2.498        2.644      3.352    2.181     3.26 
    {'1/8/2006'  }    1.815    1.932       1.822        2.046      1.969      1.805        2.189      2.132    1.717    2.613 
    {'1/15/2006' }    1.541    1.695       1.581        2.008      1.718      1.662        2.156      1.694    1.351    2.247 
    {'1/22/2006' }    1.632    1.758       1.711        2.217      1.866      2.194        2.268      1.826    1.384    2.352 
      ⋮

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

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

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

ランダムな切片の線形混合効果モデルの計画行列を定義します。ここで、切片は 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).$$

y = flu.FluRate;
X = [ones(468,1) flu.WtdILI];
Z = [ones(468,1)];
G = flu.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

変量効果項の標準偏差 $${\sigma }_b$$ の信頼限界 (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)*flu.WtdILI(28) + STATS.Estimate(STATS.Level=='10/30/2005')

y_hat = 
1.4674

fitted(lme) メソッドを使用して、近似値を簡単に表示できます。

F = fitted(lme);
F(28)

ans = 
1.4674

乱塊法

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

load('shift.mat');

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

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

$$\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).$$

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 = 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

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

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

load carbig

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

$${MPG}_{im}={\beta }_0+{\beta }_1{Acc}_i+{\beta }_{2}HP+b_{0m}+{b_1}_m{Acc}_{im}+{\epsilon }_{im},m=1,2,3,$$

変量効果項の事前分布は次のようになります。

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

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

ここで、$$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.0929e-07    NaN      NaN  

Group: Model_Year (13 Levels)
    Name1                   Name2                   Type           Estimate    Lower      Upper 
    {'Acceleration'}        {'Acceleration'}        {'std'}        0.18783     0.17509    0.2015

Group: Error
    Name               Estimate    Lower     Upper 
    {'Res Std'}        3.7258      3.4731    3.9969

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

切片と加速度について相関する可能性がある変量効果でモデルを当てはめます。この場合、変量効果項の事前分布は次のようになります。

$$b_m=\left(\begin{array}{l}{b}_{0m}\\ b_{1m}\end{array}\right)\sim N\left(0,\left(\begin{array}{cc}{\sigma }_0^2& {\sigma }_{0,1}\\ {\sigma }_{0,1}& {\sigma }_1^2\end{array}\right)\right),$$

ここで、$$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 つのモデルを比較できます。

共分散パターンの指定

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

load('weight.mat');

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;

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

$$\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\\ +b{}_{0m}+b_{1m}Week{2}_{im}+b_{2m}Week{4}_{im}+b_{3m}Week{6}_{im}+b_{4m}Week{8}_{im}\\ +b_{5m}Week10_{im}+b_{6m}Week12_{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[.]$$ はプログラムの種類を表すダミー変数です。たとえば、$$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).$$

定義した計画行列とグループ化変数で 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) 列の行列の cell 配列、r = 1, 2, ..., R

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

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

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

データ型: single | double | cell

G — グループ化変数

n 行 1 列のベクトル | R n 行 1 列のベクトルの cell 配列

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

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

  G は、categorical ベクトル、logical ベクトル、数値ベクトル、文字配列、string 配列、または文字ベクトルの cell 配列が可能です。

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

  G{r} は、categorical ベクトル、logical ベクトル、数値ベクトル、文字配列、string 配列、または文字ベクトルの cell 配列が可能です。

データ型: categorical | logical | single | double | char | string | cell

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

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

FixedEffectPredictors — 固定効果の計画行列の列名

{'x1','x2',...,'xP'} (既定値) | 長さ p の string 配列または cell 配列

固定効果の計画行列 X の列名。'FixedEffectPredictors' と長さ p の string 配列または cell 配列から構成されるコンマ区切りのペアとして指定します。

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

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

データ型: string | cell

RandomEffectPredictors — 変量効果の計画行列の列名または cell 配列

長さ q の string 配列または cell 配列 | 長さ q(r) の要素をもつ長さ R の cell 配列、r = 1, 2, ..., R

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

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

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

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

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

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

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

データ型: string | cell

ResponseVarName — 応答変数の名前

'y' (既定値) | 文字ベクトル | string スカラー

応答変数の名前。'ResponseVarName' と文字ベクトルまたは string スカラーから構成されるコンマ区切りのペアとして指定します。

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

例: 'ResponseVarName','score'

データ型: char | string

RandomEffectGroups — 変量効果のグループ化変数の名前

'g' または {'g1','g2',...,'gR'} (既定値) | 文字ベクトル | string スカラー | string 配列 | 文字ベクトルの cell 配列

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

• 文字ベクトルまたは string スカラー — 変量効果項が 1 つしかない場合、つまり G がベクトルである場合、'RandomEffectGroups' の値はグループ化変数 G の名前です。既定の設定は 'g' です。

• string 配列または文字ベクトルの cell 配列 — 複数の変量効果項がある場合、つまり G が長さ R の cell 配列である場合、'RandomEffectGroups' の値は長さ R の string 配列または cell 配列です。各要素はグループ化変数 G{r} の名前です。既定の設定は {'g1','g2',...,'gR'} です。

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

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

データ型: char | string | cell

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' (既定の設定)	fitlmematrix は、参照グループを使用してダミー変数を作成します。この方式では、最初のカテゴリを参照グループとして扱い、カテゴリの数よりも 1 つ少ないダミー変数を作成します。カテゴリカル変数のカテゴリの順序は、関数 categories を使用してチェックできます。順序を変更するには、関数 reordercats を使用します。
'effects'	fitlmematrix は、エフェクト コーディングを使用してダミー変数を作成します。この方式では、–1 を使用して最後のカテゴリを表します。この方式では、カテゴリの数よりも 1 つ少ないダミー変数を作成します。
'full'	fitlmematrix は、完全なダミー変数を作成します。この方式では、各カテゴリに対して 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 オブジェクトとして返されます。

詳細

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

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

代替機能

線形混合効果モデルの当てはめには、fitlme(tbl,formula) も使用できます。ここで、tbl は応答 y、予測子変数 X、およびグループ化変数が含まれている table であり、formula は 'y ~ fixed + (random1|g1) + ... + (randomR|gR)' という形式の式です。

バージョン履歴

R2013b で導入

参考

fitlme | LinearMixedModel | compare