coefTest

線形混合効果モデルの固定効果と変量効果についての仮説検定

構文

pVal = coefTest(lme)

pVal = coefTest(lme,H)

pVal = coefTest(lme,H,C)

pVal = coefTest(lme,H,C,Name,Value)

[pVal,F,DF1,DF2] = coefTest(___)

説明

pVal = coefTest(lme) は、切片を除くすべての固定効果係数が 0 であるという F 検定の p 値を返します。

pVal = coefTest(lme,H) は、対比行列 H を使用して、線形混合効果モデル lme の固定効果係数に対する F 検定の p 値を返します。これにより、帰無仮説 H₀: Hβ (β は固定効果ベクトル) が検定されます。

pVal = coefTest(lme,H,C) は、対比行列 H を使用して、線形混合効果モデル lme の固定効果係数に対する F 検定の p 値を返します。これにより、帰無仮説 H₀: Hβ = C (β は固定効果ベクトル) が検定されます。

pVal = coefTest(lme,H,C,Name,Value) は、1 つ以上の名前と値のペアの引数で指定された追加オプションを使用して、線形混合効果モデル lme の固定効果および/または変量効果係数に対する F 検定の p 値を返します。たとえば、'REContrast',K は帰無仮説 H₀: Hβ + KB = C (β は固定効果ベクトル、B は変量効果ベクトル) を検定するよう coefTest に指示します。

[pVal,F,DF1,DF2] = coefTest(___) はさらに F 統計量 F、F の分子と分母の自由度 DF1 と DF2 を返します。

例

カテゴリカル データの固定効果係数の検定

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

load('shift.mat')

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

Shift および Operator はノミナル変数です。

shift.Shift = nominal(shift.Shift);
shift.Operator = nominal(shift.Operator);

シフトの時間によってパフォーマンスに有意差があるかどうかを評価するために、作業者別のランダムな切片をもつ線形混合効果モデルを当てはめます。

lme = fitlme(shift,'QCDev ~ Shift + (1|Operator)')

lme = 
Linear mixed-effects model fit by ML

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
    59.012    62.552    -24.506          49.012  

Fixed effects coefficients (95% CIs):
    Name                     Estimate    SE         tStat       DF    pValue       Lower      Upper  
    {'(Intercept)'  }         3.1196     0.88681      3.5178    12    0.0042407     1.1874     5.0518
    {'Shift_Morning'}        -0.3868     0.48344    -0.80009    12      0.43921    -1.4401    0.66653
    {'Shift_Night'  }         1.9856     0.48344      4.1072    12    0.0014535    0.93227     3.0389

Random effects covariance parameters (95% CIs):
Group: Operator (5 Levels)
    Name1                  Name2                  Type           Estimate    Lower      Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        1.8297      0.94915    3.5272

Group: Error
    Name               Estimate    Lower      Upper 
    {'Res Std'}        0.76439     0.49315    1.1848

切片を除くすべての固定効果係数が 0 であるか検定します。

pVal = coefTest(lme)

pVal = 7.5956e-04

小さい $$p$$ 値は、0 ではない固定効果係数があることを示しています。

対比行列を使用して Shift の項の有意性を検定します。

H = [0 1 0; 0 0 1];
pVal = coefTest(lme,H)

pVal = 7.5956e-04

anova 法を使用して Shift の項の有意性を検定します。

anova(lme)

ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'RESIDUAL'

    Term                   FStat     DF1    DF2    pValue    
    {'(Intercept)'}        12.375    1      12      0.0042407
    {'Shift'      }        13.864    2      12     0.00075956

Shift の $$p$$ 値 0.00075956 は、前の仮説検定の $$p$$ 値と同じです。

夕方と朝のシフトに違いがあるか検定します。

pVal = coefTest(lme,[0 1 -1])

pVal = 3.6147e-04

この小さい $$p$$ 値は、朝のシフトと夕方のシフトでは作業者のパフォーマンスが異なることを示しています。

固定効果係数の仮説検定

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

load('weight.mat')

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

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

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

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

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

Program と Week の交互作用の有意性について検定します。

H = [0 0 0 0 0 0 1 0 0; 
     0 0 0 0 0 0 0 1 0;
     0 0 0 0 0 0 0 0 1];
pVal = coefTest(lme,H)

pVal = 0.9775

大きい $$p$$ 値は、Program と Week の交互作用が統計的に有意ではないことを示しています。

次に Program を含むすべての係数が 0 であるかを検定します。

H = [0 0 1 0 0 0 0 0 0;
     0 0 0 1 0 0 0 0 0;
     0 0 0 0 1 0 0 0 0;
     0 0 0 0 0 0 1 0 0; 
     0 0 0 0 0 0 0 1 0;
     0 0 0 0 0 0 0 0 1];
C = [0;0;0;0;0;0];
pVal = coefTest(lme,H,C)

pVal = 0.0274

0.0274 という $$p$$ 値は、Program に関連する係数にゼロではないものがあることを示しています。

固定効果係数および変量効果係数の仮説検定

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

load flu

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

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

flu2 = stack(flu,2:10,'NewDataVarName','FluRate',...
    'IndVarName','Region');
flu2.Date = nominal(flu2.Date);

地域に対する固定効果と、Date で変化するランダム切片で、線形混合効果モデルを当てはめます。

lme = fitlme(flu2,'FluRate ~ 1 + Region + (1|Date)')

lme = 
Linear mixed-effects model fit by ML

Model information:
    Number of observations             468
    Fixed effects coefficients           9
    Random effects coefficients         52
    Covariance parameters                2

Formula:
    FluRate ~ 1 + Region + (1 | Date)

Model fit statistics:
    AIC       BIC       LogLikelihood    Deviance
    318.71    364.35    -148.36          296.71  

Fixed effects coefficients (95% CIs):
    Name                        Estimate    SE          tStat      DF     pValue        Lower        Upper    
    {'(Intercept)'     }          1.2233    0.096678     12.654    459     1.085e-31       1.0334       1.4133
    {'Region_MidAtl'   }        0.010192    0.052221    0.19518    459       0.84534    -0.092429      0.11281
    {'Region_ENCentral'}        0.051923    0.052221     0.9943    459        0.3206    -0.050698      0.15454
    {'Region_WNCentral'}         0.23687    0.052221     4.5359    459    7.3324e-06      0.13424      0.33949
    {'Region_SAtl'     }        0.075481    0.052221     1.4454    459       0.14902     -0.02714       0.1781
    {'Region_ESCentral'}         0.33917    0.052221      6.495    459    2.1623e-10      0.23655      0.44179
    {'Region_WSCentral'}           0.069    0.052221     1.3213    459       0.18705    -0.033621      0.17162
    {'Region_Mtn'      }        0.046673    0.052221    0.89377    459       0.37191    -0.055948      0.14929
    {'Region_Pac'      }        -0.16013    0.052221    -3.0665    459     0.0022936     -0.26276    -0.057514

Random effects covariance parameters (95% CIs):
Group: Date (52 Levels)
    Name1                  Name2                  Type           Estimate    Lower     Upper  
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0.6443      0.5297    0.78368

Group: Error
    Name               Estimate    Lower      Upper
    {'Res Std'}        0.26627     0.24878    0.285

10/9/2005 の週の変量効果の項が 0 であるという仮説を検定します。

[~,~,STATS] = randomEffects(lme); % Compute the random-effects statistics (STATS)
STATS.Level = nominal(STATS.Level);
K = zeros(length(STATS),1);
K(STATS.Level == '10/9/2005') = 1;
pVal = coefTest(lme,[0 0 0 0 0 0 0 0 0],0,'REContrast',K')

pVal = 0.1692

次にこのモデルをランダムな切片と勾配のモデルに再度当てはめます。

lme = fitlme(flu2,'FluRate ~ 1 + Region + (1 + Region|Date)');

10/9/2005 の週の WNCentral 地域の結合された係数が 0 であるという仮説を検定します。

[~,~,STATS] = randomEffects(lme); STATS.Level = nominal(STATS.Level);
K = zeros(length(STATS),1);
K(STATS.Level == '10/9/2005' & flu2.Region == 'WNCentral') = 1;
pVal = coefTest(lme,[0 0 0 1 0 0 0 0 0],0,'REContrast',K')

pVal = 1.4770e-12

また、$$F$$ 統計量を分子および分母の自由度と共に取得します。

[pVal,F,DF1,DF2] = coefTest(lme,[0 0 0 1 0 0 0 0 0],0,'REContrast',K')

pVal = 1.4770e-12

F = 52.9730

DF1 = 1

DF2 = 459

サタースウェイトの近似法を分母の自由度に使用してこの検定を繰り返します。

[pVal,F,DF1,DF2] = coefTest(lme,[0 0 0 1 0 0 0 0 0],0,'REContrast',K',...
			'DFMethod','satterthwaite')

pVal = NaN

F = 52.9730

DF1 = 1

DF2 = 0

入力引数

lme — 線形混合効果モデル

LinearMixedModel オブジェクト

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

H — 固定効果の対比

m 行 p 列の行列

固定効果の対比。m 行 p 列の行列として指定します。ここで p は lme 内の固定効果係数の数です。H の各行が 1 つの対比を表します。H の各列 (左から右) は、fixedEffects メソッドによって返される p 行 1 列の固定効果ベクトル beta (上から下) の各行に対応します。

例: pVal = coefTest(lme,H)

データ型: single | double

C — 仮定された値

m 行 1 列のベクトル

帰無仮説 H*beta = C の検定に使用する仮定の値。m 行 1 列の行列を指定します。ここで beta は fixedEffects メソッドによって返される固定効果の推定値のベクトルです。

データ型: single | double

名前と値の引数

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

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

例: pVal = coefTest(lme,H,C,'DFMethod','satterthwaite')

DFMethod — 分母の自由度の近似の計算方法

'residual' (既定値) | 'satterthwaite' | 'none'

F 検定における分母の自由度の近似の計算方法です。'DFMethod' と次のいずれかの値で構成されるコンマ区切りのペアで指定します。

'residual'	既定の設定。自由度は定数で n – p に等しいと仮定されます。ここで n は観測値の数、p は固定効果の数です。
'satterthwaite'	サタースウェイトの近似法。
'none'	すべての自由度は無限大に設定されます。

たとえば、次のようにサタースウェイトの近似法を指定できます。

例: 'DFMethod','satterthwaite'

REContrast — 変量効果の対比

m 行 q 列の行列

変量効果の対比。'REContrast' と m 行 q 列の行列 K から構成されるコンマ区切りのペアとして指定します。q は lme 内の変量効果パラメーターの個数です。K の各列 (左から右) は、randomEffects メソッドによって返される変量効果の最良線形不偏予測量ベクトル B (上から下) の各行に一致します。

データ型: single | double

出力引数

pVal — p 値

スカラー値

線形混合効果モデル lme の固定効果および/または変量効果係数に対する F 検定の p 値です。スカラー値として返されます。

F — F 統計量

スカラー値

F 統計量です。スカラー値として返されます。

DF1 — F の分子の自由度

スカラー値

F の分子の自由度です。スカラー値として返されます。

• 帰無仮説 H₀: Hβ = 0 または H₀: Hβ = C を検定すると、DF1 は H にある線形的に独立した行の数と等しくなります。

• 帰無仮説 H₀: Hβ + KB = C を検定する場合、DF1 は [H,K] にある線形的に独立した行の数と等しくなります。

DF2 — F の分母の自由度

スカラー値

F の分母の自由度です。スカラー値として返されます。DF2 の値は DFMethod で選択したオプションに依存します。

バージョン履歴

R2013b で導入

参考

LinearMixedModel | anova | coefCI