randomEffects

クラス: GeneralizedLinearMixedModel

変量効果と関連する統計の推定

構文

B = randomEffects(glme)

[B,BNames] = randomEffects(glme)

[B,BNames,stats] = randomEffects(glme)

[B,BNames,stats] = randomEffects(glme,Name,Value)

説明

B = randomEffects(glme) は、推定された共分散パラメーターと観測された応答を条件とする一般化線形混合効果モデル glme における変量効果の EBP (経験的ベイズ予測子) の推定値を返します。

[B,BNames] = randomEffects(glme) はまた、係数 BNames の名前を返します。それぞれの名前は、B の係数に相当します。

[B,BNames,stats] = randomEffects(glme) はまた、glme の変量効果の推定 EBP の関連する統計 stats を返します。

[B,BNames,stats] = randomEffects(glme,Name,Value) は、1 つ以上の Name,Value のペアの引数により指定される追加オプションを使用し、上記から任意の出力引数を返します。たとえば、信頼区間レベルまたは近似の自由度を計算するメソッドを指定できます。

入力引数

glme — 一般化線形混合効果モデル

GeneralizedLinearMixedModel オブジェクト

一般化線形混合効果モデル。GeneralizedLinearMixedModel オブジェクトとして指定します。このオブジェクトのプロパティとメソッドについては、GeneralizedLinearMixedModel を参照してください。

名前と値の引数

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

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

Alpha — 有意水準

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

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

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

例: 'Alpha',0.01

データ型: single | double

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

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

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

値	説明
'residual'	自由度の値は n - p に等しい定数であると仮定されます。n は観測値の個数、p は固定効果の個数です。
'none'	自由度は無限大に設定されます。

例: 'DFMethod','none'

出力引数

B — 変量効果の経験的ベイズ予測子を推定します。

列ベクトル

一般化線形混合効果モデル glme の変量効果の EBP (経験的ベイズ予測子)。列ベクトルとして返します。B の EBP は、推定された共分散パラメーターと観測された応答が与えられた、変量効果の経験事後分布の最頻値により近似されます。

glme には R グループ化変数 g₁、g₂、...、g_R があり、それぞれ水準が m₁、m₂、...、m_R であると仮定します。また、q₁、q₂、...、q_R は、それぞれ g₁, g₂、...、g_R に関連付けられている変量効果の長さであるとします。B は、次の長さの列ベクトルです。q₁*m₁ + q₂*m₂ + ... + q_R*m_R。

randomEffects によって作成される B は、各グループ化変数の各水準に対応する変量効果ベクトルの経験的ベイズ予測子を、[g1level1; g1level2; ...; g1levelm1; g2level1; g2level2; ...; g2levelm2; ...; gRlevel1; gRlevel2; ...; gRlevelmR]' として連結します。

BNames — 変量効果係数の名前

table

B の変量効果係数の名前。table として返します。

stats — 経験的ベイズ予測子と関連する統計を推定します。

table

一般化線形混合効果モデル glme の変量効果の EBP (経験的ベイズ予測) と関連する統計。テーブルとして返されます。stats は各変量効果に対して 1 行、以下の各統計に対して 1 列をもちます。

列名	説明
Group	変量効果に関連付けられているグループ化変数
Level	変量効果に対応するグループ化変数内の水準
Name	変量効果係数の名前
Estimate	変量効果の EBP (経験的ベイズ予測)
SEPred	共分散パラメーターおよび応答が与えられた CMSEP (予測の条件付き平均二乗誤差) の平方根
tStat	変量効果係数が 0 に等しい検定の t 統計
DF	t 統計の推定自由度
pValue	t 統計の p 値
Lower	変量効果係数の 95% 信頼区間の下限
Upper	変量効果係数の 95% の信頼区間の上限

randomEffects は推定された共分散パラメーターおよび観測された応答での CMSEP (予測の条件付き平均二乗誤差) メソッドを使用して、信頼区間を計算します。信頼区間の別の解釈は、信頼区間が共分散パラメーターと観測された応答を条件とする近似のベイズの信頼できる区間であるということです。

fitglme と疑似尤度の近似メソッド ('MPL' または 'REMPL') のいずれかを使用して GLME モデルを当てはめる場合、randomEffects は、疑似尤度の最後の反復から当てはめた線形混合効果モデルに基づき、信頼区間と関係する統計を計算します。

例

推定変量効果の計算とプロット

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

load mfr

このシミュレーションされたデータは、世界中で 50 の工場を操業している製造企業から取得しており、各工場が完成品の生産のためにバッチ処理を実行しています。同社は各バッチの欠陥数を減少させるために新たな製造プロセスを開発しました。新しいプロセスの効果をテストするため、同社は実験に参加させる 20 工場を無作為に選びました。10 工場では新プロセスを実施しますが、残りの 10 工場では旧プロセスの実行を続けます。各 20 工場で、同社は 5 つのバッチ (合計 100 バッチ) を実行し以下のデータを記録しました。

• 新しいプロセスがバッチに使用されたかどうかを示すフラグ (newprocess)

• 各バッチの処理時間。時間単位 (time)

• バッチの温度。摂氏 (temp)

• バッチで使用する化学薬品の供給業者 (A、B または C) を示すカテゴリカル変数 (supplier)

• バッチ内の欠陥数 (defects)

またデータに含まれる time_dev と temp_dev は、摂氏 20 度で 3 時間の標準プロセスから得られる時間と温度の絶対偏差をそれぞれ表します。

固定効果予測子として newprocess、time_dev、temp_dev および supplier を使用して一般化線形混合効果モデルを当てはめます。工場特有の変動に起因して品質に差がある可能性を考慮するために、factory 別にグループ化された切片の変量効果項を含めます。応答変数 defects はポアソン分布であり、このモデルの適切なリンク関数は対数です。係数の予測にラプラス近似メソッドを使用します。ダミー変数エンコードを 'effects' として指定すると、ダミー変数の係数の合計が 0 になります。

欠陥数はポアソン分布を使用してモデル化できます

$${\text{defects}}_{ij}\sim \text{Poisson}({\mu }_{ij})$$

これは一般化線形混合効果モデルに対応します

$$\log ({\mu }_{ij})={\beta }_0+{\beta }_1{\text{newprocess}}_{ij}+{\beta }_2{\text{time}\text{\_}\text{dev}}_{ij}+{\beta }_3{\text{temp}\text{\_}\text{dev}}_{ij}+{\beta }_4{\text{supplier}\text{\_}\text{C}}_{ij}+{\beta }_5{\text{supplier}\text{\_}\text{B}}_{ij}+b_i,$$

ここで

• $${\text{defects}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ で実行されたバッチで観測された欠陥数です。

• $${\mu }_{ij}$$ は、バッチ $$j$$ ($$j=1,2,...,5$$) 処理中の工場 $$i$$ ($$i=1,2,...,20$$) に対応する欠陥の平均数です。

• $${\text{newprocess}}_{ij}$$、$${\text{time}\text{\_}\text{dev}}_{ij}$$ および $${\text{temp}\text{\_}\text{dev}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ に対応する各変数の測定値です。たとえば $${\text{newprocess}}_{ij}$$ は、工場 $$i$$ で実行されたバッチ $$j$$ 処理中に新プロセスが使用されたかどうかを示します。

• $${\text{supplier}\text{\_}\text{C}}_{ij}$$ および $${\text{supplier}\text{\_}\text{B}}_{ij}$$ はエフェクト (ゼロサム) コーディングを使用するダミー変数であり、バッチ $$j$$ 処理中に工場 $$i$$ で実行されたバッチに対して、それぞれ会社 C または B が加工化学薬品を供給したかどうかを示します。

• $$b_i\sim N(0,{\sigma }_b^2)$$ は、工場特有の品質変動に相当する、各工場 $$i$$ の変量効果の切片です。

glme = fitglme(mfr,'defects ~ 1 + newprocess + time_dev + temp_dev + supplier + (1|factory)','Distribution','Poisson','Link','log','FitMethod','Laplace','DummyVarCoding','effects');

変量効果の EBP (経験的ベイズ予測) の名前と推定値を計算し、表示します。

[B,BNames] = randomEffects(glme)

B = 20×1

    0.2913
    0.1542
   -0.2633
   -0.4257
    0.5453
   -0.1069
    0.3040
   -0.1653
   -0.1458
   -0.0816
    0.0145
    0.1771
    0.2487
    0.2115
    0.2777
      ⋮

BNames=20×3 table
       Group       Level          Name      
    ___________    ______    _______________

    {'factory'}    {'1' }    {'(Intercept)'}
    {'factory'}    {'2' }    {'(Intercept)'}
    {'factory'}    {'3' }    {'(Intercept)'}
    {'factory'}    {'4' }    {'(Intercept)'}
    {'factory'}    {'5' }    {'(Intercept)'}
    {'factory'}    {'6' }    {'(Intercept)'}
    {'factory'}    {'7' }    {'(Intercept)'}
    {'factory'}    {'8' }    {'(Intercept)'}
    {'factory'}    {'9' }    {'(Intercept)'}
    {'factory'}    {'10'}    {'(Intercept)'}
    {'factory'}    {'11'}    {'(Intercept)'}
    {'factory'}    {'12'}    {'(Intercept)'}
    {'factory'}    {'13'}    {'(Intercept)'}
    {'factory'}    {'14'}    {'(Intercept)'}
    {'factory'}    {'15'}    {'(Intercept)'}
    {'factory'}    {'16'}    {'(Intercept)'}
      ⋮

B の各行には、Bnames の対応する行で名前を指定された変量効果係数の推定の EPB が含まれます。たとえば、B の行 3 の値 –0.2633 は、factory の水準 '3' の '(Intercept)' に推定される EPB です。

変量効果に対する 99% の信頼区間を計算します。

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

load mfr

このシミュレーションされたデータは、世界中で 50 の工場を操業している製造企業から取得しており、各工場が完成品の生産のためにバッチ処理を実行しています。同社は各バッチの欠陥数を減少させるために新たな製造プロセスを開発しました。新しいプロセスの効果をテストするため、同社は実験に参加させる 20 工場を無作為に選びました。10 工場では新プロセスを実施しますが、残りの 10 工場では旧プロセスの実行を続けます。各 20 工場で、同社は 5 つのバッチ (合計 100 バッチ) を実行し以下のデータを記録しました。

• 新しいプロセスがバッチに使用されたかどうかを示すフラグ (newprocess)

• 各バッチの処理時間。時間単位 (time)

• バッチの温度。摂氏 (temp)

• バッチで使用する化学薬品の供給業者 (A、B または C) を示すカテゴリカル変数 (supplier)

• バッチ内の欠陥数 (defects)

またデータに含まれる time_dev と temp_dev は、摂氏 20 度で 3 時間の標準プロセスから得られる時間と温度の絶対偏差をそれぞれ表します。

固定効果予測子として newprocess、time_dev、temp_dev および supplier を使用して一般化線形混合効果モデルを当てはめます。工場特有の変動に起因して品質に差がある可能性を考慮するために、factory 別にグループ化された切片の変量効果項を含めます。応答変数 defects はポアソン分布であり、このモデルの適切なリンク関数は対数です。係数の予測にラプラス近似メソッドを使用します。ダミー変数エンコードを 'effects' として指定すると、ダミー変数の係数の合計が 0 になります。

欠陥数はポアソン分布を使用してモデル化できます

$${\text{defects}}_{ij}\sim \text{Poisson}({\mu }_{ij})$$

これは一般化線形混合効果モデルに対応します

$$\log ({\mu }_{ij})={\beta }_0+{\beta }_1{\text{newprocess}}_{ij}+{\beta }_2{\text{time}\text{\_}\text{dev}}_{ij}+{\beta }_3{\text{temp}\text{\_}\text{dev}}_{ij}+{\beta }_4{\text{supplier}\text{\_}\text{C}}_{ij}+{\beta }_5{\text{supplier}\text{\_}\text{B}}_{ij}+b_i,$$

ここで

• $${\text{defects}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ で実行されたバッチで観測された欠陥数です。

• $${\mu }_{ij}$$ は、バッチ $$j$$ ($$j=1,2,...,5$$) 処理中の工場 $$i$$ ($$i=1,2,...,20$$) に対応する欠陥の平均数です。

• $${\text{newprocess}}_{ij}$$、$${\text{time}\text{\_}\text{dev}}_{ij}$$ および $${\text{temp}\text{\_}\text{dev}}_{ij}$$ は、バッチ $$j$$ 処理中の工場 $$i$$ に対応する各変数の測定値です。たとえば $${\text{newprocess}}_{ij}$$ は、工場 $$i$$ で実行されたバッチ $$j$$ 処理中に新プロセスが使用されたかどうかを示します。

• $${\text{supplier}\text{\_}\text{C}}_{ij}$$ および $${\text{supplier}\text{\_}\text{B}}_{ij}$$ はエフェクト (ゼロサム) コーディングを使用するダミー変数であり、バッチ $$j$$ 処理中に工場 $$i$$ で実行されたバッチに対して、それぞれ会社 C または B が加工化学薬品を供給したかどうかを示します。

• $$b_i\sim N(0,{\sigma }_b^2)$$ は、工場特有の品質変動に相当する、各工場 $$i$$ の変量効果の切片です。

glme = fitglme(mfr,'defects ~ 1 + newprocess + time_dev + temp_dev + supplier + (1|factory)',...
    'Distribution','Poisson','Link','log','FitMethod','Laplace','DummyVarCoding','effects');

変量効果係数の 99% の信頼区間を計算し、表示します。

[B,BNames,stats] = randomEffects(glme,'Alpha',0.01);
stats

stats = 
    Random effect coefficients: DFMethod = 'residual', Alpha = 0.01

    Group              Level         Name                   Estimate     SEPred     tStat       DF    pValue       Lower       Upper  
    {'factory'}        {'1' }        {'(Intercept)'}          0.29131    0.19163      1.5202    94      0.13182    -0.21251    0.79514
    {'factory'}        {'2' }        {'(Intercept)'}          0.15423    0.19216     0.80259    94      0.42423      -0.351    0.65946
    {'factory'}        {'3' }        {'(Intercept)'}         -0.26325    0.21249     -1.2389    94      0.21846    -0.82191    0.29541
    {'factory'}        {'4' }        {'(Intercept)'}         -0.42568    0.21667     -1.9646    94     0.052408    -0.99534    0.14398
    {'factory'}        {'5' }        {'(Intercept)'}           0.5453    0.17963      3.0356    94    0.0031051    0.073019     1.0176
    {'factory'}        {'6' }        {'(Intercept)'}         -0.10692    0.20133    -0.53105    94      0.59664    -0.63625    0.42241
    {'factory'}        {'7' }        {'(Intercept)'}          0.30404    0.18397      1.6527    94      0.10173    -0.17964    0.78771
    {'factory'}        {'8' }        {'(Intercept)'}         -0.16527    0.20505    -0.80597    94      0.42229    -0.70438    0.37385
    {'factory'}        {'9' }        {'(Intercept)'}         -0.14577      0.203    -0.71806    94       0.4745    -0.67949    0.38795
    {'factory'}        {'10'}        {'(Intercept)'}        -0.081632    0.20256      -0.403    94      0.68786    -0.61419    0.45093
    {'factory'}        {'11'}        {'(Intercept)'}         0.014529    0.21421    0.067826    94      0.94607    -0.54866    0.57772
    {'factory'}        {'12'}        {'(Intercept)'}          0.17706    0.20721     0.85446    94      0.39502    -0.36774    0.72185
    {'factory'}        {'13'}        {'(Intercept)'}          0.24872    0.20522       1.212    94      0.22857    -0.29083    0.78827
    {'factory'}        {'14'}        {'(Intercept)'}          0.21145    0.20678      1.0226    94      0.30913    -0.33221    0.75511
    {'factory'}        {'15'}        {'(Intercept)'}           0.2777    0.20345       1.365    94      0.17552    -0.25719    0.81259
    {'factory'}        {'16'}        {'(Intercept)'}         -0.25175    0.22568     -1.1156    94      0.26746    -0.84509    0.34158
    {'factory'}        {'17'}        {'(Intercept)'}         -0.13507    0.22301    -0.60568    94      0.54619     -0.7214    0.45125
    {'factory'}        {'18'}        {'(Intercept)'}          -0.1627    0.22269    -0.73061    94      0.46684    -0.74817    0.42278
    {'factory'}        {'19'}        {'(Intercept)'}         -0.32083    0.23294     -1.3773    94      0.17168    -0.93325    0.29159
    {'factory'}        {'20'}        {'(Intercept)'}         0.058418    0.21481     0.27195    94      0.78626    -0.50635    0.62319

stats の最初の 3 列には、グループ名、水準、変量効果係数の名前が含まれます。列 4 には変量効果係数の推定 EBP が含まれます。stats、Lower および Upper の最後の 2 列には、それぞれ、99% 信頼区間の下限と上限が含まれます。たとえば、factory の水準 3 における '(Intercept)' の係数について、推定 EBP は -0.26325 で、99% の信頼区間は [-0.82191,0.29541] です。

参照

[1] Booth, J.G., and J.P. Hobert. “Standard Errors of Prediction in Generalized Linear Mixed Models.” Journal of the American Statistical Association, Vol. 93, 1998, pp. 262–272.

参考

GeneralizedLinearMixedModel | coefCI | coefTest | fixedEffects