covarianceParameters

線形混合効果モデルの共分散パラメーターの抽出

構文

psi = covarianceParameters(lme)

[psi,mse] = covarianceParameters(lme)

[psi,mse,stats] = covarianceParameters(lme)

[psi,mse,stats] = covarianceParameters(lme,Alpha=alpha)

説明

psi = covarianceParameters(lme) は、変動効果の前の共分散をパラメーター表現する、推定された共分散パラメーターを返します。

[psi,mse] = covarianceParameters(lme) は、残差分散の推定も返します。

[psi,mse,stats] = covarianceParameters(lme) は、共分散パラメーターおよび関連する統計を含む cell 配列 stats も返します。

[psi,mse,stats] = covarianceParameters(lme,Alpha=alpha) は、共分散パラメーターの信頼限界の有意水準を指定します。

例

切片に対する 2 つの変量効果

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

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

$$y_{ijk}={\beta }_0+\sum _{j=2}^5{\beta }_{2j}I[T{]}_{ij}+b_{0jk}(S*T{)}_{jk}+{ϵ}_{ijk},$$

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

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

lme = fitlme(tbl,"Yield ~ Fertilizer + (1|Soil) + (1|Soil:Tomato)");

変量効果項の共分散パラメーターの推定値 ($${\sigma }_S^2$$ と $${\sigma }_{S*T}^2$$ の推定値) を計算します。

psi = covarianceParameters(lme)

psi=2×1 cell array
    {[2.7730e-17]}
    {[  352.8481]}

残差分散 ($${\sigma }^2$$) を計算します。

[~,mse] = covarianceParameters(lme)

mse = 
151.9007

相関する可能性がある変量効果の項

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

load weight.mat;

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

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

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

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

ダミー変数の符号化が "reference" である場合、fitlme はプログラム A を基準として使用し、必要なダミー変数 $$I[.]$$ を作成します。このモデルは以下に対応します。

$$\begin{array}{rrl}{y}_{im}& =& {\beta }_0+{\beta }_{1}I{W}_i+{\beta }_2{\text{Week}}_i+{\beta }_{3}I[PB{]}_I+{\beta }_{4}I[PC{]}_i\\ & & +{\beta }_{5}I[PD{]}_i+b_{0m}+b_{1m}{\text{Week}}_{im}+{ϵ}_{im}\end{array}$$

ここで、$$i$$ は観測番号 $$i=1,2,...,120$$ に、$$m$$ は被験者番号 $$m=1,2,...,20$$ に対応します。$${\beta }_j$$ は固定効果係数、$$j=0,1,...,8$$、および $$b_{0m}$$、$$b_{1m}$$ は変量効果です。$$IW$$ は初期重量を意味し、$$I[.]$$ はプログラムの種類を表すダミー変数です。たとえば、$$I[PB{]}_i$$ はプログラム B を表すダミー変数です。

変量効果と観測値の誤差の事前分布は次のとおりです。

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

および

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

lme = fitlme(tbl,"y ~ InitialWeight + Program + (Week|Subject)");

変量効果の共分散パラメーターの推定を計算します。

[psi,mse,stats] = covarianceParameters(lme)

psi = 1×1 cell array
    {2×2 double}

mse = 
0.0105

stats=2×1 cell array
    {3×7 classreg.regr.lmeutils.titleddataset}
    {1×5 classreg.regr.lmeutils.titleddataset}

mse は推定された残差分散です。これは、$${\sigma }^2$$ の推定値です。

変量効果項 ($${\sigma }_0^2$$、$${\sigma }_1^2$$ および $${\sigma }_{0,1}$$) に対する共分散パラメーターの推定値を確認するには、psi のインデックスを指定します。

psi{1}

ans = 2×2

    0.0572    0.0490
    0.0490    0.0624

切片に対する変量効果項の分散の推定値 $${\sigma }_0^2$$ は 0.0572 です。週に対する変量効果項の分散の推定値 $${\sigma }_1^2$$ は 0.0624 です。切片と週に対する変量効果項の共分散の推定値 $${\sigma }_{0,1}$$ は 0.0490 です。

stats は 2 行 1 列の cell 配列です。stats の最初の cell には、変量効果の標準偏差の信頼区間と、切片と週に関する変量効果間の相関が含まれます。これらを表示するには、stats のインデックスを指定します。

stats{1}

ans = 
    Covariance Type: FullCholesky

    Group      Name1                  Name2                  Type            Estimate    Lower      Upper  
    Subject    {'(Intercept)'}        {'(Intercept)'}        {'std' }        0.23927     0.14364    0.39854
    Subject    {'Week'       }        {'(Intercept)'}        {'corr'}        0.81971     0.38662    0.95658
    Subject    {'Week'       }        {'Week'       }        {'std' }         0.2497     0.18303    0.34067

表示では、グループ パラメーター (Group)、変量効果変数 (Name1、Name2)、共分散パラメーターの種類 (Type)、各パラメーターの推定 (Estimate)、パラメーターの 95% 信頼区間 (Lower, Upper) が示されています。この表の推定は、次のようにして psi の推定に関連付けられます。

切片の変量効果の項の標準偏差は、0.23927 = sqrt(0.0527) です。同様に、週の変量効果の項の標準偏差は、0.2497 = sqrt(0.0624) です。最終的に、切片と週の変量効果の項の間にある相関は、0.81971 = 0.0490/(0.23927*0.2497) です。

この表示では、モデルを当てはめるときに使用された共分散パターンも示されていることに注意してください。この例では、共分散パターンは FullCholesky です。変量効果の項の共分散パターンを変更するには、モデルを当てはめるときに名前と値のペアの引数 CovariancePattern を使用しなければなりません。

stats の 2 番目の cell には、残差標準偏差に関する同様の統計が含まれます。2 番目の cell の内容を表示します。

stats{2}

ans = 
    Group    Name               Estimate    Lower       Upper  
    Error    {'Res Std'}        0.10261     0.087882    0.11981

残差標準偏差の推定は、mse の平方根 0.10261 = sqrt(0.0105) です。

2 つのグループ化変数

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

load carbig

変数 MPG、Acceleration、Weight、Model_Year、および Origin には、自動車の燃費、加速度、重量、製造年、および自動車が製造された国が格納されています。

MPG を応答変数、Acceleration と Weight を固定効果として使用して線形混合効果モデルを当てはめます。切片と Acceleration についての Model_Year でグループ化された変量効果、および Weight についての Origin でグループ化された独立した変量効果を含めます。このモデルの式は次のようになります。

$$z_{imk}={\beta }_0+{\beta }_1{x}_i+{\beta }_2{\text{y}}_i+b_{10m}+b_{11m}{x}_i+b_{2k}{y}_i+{ϵ}_{imk},$$

ここで、$$x$$、$\mathit{y}$、$\mathit{z}$ は、それぞれ Acceleration、Weight、MPG を表します。添字 $\mathit{i}$ は MPG の観測値の行に対応し、添字 $$m=1,2,...,13$$ と $$k=1,2,...,8$$ は Model_Year と Origin の水準にそれぞれ対応します。変量効果係数と観測誤差の事前分布は次のとおりです。

$$b_{1m}=\left(\begin{array}{c}{b}_{10m}\\ b_{11m}\end{array}\right)\sim N(0,\left(\begin{array}{cc}{\sigma }_{10}^2& {\sigma }_{10,11}\\ {\sigma }_{10,11}& {\sigma }_{11}^2\end{array}\right)),$$

$$b_{2k}\sim N(0,{\sigma }_2^2),$$

$${ϵ}_{imk}\sim N(0,{\sigma }^2).$$

係数ベクトル $$b_{1m}$$ は、Model_Year の水準 $$m$$ における変量効果を表します。

• $$b_{10m}$$ は、水準 $$m$$ におけるランダム切片です。

• $$b_{11m}$$ は、水準 $$m$$ における Acceleration の変量効果係数です。

同様に、係数 $$b_{2k}$$ は、Origin の水準 $$k$$ における Weight の変量効果係数です。

$${\sigma }_{10}^2$$ は $$b_{10m}$$ の分散、$${\sigma }_{11}^2$$ は変量効果係数 $$b_{11m}$$ の分散、$${\sigma }_{10,11}$$ は $$b_{10m}$$ と $$b_{11m}$$ の共分散です。$${\sigma }_2^2$$ は $$b_{2k}$$ の変量効果係数の分散、$${\sigma }^2$$ は残差分散です。

線形混合効果モデルを当てはめるための計画行列を作成します。

F = [ones(406,1) Acceleration Weight];
R = {[ones(406,1) Acceleration],Weight};
Model_Year = nominal(Model_Year);
Origin = nominal(Origin);
G = {Model_Year,Origin};

F を固定効果、MPG を応答、R を変量効果、G をグループ化変数として使用してモデルを当てはめます。

lme = fitlmematrix(F,MPG,R,G,FixedEffectPredictors=....
["Intercept","Acceleration","Weight"],RandomEffectPredictors=...
{["Intercept","Acceleration"],"Weight"},RandomEffectGroups=["Model_Year","Origin"]);

変量効果についての共分散パラメーターの推定値と 99% 信頼区間を計算します。残差分散の平均二乗誤差を表示します。

[psi,mse,stats] = covarianceParameters(lme,Alpha=0.01);
mse

mse = 
9.0753

残差分散 mse は 9.0753 です。psi と stats は、共分散パラメーターの推定値とそれに関連する統計が格納された cell 配列です。

psi の最初の cell を表示します。

psi{1}

ans = 2×2

    8.2649   -0.8698
   -0.8698    0.1157

psi の最初の cell には、相関がある変量効果係数についての共分散パラメーターの推定値が格納されています。切片に対応する分散の推定値は 8.2649、Acceleration に対応する分散の推定値は 0.1157 です。切片と Acceleration に対応する共分散の推定値は -0.8698 です。

psi の 2 番目の cell を表示します。

psi{2}

ans = 
6.6770e-08

psi の 2 番目の cell には、Weight に対応する変量効果係数についての分散の推定値が格納されています。

stats の最初の cell を表示します。

stats{1}

ans = 
    Covariance Type: FullCholesky

    Group         Name1                   Name2                   Type            Estimate    Lower       Upper    
    Model_Year    {'Intercept'   }        {'Intercept'   }        {'std' }          2.8749     0.76378       10.821
    Model_Year    {'Acceleration'}        {'Intercept'   }        {'corr'}        -0.88949    -0.99322    0.0026856
    Model_Year    {'Acceleration'}        {'Acceleration'}        {'std' }         0.34015     0.16213      0.71364

切片と Acceleration に対応する変量効果係数についての標準偏差の推定値と 99% 信頼限界が出力に表示されます。また、グループ化変数 Model_Year の名前も出力に表示されています。標準偏差の推定値は psi の最初の cell の対角要素の平方根であることに注意してください。

stats の 2 番目の cell を表示します。

stats{2}

ans = 
    Covariance Type: FullCholesky

    Group     Name1             Name2             Type           Estimate     Lower         Upper    
    Origin    {'Weight'}        {'Weight'}        {'std'}        0.0002584    6.5446e-05    0.0010202

stats の 2 番目の cell には、Weight に対応する変量効果係数についての標準偏差の推定値と 99% 信頼限界が格納されています。

stats の 3 番目の cell を表示します。

stats{3}

ans = 
    Group    Name               Estimate    Lower     Upper 
    Error    {'Res Std'}        3.0125      2.7395    3.3127

stats の 3 番目の cell には、残差標準偏差の推定値と対応する 99% 信頼限界が格納されています。残差標準偏差の推定値は mse の平方根であることに注意してください。

入力引数

lme — 線形混合効果モデル

LinearMixedModel オブジェクト

線形混合効果モデル。fitlme または fitlmematrix を使用して構築した LinearMixedModel オブジェクトとして指定します。

alpha — 有意水準

0.05 (既定値) | 0 ～ 1 の範囲のスカラー値

有意水準。0 ～ 1 の範囲のスカラー値として指定します。値が α の場合、信頼水準は 100 × (1 – α)% です。

たとえば、99% の信頼区間の場合は、次のように信頼水準を指定できます。

例: Alpha=0.01

データ型: single | double

出力引数

psi — 共分散パラメーターの推定

cell 配列

変動効果の前の共分散をパラメーター表現する共分散パラメーターの推定。長さ R の cell 配列として返されます。psi{r} には、グループ化変数 g_r (r = 1, 2, ..., R) に関連付けられている変動効果の共分散行列が含まれます。グループ化変数の順序は、モデルを当てはめるときに入力する順序と同じです。

mse — 残差分散の推定

スカラー値

残差分散の推定。スカラー値として返します。

stats — 共分散パラメーターの推定および関連する統計

cell 配列

共分散パラメーターの推定および関連する統計。以下の列で構成されるデータセット配列を含む長さ (R + 1) の cell 配列として返します。

Group	グループ化変数名
Name1	最初の予測子変数の名前
Name2	2 番目の予測子変数の名前
Type	std (標準偏差)。Name1 と Name2 が同じ場合 corr (相関)。Name1 と Name2 が異なる場合
Estimate	Name1 と Name2 が同じ場合は、予測子 Name1 または Name2 と関連付けられている変量効果の標準偏差 Name1 と Name2 が異なる場合は、予測子 Name1 および Name2 と関連付けられている変量効果の間の相関
Lower	共分散パラメーターの 95% 信頼区間の下限
Upper	共分散パラメーターの 95% 信頼区間の上限

stats{r} は、r 番目のグループ化変数 (r = 1、2、...、R) の共分散パラメーターに対する統計を含むデータセット配列です。stats{R+1} には、残差標準偏差の統計値が格納されます。残差誤差のデータセット配列には、Group、Name、Estimate、Lower、Upper の各フィールドがあります。

バージョン履歴

R2013b で導入

参考

LinearMixedModel | compare | fixedEffects | randomEffects