Main Content

covarianceParameters

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

説明

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

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

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

[psi,mse,stats] = covarianceParameters(lme,Name,Value) は、1 つ以上の Name,Value のペアの引数によって指定された追加オプションで stats の共分散パラメーターおよび関連する統計を返します。

たとえば、共分散パラメーターの信頼限界に対する信頼水準を指定できます。

すべて折りたたむ

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

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

yijk=β0+j=25β2jI[T]ij+b0jk(S*T)jk+ϵijk,

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

変量効果と観測誤差の事前分布は、b0kN(0,σS2)b0jkN(0,σS*T2) および ϵijkN(0,σ2) です。

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

変量効果項の共分散パラメーターの推定値 (σS2σS*T2 の推定値) を計算します。

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

残差分散 (σ2) を計算します。

[~,mse] = covarianceParameters(lme)
mse = 151.9007

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

load('weight.mat');

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

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

ds = dataset(InitialWeight,Program,Subject,Week,y);
ds.Subject = nominal(ds.Subject);
ds.Program = nominal(ds.Program);

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

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

yim=β0+β1IWi+β2Weeki+β3I[PB]I+β4I[PC]i+β5I[PD]i+b0m+b1mWeekim+ϵim

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

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

(b0mb1m)N(0,(σ02σ0,1σ0,1σ12))

および

ϵimN(0,σ2).

lme = fitlme(ds,'y ~ InitialWeight + Program + (Week|Subject)');

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

[psi,mse,stats] = covarianceParameters(lme)
psi = 1x1 cell array
    {2x2 double}

mse = 0.0105
stats=2×1 cell array
    {3x7 classreg.regr.lmeutils.titleddataset}
    {1x5 classreg.regr.lmeutils.titleddataset}

mse は推定された残差分散です。これは、σ2 の推定値です。

変量効果項 (σ02σ12 および σ0,1) に対する共分散パラメーターの推定値を確認するには、psi のインデックスを指定します。

psi{1}
ans = 2×2

    0.0572    0.0490
    0.0490    0.0624

切片に対する変量効果項の分散の推定値 σ02 は 0.0572 です。週に対する変量効果項の分散の推定値 σ12 は 0.0624 です。切片と週に対する変量効果項の共分散の推定値 σ0,1 は 0.0490 です。

stats は 2 行 1 列の cell 配列です。stats の最初のセルには、変量効果の標準偏差の信頼区間と、切片と週に関する変量効果間の相関が含まれます。これらを表示するには、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)、変量効果変数 (Name1Name2)、共分散パラメーターの種類 (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 番目のセルには、残差標準偏差に関する同様の統計が含まれます。2 番目のセルの内容を表示します。

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

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

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

load carbig

ガロンあたりの走行マイル数 (MPG) の線形混合効果モデルを当てはめます。加速度および重量は固定効果で、モデル年度によってグループ化される切片と加速度は相関された変量効果の可能性があり、重量は独立変量効果であり、自動車の原点でグループ化されます。このモデルは以下に対応します。

MPGimk=β0+β1Acci+β2Weighti+b10m+b11mAcci+b21kWeighti+ϵimk

ここで、m=1,2,...,13 は変数 Model_Year の水準を表し、k=1,2,...,8 は変数 Origin の水準を表します。MPGimk は、i 番目の観測値に対応する m 番目のモデル年および k 番目の生産国のガロンあたりの走行マイル数です。変量効果項と観測値の誤差の事前分布は次のとおりです。

b1m=(b10mb11m)N(0,(σ102σ10,11σ10,11σ112)),

b2kN(0,σ22),

ϵimkN(0,σ2).

ここで、変量効果項 b1m は 1 番目のグループ化変数の水準 m における 1 番目の変量効果を表します。変量効果項 b10m は、1 番目のグループ化変数の m 番目の水準 (m) における切片 (0) に対する 1 番目の変量効果項 (1) に対応します。同様に、b11m は 1 番目の変量効果項 (1) における 1 番目の予測子 (1) の水準 m です。

同様に、b2k は 2 番目のグループ化変数の水準 k における 2 番目の変量効果項を表します。

σ102 は切片に対する変量効果項の分散、σ112 は予測子 Acceleration に対する変量効果項の分散、σ10,11 は切片と予測子 Acceleration に対する変量効果項の共分散です。σ22 は 2 番目の変量効果項の分散、σ2 は残差分散です。

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

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

計画行列を使用してモデルを当てはめます。

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

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

[psi,mse,stats] = covarianceParameters(lme)
psi=2×1 cell array
    {2x2 double  }
    {[6.6770e-08]}

mse = 9.0753
stats=3×1 cell array
    {3x7 classreg.regr.lmeutils.titleddataset}
    {1x7 classreg.regr.lmeutils.titleddataset}
    {1x5 classreg.regr.lmeutils.titleddataset}

残差分散 mse は 9.0755 です。psi は 2 行 1 列の cell 配列であり、stats は 3 行 1 列の cell 配列です。内容を確認するには、これらの cell 配列のインデックスを指定しなければなりません。

まず、psi の最初のセルのインデックスを指定します。

psi{1}
ans = 2×2

    8.2649   -0.8698
   -0.8698    0.1157

psi の 1 番目のセルには、切片に対する相関性がある変量効果の共分散パラメーター σ102 として 8.5160、加速度に対する σ112 として 0.1087 が格納されています。切片と加速度に対する変量効果項の共分散の推定値 σ10,11 は -0.8387 です。

次に、psi の 2 番目のセルのインデックスを指定します。

psi{2}
ans = 6.6770e-08

psi の 2 番目のセルには、重量に対する変量効果項の分散の推定値 σ22 が格納されています。

stats の最初のセルのインデックスを指定します。

stats{1}
ans = 
    COVARIANCE TYPE: FULLCHOLESKY

    Group         Name1                   Name2                   Type            Estimate    Lower       Upper  
    Model_Year    {'Intercept'   }        {'Intercept'   }        {'std' }          2.8749      1.0486      7.882
    Model_Year    {'Acceleration'}        {'Intercept'   }        {'corr'}        -0.88949    -0.98666    -0.3251
    Model_Year    {'Acceleration'}        {'Acceleration'}        {'std' }         0.34015     0.19355    0.59777

この表では、切片と加速度に関する変量効果の項の標準偏差が推定して示されています。標準偏差の推定は、psi の最初のセルの対角要素の平方根であることに注意してください。具体的には、2.9182 = sqrt(8.5160) および 0.32968 = sqrt(0.1087) です。相関は、切片と加速度の共分散および切片と加速度の標準偏差の関数です。切片と加速度の共分散は、psi の 1 番目のセルの非対角値 -0.8387 です。したがって、相関は -.8387/(0.32968*2.92182) = -0.87 です。

切片と加速度のグループ化変数は Model_Year です。

stats の 2 番目のセルのインデックスを指定します。

stats{2}
ans = 
    COVARIANCE TYPE: FULLCHOLESKY

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

stats の 2 番目のセルには、標準偏差の推定と、Weight の変量効果の項の標準偏差の 95% 信頼限界が含まれます。グループ化変数は Origin です。

stats の 3 番目のセルのインデックスを指定します。

stats{3}
ans = 
    Group    Name               Estimate    Lower     Upper 
    Error    {'Res Std'}        3.0125      2.8025    3.2383

stats の 3 番目のセルには、残差標準偏差の推定および 95% の信頼限界が含まれます。残差標準偏差の推定は、mse の平方根 sqrt(9.0755) = 3.0126 です。

共分散パラメーターに対する 99% 信頼区間を作成します。

[~,~,stats] = covarianceParameters(lme,'Alpha',0.01);
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

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{3}
ans = 
    Group    Name               Estimate    Lower     Upper 
    Error    {'Res Std'}        3.0125      2.7395    3.3127

入力引数

すべて折りたたむ

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

名前と値の引数

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

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

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

有意水準。'Alpha' と 0 ~ 1 の範囲にあるスカラー値から構成されるコンマ区切りのペアとして指定します。値が α の場合、信頼水準は 100 × (1 – α)% です。

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

例: 'Alpha',0.01

データ型: single | double

出力引数

すべて折りたたむ

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

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

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

Groupグループ化変数名
Name1最初の予測子変数の名前
Name22 番目の予測子変数の名前
Type

std (標準偏差)。Name1Name2 が同じ場合

corr (相関)。Name1Name2 が異なる場合

Estimate

Name1Name2 が同じ場合は、予測子 Name1 または Name2 と関連付けられている変量効果の標準偏差

Name1Name2 が異なる場合は、予測子 Name1 および Name2 と関連付けられている変量効果の間の相関

Lower共分散パラメーターの 95% 信頼区間の下限
Upper共分散パラメーターの 95% 信頼区間の上限

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

バージョン履歴

R2013b で導入