Main Content

anova

線形混合効果モデルの分散分析

ページ内をすべて折りたたむ

構文

stats = anova(lme)
stats = anova(lme,Name,Value)

説明

例

stats = anova(lme) はデータセット配列 stats を返します。これには、線形混合効果モデル lme 内における固定効果の項ごとの F 検定結果が含まれています。

例

stats = anova(lme,Name,Value) も、1 つ以上の Name,Value の引数ペアで指定された追加オプションを使用して、データセット配列 stats を返します。

例

すべて折りたたむ

ライブ スクリプトを開く

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

load('shift.mat')

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

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

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

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

'effects' の対比では、係数の合計が 0 になることを指定します。fitlme は、2 つの対比コード化された変数 $X$1 および $X$2 を固定効果計画行列内に作成します。ここで、次のようになります。

Shift_Evening={0,ifMorning1,ifEvening-1,ifNight

および

Shift_Morning={1,ifMorning0,ifEvening-1,ifNight.

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

MorningShift:QCDevim=β0+β2Shift_Morningi+b0m+εim,m=1,2,...,5,EveningShift:QCDevim=β0+β1Shift_Eveningi+b0m+εim,NightShift:QCDevim=β0-β1Shift_Eveningi-β2Shift_Morningi+b0m+εim,

ここで、b ~ N(0, σb2 ) および ϵ ~ N(0, σ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       Lower      Upper   
    {'(Intercept)'  }          3.6525    0.94109     3.8812    12    0.0021832     1.6021       5.703
    {'Shift_Evening'}        -0.53293    0.31206    -1.7078    12      0.11339    -1.2129     0.14699
    {'Shift_Morning'}        -0.91973    0.31206    -2.9473    12     0.012206    -1.5997    -0.23981

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

すべての固定効果係数が 0 であるかどうか判定するため、F 検定を実行します。 

anova(lme)
ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'RESIDUAL'

    Term                   FStat     DF1    DF2    pValue   
    {'(Intercept)'}        15.063    1      12     0.0021832
    {'Shift'      }        11.091    2      12     0.0018721

定数項の p 値 0.0021832 は、lme の係数表の値と同じです。Shift の p 値 0.0018721 は、Shift を表す両方の係数の結合有意性の尺度となります。

ライブ スクリプトを開く

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

load('fertilizer.mat')

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

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

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

線形混合効果モデルを当てはめます。Fertilizer および Tomato は固定効果変数であり、平均収穫量はブロック (土壌の種類) とブロック内のプロット (土壌の種類の中のトマトの種類) によって独立して変化します。データをタイプ III 二乗和に近似するときに 'effects' 対比を使用します。 

lme = fitlme(ds,'Yield ~ Fertilizer * Tomato + (1|Soil) + (1|Soil:Tomato)',...
'DummyVarCoding','effects')
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    pValue        Lower      Upper  
    {'(Intercept)'                 }           104.6    3.3008       31.69    40    5.9086e-30     97.929     111.27
    {'Tomato_Cherry'               }             1.4    5.9353     0.23588    40       0.81473    -10.596     13.396
    {'Tomato_Grape'                }         -7.7667    5.9353     -1.3085    40       0.19816    -19.762     4.2291
    {'Tomato_Heirloom'             }         -11.183    5.9353     -1.8842    40      0.066821    -23.179    0.81242
    {'Tomato_Plum'                 }          30.233    5.9353      5.0938    40     8.777e-06     18.238     42.229
    {'Fertilizer_1'                }         -28.267    2.3475     -12.041    40    7.0265e-15    -33.011    -23.522
    {'Fertilizer_2'                }         -1.9333    2.3475    -0.82356    40       0.41507    -6.6779     2.8112
    {'Fertilizer_3'                }          10.733    2.3475      4.5722    40     4.577e-05     5.9888     15.478
    {'Tomato_Cherry:Fertilizer_1'  }        -0.73333    4.6951    -0.15619    40       0.87667    -10.222     8.7558
    {'Tomato_Grape:Fertilizer_1'   }         -7.5667    4.6951     -1.6116    40       0.11491    -17.056     1.9224
    {'Tomato_Heirloom:Fertilizer_1'}          5.1833    4.6951       1.104    40       0.27619    -4.3058     14.672
    {'Tomato_Plum:Fertilizer_1'    }          2.7667    4.6951     0.58927    40       0.55899    -6.7224     12.256
    {'Tomato_Cherry:Fertilizer_2'  }             7.6    4.6951      1.6187    40       0.11337    -1.8891     17.089
    {'Tomato_Grape:Fertilizer_2'   }            -1.9    4.6951    -0.40468    40       0.68787    -11.389     7.5891
    {'Tomato_Heirloom:Fertilizer_2'}          5.5167    4.6951       1.175    40       0.24695    -3.9724     15.006
    {'Tomato_Plum:Fertilizer_2'    }            -3.9    4.6951    -0.83066    40        0.4111    -13.389     5.5891
    {'Tomato_Cherry:Fertilizer_3'  }         -6.0667    4.6951     -1.2921    40       0.20373    -15.556     3.4224
    {'Tomato_Grape:Fertilizer_3'   }          3.7667    4.6951     0.80226    40       0.42714    -5.7224     13.256
    {'Tomato_Heirloom:Fertilizer_3'}          3.1833    4.6951     0.67802    40       0.50167    -6.3058     12.672
    {'Tomato_Plum:Fertilizer_3'    }             1.1    4.6951     0.23429    40       0.81596    -8.3891     10.589

Random effects covariance parameters (95% CIs):
Group: Soil (3 Levels)
    Name1                  Name2                  Type           Estimate    Lower       Upper 
    {'(Intercept)'}        {'(Intercept)'}        {'std'}        2.5028      0.027711    226.05

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

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

固定効果の検定のために分散分析を実施します。 

anova(lme)
ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'RESIDUAL'

    Term                         FStat     DF1    DF2    pValue    
    {'(Intercept)'      }        1004.2     1     40     5.9086e-30
    {'Tomato'           }        7.1663     4     40     0.00018935
    {'Fertilizer'       }        58.833     3     40     1.0024e-14
    {'Tomato:Fertilizer'}        1.4182    12     40        0.19804

定数項の p 値 5.9086e-30 は、lme の係数表の値と同じです。Tomato、Fertilizer および Tomato:Fertilizer の p 値 0.00018935、1.0024e-14 および 0.19804 はそれぞれ、すべてのトマト係数、肥料係数、およびトマトと肥料の交互作用を表す係数の結合有意性を表します。p 値 0.19804 は、トマトと肥料の間の交互作用が有意ではないことを示しています。

ライブ スクリプトを開く

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

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);

'effects' 対比を使用してモデルを当てはめます。 

lme = fitlme(tbl,'y ~ InitialWeight + Program*Week + (Week|Subject)',...
		'DummyVarCoding','effects')
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.77122      0.24309      3.1725    111     0.0019549       0.28951       1.2529
    {'InitialWeight' }        0.0031879    0.0013814      2.3078    111      0.022863    0.00045067    0.0059252
    {'Program_A'     }         -0.11017     0.080377     -1.3707    111       0.17323      -0.26945       0.0491
    {'Program_B'     }          0.25061      0.08045      3.1151    111     0.0023402      0.091195      0.41003
    {'Program_C'     }         -0.14344     0.080475     -1.7824    111      0.077424       -0.3029     0.016031
    {'Week'          }          0.19881     0.033727      5.8946    111    4.1099e-08       0.13198      0.26564
    {'Program_A:Week'}        -0.025607     0.058417    -0.43835    111       0.66198      -0.14136     0.090149
    {'Program_B:Week'}         0.013164     0.058417     0.22535    111       0.82212      -0.10259      0.12892
    {'Program_C:Week'}        0.0049357     0.058417    0.084492    111       0.93282      -0.11082      0.12069

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.21076    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

p 値 0.022863 および 4.1099e-08 は、被験者の初期体重と時間因子が体重の減少量に対して有意な影響を与えることを示しています。プログラム B の被験者の体重の減少は、プログラム A の被験者の体重の減少に対して有意差があります。変量効果の共分散パラメーターの下限または上限に 0 は含まれないため、この差は有意です。

すべての固定効果係数が 0 である F 検定を実施します。 

anova(lme)
ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'RESIDUAL'

    Term                     FStat       DF1    DF2    pValue    
    {'(Intercept)'  }          10.065    1      111     0.0019549
    {'InitialWeight'}           5.326    1      111      0.022863
    {'Program'      }          3.6798    3      111      0.014286
    {'Week'         }          34.747    1      111    4.1099e-08
    {'Program:Week' }        0.066648    3      111       0.97748

定数項、初期体重および週に関する p 値は、前の lme の出力で表示された係数表の値と同じです。Program の p 値 0.014286 は、すべてのプログラムの係数の結合有意性を表します。同様に、プログラムと週の交互作用 (Program:Week) の p 値は、この交互作用を表すすべての係数の結合有意性の尺度となります。

次にサタースウェイト法を使用して自由度を計算します。 

anova(lme,'DFMethod','satterthwaite')
ans = 
    ANOVA MARGINAL TESTS: DFMETHOD = 'SATTERTHWAITE'

    Term                     FStat       DF1    DF2       pValue    
    {'(Intercept)'  }          10.065    1      20.445      0.004695
    {'InitialWeight'}           5.326    1          20      0.031827
    {'Program'      }          3.6798    3       19.14      0.030233
    {'Week'         }          34.747    1          20    9.1347e-06
    {'Program:Week' }        0.066648    3          20       0.97697

サタースウェイト法では、分母の自由度が小さくなり、p 値が少しだけ大きくなります。

入力引数

すべて折りたたむ

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

名前と値の引数

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

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

例: stats = anova(lme,'DFMethod','satterthwaite')

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

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

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

例: 'DFMethod','satterthwaite'

出力引数

すべて折りたたむ

固定効果項の F 検定の結果です。次の列をもつデータセット配列として返されます。

Term固定効果項の名前
Fstat項の F 統計量
DF1F 統計量の分子の自由度
DF2F 統計量の分母の自由度
pValue項の検定の p 値

固定効果の項ごとに 1 つの行があります。それぞれの項は連続変数、グループ化変数、2 つ以上の連続変数またはグループ化変数間の交互作用のいずれかです。anova は固定効果の項ごとに、固定効果の項を表すすべての係数が 0 であるかを判定する F 検定 (周辺検定) を行います。タイプ III 仮説の検定を行うには、線形混合効果モデルを当てはめながら、'effects' の対比を使用しなければなりません。

ヒント

  • anova は固定効果の項ごとに、固定効果の項を表すすべての係数が 0 であることについての F 検定 (周辺検定) を行います。タイプ III 仮説の検定を行うには、線形混合効果モデルを当てはめながら、名前と値のペアの引数 'DummyVarCoding' を 'effects' の対比に設定しなければなりません。

バージョン履歴

R2013b で導入

参考

| |